BFKL Spectrum of N=4 SYM: non-Zero Conformal Spin

We developed a general non-perturbative framework for the BFKL spectrum of planar N=4 SYM, based on the Quantum Spectral Curve (QSC). It allows one to study the spectrum in the whole generality, extending previously known methods to arbitrary values of conformal spin $n$. We show how to apply our approach to reproduce all known perturbative results for the Balitsky-Fadin-Kuraev-Lipatov (BFKL) Pomeron eigenvalue and get new predictions. In particular, we re-derived the Faddeev-Korchemsky Baxter equation for the Lipatov spin chain with non-zero conformal spin reproducing the corresponding BFKL kernel eigenvalue. We also get new non-perturbative analytic results for the Pomeron eigenvalue in the vicinity of $|n|=1,\;\Delta=0$ point and we obtained an explicit formula for the BFKL intercept function for arbitrary conformal spin up to the 3-loop order in the small coupling expansion and partial result at the 4-loop order. In addition, we implemented the numerical algorithm of arXiv:1504.06640 as an auxiliary file to this arXiv submission. From the numerical result we managed to deduce an analytic formula for the strong coupling expansion of the intercept function for arbitrary conformal spin.

This article is dedicated to the memory of Lev Nikolaevich Lipatov, who was a constant source of inspiration for us and who deeply influenced our research. He will be greatly missed.
1 Introduction N = 4 Super-Yang-Mills theory has been playing an important role in our understanding of Quantum Field Theories, especially in an AdS/CFT context. Due to the Kotikov-Lipatov maximal transcendentality principle [2,3] some of the results obtained in this theory can be directly exported to more realistic planar QCD. In this paper we describe how to efficiently perform calculations in this theory for one of the key QCD observables -BFKL spectrum, using integrability at any value of the 't Hooft coupling λ, which was discovered initially by Lipatov in the LO BFKL spectrum [4], and developed far beyond the perturbative regime in the N = 4 SYM in recent years. Lev Nikolaevich was one of the main driving forces behind this progress and it is deeply saddening for us to know that he left us in September 2017.
In the beginning we are going to briefly describe the meaning of the quantities studied in the present work in the context of high energy scattering. The total cross-section σ(s) for the highenergy scattering of two colorless particles A and B in the next-to-leading logarithmic approximation can be written as [5] σ(s) = d 2 qd 2 q (2π) 2 q 2 q 2 Φ A (q)Φ B (q ) a+i∞ a−i∞ dω 2πi mwhere Φ i (q i ) are the impact factors, G ω (q, q ) is the t-channel partial wave for the gluon-gluon scattering, s 0 = |q||q | and depend on the transverse momenta and s = 2p A p B , where p A and p B are the 4-momenta of the particles A and B respectively. For the t-channel partial wave there holds the Bethe-Salpeter equation ωG ω (q, q 1 ) = δ D−2 (q − q 1 ) + d D−2 q 2 K(q, q 2 )G ω (q 2 , q 1 ) , (1.2) where K(q, q 2 ) is called the BFKL integral kernel. It appears to be possible to classify the eigenvalues ω of this BFKL kernel using two quantum numbers: integer n (conformal spin) and real ν ω = ω(n, ν) . (1. 3) The function ω(n, ν) is called the Pomeron eigenvalue of the BFKL kernel or just the BFKL Pomeron eigenvalue and its values for different n and ν constitute the BFKL spectrum. For the phenomenological applications of the BFKL kernel eigenvalues with non-zero conformal spin see [6]. The object ω(n, ν) 1 in the planar N = 4 SYM will be studied in this work by means of integrability.
The study of integrable structures in 4d gauge theory has long and interesting history of development. Integrability in QCD and supersymmetric Yang-Mills theories appeared in two contexts. First, in the gauge theory, namely QCD, the Bartels-Kwiecinski-Praszalowicz (BKP) equation [7,8] for multi-reggeon states was reformulated by L.N. Lipatov [4] as the model with holomorphic and antiholomorphic hamiltonians, which has a set of mutually commuting operators originating from the monodromy matrix satisfying the Yang-Baxter equation. After that L.D. Faddeev and G.P. Korchemsky in [9] proved this model to be completely integrable and equivalent to the spectral problem for SL(2, C) XXX Heisenberg spin chain. Then, in the context of high-energy scattering there was considered a certain class of light-cone operators in QCD and supersymmetric Yang-Mills theories and in [10][11][12][13] the problem of finding the anomalous dimensions of light-cone operators was formulated in terms of SL(2, R) Heisenberg spin chain.
The other achievement was that the maximally supersymmetric N = 4 Yang-Mills theory in 4 dimensions, which is dual to AdS 5 × S 5 type IIB superstring theory was shown to be integrable [14,15]. The study of the integrability structure of the latter theory allowed to explore its spectrum in the non-perturbative regime. The solution to the spectral problem was formulated in terms of the Quantum Spectral Curve (QSC) [16,17] (for the recent reviews see [18] and [19]). Nevertheless, until recently it was not known how to build the bridge between the integrability in the BFKL limit and integrability found in the AdS/CFT framework. In [20] the 4-loop Asymptotic Bethe Ansatz (ABA) contribution to the anomalous dimension of the twist-2 sl(2) operators was analytically continued to the non-integer spins and compared with the corresponding prediction from the BFKL Pomeron eigenvalues. This analytic continuation to non-integer spins was incorporated into the QSC formalism in [21] for twist-2 operators from the sl(2) sector and then in [22] the Faddeev-Korchemsky Baxter equation [9] for Lipatov SL(2, C) spin chain was derived correctly reproducing the leading order (LO) BFKL Pomeron eigenvalue. In addition, QSC allowed to calculate analytically [23] the previously unknown next-to-next-to-leading order (NNLO) BFKL eigenvalue in the N = 4 supersymmetric Yang-Mills theory. At the same time, a very efficient numerical algorithm was constructed in [1], which allows to study not only the BFKL limit of the spectrum of the theory, but the whole anomalous dimension of a given operator for arbitrary values of the charges.
In [22] the twist-2 sl(2) operators of the form O = trZD S + Z + (permutations) , (1.4) were considered and from the perturbative calculations in the gauge theory for the case of even integer S we know the dimension of these operators ∆ as a function of S up to several loops order.
In the QSC framework the solution of the Baxter equation for the spectrum of such operators in the case of zero conformal spin n and integer even spins S was obtained in [16]. Then in [24,25] there was found the solution of this Baxter equations valid for arbitrary spin S, which leads to the anomalous dimension of the twist-2 sl(2) operators analytically continued for non-integer spin S. After making this analytic continuation in the BFKL regime we are able to exchange the roles of ∆ and S obtaining S + 1 = ω(n = 0, ν), where ν = −i∆/2 and ∆ is the dimension of the operator in question.
In the present work we consider the generalization allowing for an arbitrary value of the conformal spin. Namely, we consider the length-2 operators O = trZD S1 + ∂ S2 ⊥ Z + (permutations) . (1.5) For the operators (1.5) we follow the same strategy as for the case of zero conformal spin. Analogously to that case we build the analytic continuation in the spins S 1 and S 2 , which are identified with the spin S and conformal spin n respectively. Let us illustrate this analytic continuation with the Figure 1. The physical operators, for which the sum of non-negative integer S = S 1 and n = S 2 is even, are depicted with the dots. Then, flipping the roles of the dimension ∆ and S = S 1 we can reach the BFKL regime described by the quantity ω(n = S 2 , ν) = S 1 + 1, where ν = −i∆/2. The way to proceed with the problem in question is to first generalize the QSC approach to non-integer values of S 1 (as was already done in [21]) and then also to non-integer values of S -1 |n|+1 -|n|-1 -|n|-2 Figure 1. Trajectory of the length-2 operator for conformal spin n = S2 as a function of the full dimension ∆. The dots correspond to local operators trZD S + ∂ n ⊥ Z. For the local operators S + n is restricted to be even. S 2 . We describe the technical details of this procedure in the Section 2. This allows to treat ω(n, ν) as an analytic function of both its parameters which simplifies both analytical and numerical considerations. This gives a universal framework for studying the BFKL spectrum in full generality for all values of the parameters on equal footing within the extended QSC formalism.
Having formulated the problem as an extension of the initial QSC, a number of methods, initially developed for the local operators, became available for the BFKL problem. In particular we are enabled to employ a very powerful numerical algorithm [1] after some modifications. As we take the spins S 1 to be continuous variable we can consider instead of the function ∆(S 1 , S 2 ) the function S 1 (∆, S 2 ). Then, using the algorithm we build the operator trajectories for different values of conformal spin S 2 and the dependencies of the spin S 1 on the coupling constant g for different values of conformal spin and dimension ∆ (including a particular interesting intercept function corresponding to ∆ = 0). Having the numerical results for the operator trajectories we were able to fit the numerical values of the BFKL kernel eigenvalues 2 , which were confirmed in [27] using a different method.
Another method available within the QSC formalism is an efficient perturbative expansion developed in [23,[28][29][30][31]. We applied this method to find the value of the Pomeron intercept for the arbitrary value of the conformal spin up to 3 loops. Our result is in full agreement with [27] at the NNLO level, but we also give a prediction for the next NNNLO order.
Then, we found and studied in detail a particularly interesting point in the space of parameters of the BFKL Pomeron eigenvalue. This is the "BPS" point ∆ = 0 and n = 1. As we have confirmed both numerically and analytically, the operator trajectory goes through the point S = −1, n = 1 and ∆ = 0 for any value of the coupling constant g. Studying the vicinity of this point we were able to find two non-perturbative quantities: "slope-to-intercept function" and "curvature function". The first function is the first derivative of S (∆, n) with respect to n at the point ∆ = 0, n = 1 and the second function is the second derivative of S (∆, n) with respect to ∆ at the same point. We used the methods developed in [21] to compute analytically these quantities non-perturbatively to all orders in g.
Finally, we were able to identify the intercept function in the strong coupling expansion up to the 4th order. To obtain it we utilized the dependencies of the intercept on the coupling constant calculated by the QSC numerical method. By conducting the numerical fit of these dependencies for different values of conformal spin n we predict the formula for the intercept strong coupling expansion up to the 4th order for arbitrary conformal spin.
Let us present a brief summary of the quantities we calculated. They include the NNLO intercept function (4.16) and the non-rational part of the NNNLO intercept function (4.18). The other quantities we computed exactly to all orders in the 't Hooft coupling constant are the slopeto-intercept (5.43) and the curvature (5.131) functions with the strong coupling expansions of these functions given by (5.49) and (5.136) respectively. In addition, there was written the strong coupling expansion (6.2) of the intercept function for arbitrary conformal spin n. We also implemented the numerical method for finding the eigenvalues at arbitrary values of the parameters in Mathematica, the corresponding files code_for_arxiv.nb and BFKLdata.mx can be found in the attachments to this arXiv submission. See description.txt file for the description.
This work is organized as follows. In the Section 2 we give the general introduction into the QSC approach, extending it to the situation when both spins are non-integer. The Section 3 describes our numerical results. The Section 4 contains the weak coupling analysis. In the Section 5 we analyze the expansion near the BPS point to find the non-perturbative quantities such as slope of intercept and curvature functions. In Section 6 we analyze the Pomeron intercept at strong coupling.

Description of the QSC based framework
In this Section we are going to present the framework which we use to solve the QSC [16,17] and whose derivation is based on the analytic and asymptotic properties of the Q-functions. First, we reformulate the QSC in terms of gluing matrix. Namely, we start from the several axioms concerning the analytic structure of the Q-system and the symmetries which preserve the QQrelations and derive from them the so-called gluing conditions. These gluing conditions already appeared in [17,23] but our approach presented below does not utilize the notion of µ-and ωfunctions to obtain the gluing matrix. Second, using the connection between the asymptotics of the certain subset of Q-functions and the global charges together with their analytic properties, the system of constraints for the gluing matrix is derived. It appears to be possible to solve these equations in some physically interesting cases. Namely, we find the gluing matrix for the case when both AdS spins S 1 and S 2 are integers of the same parity and its form appears to be very simple and in complete agreement with the result of [17]. Then we consider a more general case of non-integer AdS spins S 1 and S 2 , which is particularly interesting for the exploration of the BFKL regime. For this case we have not found the general solution for the gluing matrix, however we found the certain subclass of solutions and it appears to be applicable to our quantities of interest. We mostly follow the original paper [17], but the discussion of the gluing matrix and the extension to the non-integer quantum numbers is new. The reader familiar with the QSC formalism could skip to Subsection 2.5.

Algebraic part of the construction
QSC consists of a set of Q-functions of the complex spectral parameter u and relations between them. We will restrict ourselves to the most essential parts of the construction but still keeping the discussion self-contained. For more detailed description of the QSC see [16,21] and for the pedagogical introduction see [17].
In total there are 256 Q-functions Q a1,...,an|i1,...,im (u) totally antisymmetric in the two groups of "bosonic" (a's) and "fermionic" (i's) indices with 1 ≤ n, m ≤ 4, however not all of them are independent. The main building blocks of the QSC construction are the 4 + 4 "elementary" Qfunctions: Q a|∅ (u), where a = 1, . . . , 4, and Q ∅|i (u), where i = 1, . . . , 4. Setting the normalization Q ∅|∅ = 1 and starting from these 8 Q-functions, one can recover the whole Q-system applying the QQ-relations written in [17]. In particular the QQ-relation for the Q-function with one "bosonic" and one "fermionic" index looks as follows and Q a|i (u) is a solution of (2.1). From now on we are going to use the shorthand notation for the shift in the variable u: In a similar way one can build all 256 Q-functions out of the basic 8 mentioned above. One should also impose the quantum unimodularity condition The Q-function Q a|i allows to write the Q-functions with one upper index in a concise form From the condition (2.2) it can be shown that In addition, the Q-system has a symmetry, which is called the H-symmetry [17] and which leaves the QQ-relations intact. It corresponds to the transformations of Q-functions by i-periodic matrices that rotate the "bosonic" and "fermionic" indices separately. Their form for all Q-functions can be found in [17], but in this Section we need the explicit form of them only for the Q-functions with one index. They are where H B (u) and H F (u) are i-periodic 4 × 4 matrices. The determinants of these matrices have to for the quantum unimodularity condition (2.2) not to change under such H-rotations. The important particular case of this symmetry is the rescaling of the Q-functions with one index. It acts as follows The equation (2.1) allows to obtain a 4th order Baxter equation for the functions Q ∅|i (u), i =, 1 . . . , 4. In [22] this equation was derived and looks as follows It is also possible to show from the same equation as (2.1) for the Q-functions with upper indices that the functions Q ∅|i (u), i = 1, . . . , 4 satisfy the 4th order Baxter equation, which looks as (2.11) but with Q a|∅ exchanged with Q a|∅ . For the sake of conciseness we do not write this Baxter equation explicitly.
After finishing the description of the algebraic structure of the QSC essential for the formulation of the QSC equations in the next Subsection we describe the analyticity properties of the Qfunctions, which constitute the crucial part of our QSC framework.

Analytic part of the construction
To describe the analytic structure of the Q-system we have to first define the analyticity properties of the basic set of Q-functions: Q a|∅ (u), a = 1, . . . , 4 and Q ∅|i (u), i = 1, . . . , 4. The only singularities of these functions are quadratic branch points which come in pairs at the positions ±2g + ik, where k ∈ Z. For each pair of branch points we can choose either short cut on the interval [−2g + ik, 2g + ik] or a long cut (−∞ + ik, −2g + ik] ∪ [2g + ik, ik + ∞), where k is some integer. In what follows the sheet of the Q-function with only the short cuts is called physical and the function on this sheet is denoted byQ(u), while on the sheet, where all the cuts are long, is called mirror and the function is designated byQ(u) on it. The continuation of any function f (u) under the cut on the real axis is denoted byf (u). The branch points of all the functions we will consider are quadratic, i.e.f (u) = f (u).
In what follows we will denote the functions Q a|∅ (u), Q a|∅ (u), Q ∅|i (u) and Q ∅|i (u), with prescribed analytical properties, as P a (u), P a (u), Q i (u) and Q i (u) respectively. To proceed let us write the asymptotics of the Q-functions with one index. We know that all the Q-functions including P a , P a , Q i and Q i have the power-like asymptotics at large u, which for the basic 8 Q-functions can be taken from [17] (2.14) As we know from the classical integrability of the dual superstring σ-model (see, for example, [18]), the Pand Q-functions at least have the quadratic branch points at u = ±2g. From the asymptotics (2.13) and (2.14) we can expect the cuts of P-functions to be short 3 . The minimal choice for the functions P a (u) and P a (u), a = 1, . . . , 4 is to have only one short cut on the real axis. From the asymptotics of the Q-functions we can see that they have a nontrivial monodromy around infinity, thus we have to assume the cuts of these functions to be long. So again the minimal choice for Q i (u) and Q i (u), i = 1, . . . , 4 would be to have only one long cut on the real axis. The analytic structure of Pand Q-functions is illustrated on the Figure 2. Notice, that because the functions Q i and Q i have the long cuts in the complex plane, their asymptotics prescribed from the large u limit of the superstring σ-model hold in the upper half-plane and in the lower half-plane they can be different, therefore the third and fourth formulas from (2.13) are valid for Im u > 0 4 . It is convenient for us to introduce some short-hand notations as in [17]: UHP -upper halfplane, LHP -lower half-plane, UHPA -upper half-plane analytic and LHPA -lower half-plane analytic. Besides, we are going to use LHS and RHS for left and right-hand side respectively.
As we have introduced the analytic structure of the basic set of Q-functions, let us proceed with the consideration of the other ones. We define the function Q a|i as an UHPA solution of the equation with the asymptotic In what follows we will denote the UHPA Q-functions of the Q-system obtained from P a , Q i and Q a|i by the application of the QQ-relations by curly Q (as it was done in [17] to underline that these Q-functions have certain analytic properties and where the corresponding Q-system was called fundamental). For the Hodge-dual Q-functions, which are UHPA as well and satisfy the same QQ-relations, we will also use curly Q to depict its Q-functions. Substitution of (2.13) and (2.16) into Q i = −Q + a|i P a and P a = −Q + a|i Q i expressed from (2.5) and (2.6) themselves leads us to the systems of equations for A a A a and B i B i respectively, which can be solved [17,18] and give the result where there is no summation over the indices a 0 and i 0 implied. For further convenience we introduce the shorthand notations Each function Q a|i (u) is analytic for Im u > −1/2 due to the fact that both P a (u) and Q i (u) are UHPA. As it was mentioned in the Subsection 2.1 with the usage of the QQ-relations we can restore the remaining Q-functions thus building the UHPA Q-system. The Hodge-duality (2.3) does not change the analytic properties, therefore the Hodge-dual Q-system with the upper indices has the same analytic properties, i.e. UHPA. Now we are going to turn to the analytic structure of the functions Q i (u), i = 1, . . . , 4. Let us remember the formulas (2.5). One can notice that the functions Q i and Q i with −Q + a|i P a and Q a|i+ P a respectively coincide only in the UHP, because their analytic structure in the LHP is different. Indeed, we can see that if we rewrite the QQ-relations for Q a|i using (2.5) it is possible to find the values of Q a|i for Im u < −1/2. In the strip −k − 1 < Im u < −k for k = 0, 1, 2, . . . the functions Q − a|i are given by From (2.19) we can see that the functions Q a|i have the infinite number of short cuts at the horizontal lines with Im u = −k − 1/2 for k = 0, 1, 2, . . . and due to (2.4) Q a|i has the same structure of cuts in the complex plane. Therefore, the functions −Q + a|i P a and Q a|i+ P a have also the infinite ladder of cuts at Im u = −k for k = 0, 1, 2, . . ., which clearly do not coincide with the analytic structure of Q i and Q i in the LHP, who are LHPA.
However, we can resolve this difficulty by interpreting the analytic continuation of −Q + a|i P a and Q a|i+ P a under their short cut on the real axis from above as Q i and Q i respectively. To formulate this clearly let us use the hats and checks introduced in the beginning of the present Subsection, which denote the values of the Pand Q-functions on the different sheets. In these notations first of all the equations (2.15) and (2.19) determineQ a|i on the physical sheet with the short cuts. Then, the values of theQ i andQ i on their physical sheet with the short cuts are given bŷ (2.20) and coincide with the Q-functions on the mirror sheet with the long cut on the real axis in the UHPQ Whereas in the LHP we interpret the analytic continuation ofQ i andQ i under their cut on the real axis asQ i andQ i in the LHP where tilde denotes the analytic continuation under the cut on the real axis. Looking at the obtained picture from above, we conclude that there is no fundamental reason to choose the generated Q-system to be UHPA. Indeed, there exists a transformation of complex conjugation, which preserves the QQ-relations but interchanges UHPA with LHPA. Its explicit form is written in [17] Q a1,...,an|i1,...,im (u) → (−1) (m+n)(m+n−1) 2Q a1,...,an|i1,...,im (u) . (2.23) The transformation (2.23) generates the Q-system which is LHPA and satisfies the same QQrelations as the initial UHPA Q-system. It should be noted that the Hodge-dual Q-system also admits such a transformation Q a1,...,an|i1,...,im (u) → (−1) (m+n)(m+n−1) 2Q a1,...,an|i1,...,im (u) .

(2.24)
Let us now remember the analyticity properties of the function Q i from (2.21) and (2.22).
We see that the functionsQ i are LHPA and therefore have the same analyticity properties as the functionsQ i andQ i . From the strong coupling limit of the superstring σ-model and its classical integrability (see the pedagogical explanation of this in [18]) we know that for the case of integer spins S 1 and S 2 each functionQ i , i = 1, . . . , 4 coincides with the certain function from the setQ j , j = 1, . . . , 4 in this limit. Thus, summarizing all this, we impose the equality of the LHPA functionŝ From now on we will call (2.25) the gluing condition and M ij (u) the gluing matrix, whose properties we will analyze below. We formulated QSC in the form (2.25) because this form is convenient to analytically continue the QSC solution for the case of non-integer spins S 1 and S 2 . The transformation (2.23) generates the Q-system which is LHPA and satisfies the same QQ-relations as the initial UHPA Q-system. As there is no principal difference between UHPA and LHPA Q-systems and they describe the same spectral problem and due to the unitarity of the N = 4 SYM theory, the UHPA and LHPA Q-systems have to be related by the symmetries of the Q-system, namely, the combination of the Hodge duality and H-symmetry 5 . Thus, we can interpret the gluing matrix M ij as an i-periodic matrix of the H-transformation, which, in particular, relates the Q-functions with one lower and one upper "fermionic" index on the mirror sheeť where −t means that the inverse matrix is transposed.
Using the analyticity properties of the functionsQ i andQ j we are able to establish some properties of the matrix M ij . Utilizing the i-periodicity of the matrix M ij , it is possible to express its elements in terms ofQ i andQ j . From the i-periodicity of M ij (u), (2.26) and remembering the QQ-relation for Q-function with 4 "fermionic" indices Q ∅|1234 = det First of all let us show that the matrix elements M ij do not have any branch points. To see this let us notice that for Im u < −3/2 the Q-function on the LHS and the determinant on the RHS of (2.27) are analytic and do not have branch points. Therefore in the same region and due to the i-periodicity M ij has to be free from the branch points in the whole complex plane.
Second, in principle, M ij can have poles. As M ij is i-periodic, the existence of a pole, for example, in the point u 0 automatically leads to the infinite number of poles in the points u 0 + ik, where k ∈ Z. However, at least in the region Im u < −3/2 the RHS of (2.27) is analytic, thus the poles of M ij have to be compensated by the zeroes ofQ ∅|1234 in the same points or, in other words, there exists k 0 such thatQ ∅|1234 (u 0 + ik) = 0 for k ≤ k 0 . In its turn this means that the number of zeroes ofQ ∅|1234 is infinite and these zeroes accumulate at infinity. We know thatQ ∅|1234 has power-like asymptotic, then there exist such a and b thať (2.28) 5 The presence of the Hodge duality can explained from the consideration of the classical limit of the superstring σ-model, which says that the analytic continuation of the Q-functions with the lower indicesQ i (u) are related to the Q-functions with the upper indicesQ i , not the lower ones.
However, evaluating the LHS of (2.28) at u = u 0 + ik for k ≤ k 0 leads to a contradiction with the RHS of (2.28). From this contradiction we conclude that M ij cannot have any poles.
Summarizing what was said we see that the gluing matrix M ij is analytic in the whole complex plain. For the physical state, which means that the spins S 1 and S 2 are integer and the asymptotics of the functionsQ i are the power-like, analyticity and i-periodicity of the gluing matrix M ij leads us to the conclusion that it is constant in this case. Now let us return back to the equations (2.26). As we know that the matrix M ij is free of any singularities, then analytically continuing both sides of the equations (2.26) to the sheet with the short cuts we obtainQ Analytic properties of the Q-functions allow us to establish one important property of the gluing matrix. Analytically continuing both sides of the first equation from (2.29), using the fact that due to the quadratic nature of the branch points analytic continuation and complex conjugation commute with each other, and then applying the second equation from (2.29) we derivê (2.30) Noticing that (2.30) is true for any point u and applying the same trick as above and using the analyticity properties of the gluing matrix, we arrive to the conclusion that and the gluing matrix is hermitianM In what follows we are going mainly to deal with the Q-functions on the physical sheet therefore from now on we omit the hats and checks above the designations of the Q-functions implying that all the Q-functions are considered on the sheet with the short cuts if the opposite is not mentioned specifically. Now we are ready to find the other constraints on the gluing matrix which follow from the conjugation and parity symmetries of the Q-system. To do this we will need the 4th order Baxter equation for the functions Q i (u) to see if the certain properties of the P a (u) and P a (u) can allow us to relate Q i (u),Q i (u) and Q i (−u). In the two subsequent Subsections we analyze the implications of the conjugation and parity properties respectively.

Complex conjugation symmetry
Let us now concentrate on the conjugation properties of the Pand Q-functions assuming the charges ∆, S 1 and S 2 to be real. In [17] from the reality of the energy, Y-and T-functions and the fact that the complex conjugation supplemented by the certain sign factor is the symmetry of the Q-system it is shown that complex conjugation is equivalent to some H-symmetry transformation already mentioned in the Subsection 2.1. The matrix h B of this transformation was proven in [17] to be constant due to analytic properties and power-like asymptotic of the P-functions. Then there was found a transformation which allows to make all P-functions with lower indices purely real and thus the P-functions with upper indices pure imaginary. However, in our calculations we use the different normalization and make the other H-rotation by multiplying P 3 and P 4 by i and thus P 1 and P 2 also by i and obtain P 1,2 = P 1,2 ,P 3,4 = −P 3,4 ,P 1,2 = −P 1,2 ,P 3,4 = P 3,4 or, in other wordsP Given the conjugation properties (2.33), we see that the Baxter equation (2.11), written for the Q-functions with prescribed analytic properties with D k andD k given by (2.12) with Q a|∅ = P a , remains the same, but forQ i (u) now. Thus, the functionsQ i satisfy the equation (2.34) too. As the functionsQ i (u) constitute a basis in the space of solutions of the Baxter equation (2.34) this means that there has to exist an i-periodic matrix As in [32] from (2.33) and (2.5) this matrix can be found to be It is i-periodic Ω j++ i = Ω j i (see Appendix A for the proof) and using this it is not hard to show that which means that Ω −1 =Ω and The matrixΩ j i also relatesQ j and Q iQ and vice versa To determine the consequences of the conjugation symmetry for the gluing matrix we substitute (2.35) into the first gluing condition from (2.29) and obtaiñ Let us analyze the matrix M ij Ω k j more closely. As it is a product of two i-periodic matrices it has also to be i-periodic. We remember that according to its definition (2.36) the matrix Ω k j (u) has an infinite ladder of short cuts. Using the result of [32] we get the discard of Ω k j (u) Multiplying both sides of (2.42) by M ij and utilizing the first gluing condition from (2.29) we derive the equation We note that the RHS of (2.43) is antisymmetric in the indices i and k, thus we conclude that the function M ij Ω k j + M kj Ω i j has no cuts on the real axis. As this function is i-periodic it follows that M ij Ω k j + M kj Ω i j is analytic in the whole complex plane. Let us introduce a new notation where ω ik are the i-periodic functions on the sheet with the short cuts. Remembering (2.41) and applying (2.44) on the sheet with the short cuts we obtain the equatioñ On the other hand, the functions ω ik are i-periodic on the sheet with the short cuts, thus on the sheet with the long cuts their analytic continuation under the cut on the real axis is given by the simple formulaω ik = ω ik ++ . If we return several steps back, we can derive from the QQ-relations the equation for the function Q ab|ij , which looks almost exactly like (2.47) Recalling the notion of µ-functions introduced in the QSC framework in [16,17], which are i-periodic on the sheet with the long cuts, we multiply both sides byμ ab , which leads us to are antisymmetric in i and j due to the antisymmetry of the Q-function Q ab|ij , it is natural to impose the constraint that ω ik is also antisymmetric and In the following Subsections we are going to exploit (2.50) to constrain the gluing matrix for different spins S 1 and S 2 .

Parity symmetry
Now we are going to describe the parity properties of the Q-system. For a large class of states the P-functions possess the certain parity. Such states include the states with the charges J 1 = 2, J 2 = J 3 = 0, which we consider in the rest of the paper and also the ground state with the charges J 1 = 3 and J 2 = J 3 = 0, which is relevant for the BFKL Odderon eigenvalue (see [33,34]). Thus, for the case J 1 = 2, J 2 = J 3 = 0 we havẽ As we understood the analytic structure of Pand Q-functions, taking into account the asymptotics of these functions expressed in terms of the charges (2.51), it is natural to assume the existence of the certain symmetry between the Q-functions with lower and upper indices, which also changes sign of ∆. This symmetry takes a particularly simple form for the following choice of the normalization of the P- and Q-functions which can be set with the usage of the rescaling symmetry (2.10). We obtain 7 and χ ab is the same matrix for the left-right symmetric states as in [16,17]. For the operators we examine (J 1 = 2 and J 2 = J 3 = 0) in the following Sections the Pfunctions have the certain parity. Their parity is dictated by the asymptotics of the P-functions The symmetry (2.56) is a symmetry of the Baxter equation (2.34), thus Q i (−u) is also a solution of (2.34). Utilizing the same logic as in the case of the complex conjugation, we conclude that there exists an i-periodic matrix Θ j i (u) (see the proof in Appendix A) such that It is possible also to find the matrix with such a property. Utilizing again (2.5), we obtain where the summation over a is implied. The matrix Θ has the property Θ j i (u)Θ k j (−u) = δ k i and thus Θ −1 (u) = Θ(−u). We can write The matrix Θ j i (−u) also relates Q j (−u) and Q i (u) Analogously to the consideration of parity symmetry, we can find the discard of Θ k j on the cut which is situated on the real axis Using the matrix Θ j i (u) from (2.58), the gluing conditions (2.41) and (2.59) we are able to introduce another gluing matrix The change of the sign ∆ → −∆ is a symmetry of the equation and it should map one solution to another solution. One can check that P a = χ ab P b , P a = χ −1 is also a solution to the QQ-relations and the gluing conditions but with ∆ flipped to −∆. As it can be seen explicitly in the Appendix B in the weak coupling limit these two solutions coincide therefore our solution is mapped onto itself. Given the starting point the recursive procedure described in the Section 4 is non-ambiguous we conclude that this property holds to all orders in the coupling constant. and as the Q-functions on the both sides of the gluing conditions are LHPA, then by the same arguments as for M ij (u) the matrix L ij (u) is also analytic in the whole complex plane. Analogously to the case of the gluing matrix M ij (u) by going under the cut twice we derive the property Together with (2.32), (2.50) and (2.64) condition (2.65) constitutes the set of equations which are used to calculate the gluing matrix for different values of the spins S 1 and S 2 .
Since for the states in question the P-functions have the certain parity, this has some consequences for the asymptotic expansion of the Q-functions. As it is explained in detail in the Section 3 with the description of the numerical algorithm the certain parity of the P-functions leads to the form (3.6) of the asymptotic expansion of Q a|i By applying analogous arguments to the QQ-relation for the Hodge dual function Q a|i , we conclude that the asymptotic expansion of Q a|i is also given by Then we remember that Q i = −Q + a|i P a and Q i = Q a|i+ P a , which after the substitution of (3.1), (2.66) and (2.67) lead us to the asymptotic expansions at infinity of the Q-functions In what follows we will regard to the Q-functions with the asymptotic expansions (2.68) as having the "pure" asymptotic, as these expansions contain the powers of u, which differ only by an integer number. The asymptotic expansions (2.68) will be important in determining the structure of the matrices Ω j i (u) and Θ j i (u), which are analyzed below.

Constraining the gluing matrix
In the present Subsection we are going to derive the set of equations for the elements of the gluing matrix originating from the conditions found in the previous Subsections. To remind briefly the QSC framework we are using let us recall the constraints on the gluing matrices known by now. The non-degenerate matrices M ij (u) and L ij (u) satisfy the following set of constraints Now we are able to consider the gluing matrix for the case of different AdS spins S 1 and S 2 solving the set of constraints (2.69).

Integer S 1 and S 2
Let us start our study from the situation when all the charges except for the dimension ∆ are integer. More precisely, in the present Subsection we address the case when the spins S 1 and S 2 have the same parity. This is motivated by the fact that for S 2 = 0 the physical states have even non-negative S 1 8 . To analyze the constraints (2.69) more closely we need to find the properties of the matrices Ω j i (u) and Θ j i (u). In what follows we will need the asymptotics of the matrices Ω j i and Θ j i (u). To analyze them let us remember (2.35) and (2.57). As the asymptotics of the Q-functions on the both sides of (2.35) and (2.57) are power-like and the Ω-matrix consists of the i-periodic functions, the series expansions of Ω j i (u) and Θ j i (u) for |Re u| 1 are given by where the signs correspond to expansion at +∞ and −∞ respectively. It should be noted that (2.70) does not have any growing terms on the RHS, because this would violate the power-like asymptotic of the Q-functions.
The asymptotics of the Q-functions are pure and the asymptotic expansion is given by (2.68). Looking at the values ofM i , i = 1, . . . , 4 one may think that becauseM 1 −M 2 andM 3 −M 4 are integers, there could potentially appear a mixing of Q 1 (u) with Q 2 (u) and Q 3 (u) with Q 4 (u) as this does not violate the purity of the asymptotics. However as the spins S 1 and S 2 have the same parity andM the functions Q 1 , Q 2 and Q 3 , Q 4 cannot mix, because their asymptotic expansions (2.68) contain only even powers of u in the round brackets. Therefore the matrices (Ω First, we consider the matrix Ω j i (u) and remember (2.35). If Re u tends to +∞, then, as But if Re u tends to −∞ the situation is a little more subtle. The functions Q i (u) have an infinite ladder of short cuts going down from the real axis, while the functionsQ i (u) have the same ladder of cuts going up. Then taking the limit of Re u to −∞ we have to go to −∞ along the semicircle in the UHP for Q i (u), i.e. Q i (u) B i e iπ(Mi−1) (−u)M i−1 and along the semicircle in the LHP forQ i (u), i.e.
To sum up, we obtain Second, analyzing the matrix Θ j i (u) from (2.57) is analogous. Thus, applying the arguments from the previous paragraph, we see that at Re u tending to +∞ we have to go around the semicircle in the UHP and As it was explained for example in [18] in the strong coupling limit the asymptotics of the functionsQ i (u) are some powers of u, then the only possible ansatz for the gluing matrix M ij (u) is to assume it to be a constant matrix. Thus for the case in question we obtain from (2.50) and (2.72) rather simple conditions Combining the two conditions (2.74) and (2.75) we obtain the following Let us see now which additional restrictions do we have in the case J 2 = J 3 = 0 and J 1 = 2. First of all from our assumptions about the asymptotics of the functions Q i (u) we understand that the gluing matrix L ij given by (2.64) has to be constant, i.e. L ij (u) = L ij is a symmetric matrix. Therefore, from the (2.64), (2.72) and (2.73) we immediately find (2.77) Then, using (2.74) and the symmetry of L ij , we derive It is easy to see that if (2.78) is true then (2.76) is also true. We have to calculate the differences between the chargesM i to determine which elements of the matrix M ij are non-vanishing. It appears that Thus for the case of integer spins S 1 and S 2 we are left with the spins S 1 and S 2 with the same parity, which is consistent with our initial setup. Therefore, only the matrix elements M 12 =M 21 and M 34 =M 43 are non-zero. Then, in the case of integer spins S 1 and S 2 of the same parity we obtain the following gluing matrix (2.80) Using also (2.74) and (2.75), which are equivalent for S 1 of S 2 of the same parity asM 1 −M 2 and M 3 −M 4 equal 1 modulo 2, we are able to fix the phases of the non-zero matrix elements of (2.80) Now let us start the consideration of the case when at least one of the spins is not integer as this is particularly interesting for the BFKL limit.

Non-integer S 1 and S 2
First of all, from the asymptotics (2.51) we immediately see that if at least one of the charges S 1 or S 2 is non-integer, then not to violate the purity of the asymptotic expansions (2.68) the matrices (Ω ± ) j i cannot mix different Q-functions and have to be diagonal. Therefore, these matrices are given by (2.70). This means, that in the case of at least one non-integer spin under the assumption that the gluing matrix is constant we obtain the same constraint (2.78). However, as S 1 or S 2 or both spins are non-integer, all the differencesM i −M j are non-integer in general, therefore we conclude that M ij = 0. Then we have to modify the ansatz for M ij (u).
The matrix M ij (u) is analytic and i-periodic, so the minimal choice would be to add the terms proportional to e 2πu and e −2πu and this is consistent with what we know from the consideration of the BFKL limit for which S 1 approaches −1 and S 2 = 0 (see [22]). From the previous conditions (2.32) it follows that the matrices M ij 1,2,3 are hermitian. Substituting (2.82) into (2.50) we obtain the following conditions for the matrix M ij (u) Let us remember that for the case in question J 2 = J 3 = 0 and J 1 = 2. The matrix L ij (u) is given by the formula (2.64). Taking the limits u → ±∞ and remembering the expansions (2.72) and (2.73) we have to assume the existence of the exponential contributions to L ij (u) where the matrix L ij 1 is symmetric and L ji 2 = L ij 3 due to (2.65) and the latter two of them are given by Exploiting the symmetry L ji 3 = L ij 2 and the relation from (2.83) we derive As we observed in the case of integer spins S 1 and S 2 the determinant of the gluing matrix is constant. According to (2.82) in the case of at least one non-integer spin this determinant is not guaranteed to be integer. However, if we assume for a moment that the determinant of (2.82) contains exponents, the form of the second gluing condition from (2.29) will contain exponents in the denominator. But as there is no preference to upper and lower indices, which get exchanged under ∆ → −∆ symmetry, we have to assume that both gluing conditions (2.29) include the exponents e 2πu only in the numerator of M ij (u), therefore we impose a new constraint Let us now show that staring from some simple ansatz for the gluing matrix M ij (u) we are able to solve the constraints (2.88) and (2.89). We saw from the implementation of the numerical algorithm described in the Section 3 that in the case when both spins S 1 and S 2 are non-integer it is sufficient for convergence of the numerical procedure to allow presence of exponents in M 13 (u) and M 14 (u) only. Thus, we have the following ansatz for the hermitian matrices M ij where the non-zero matrix elements are subject to the relations (2.83) and (2.88). From the constraint (2.50) we find the equation Application of the constraint (2.89), i.e. the demand of the absence of the powers of e 2πu in the determinant of the gluing matrix leads us to the following equations As the spins S 1 and S 2 are non-integer and we assume their difference to be non-integer too and from ( where the elements of the matrices M ij 2,3 are subject to (2.83) and (2.88). As we have the relation (2.88) it is sufficient to write only the phases of the non-zero matrix elements of the matrix M ij 2 extracted from (2.83) Let us point out that the construction presented above will provide an analytic continuation to all values of S 2 from the integer values S 2 ≥ 0. However, this analytic continuation breaks down the symmetry S 2 → −S 2 , which is naively present in the QSC, as one can see from the asymptotic (2.14). The analytic continuation, which describes perfectly positive integer S 2 will produce poles at negative integer S 2 . This could look a bit puzzling, but the resolution of this paradox is in the existence of the second solution for the mixing matrix which is obtained by relabeling indices in accordance with S 2 → −S 2 . In practice result must be even in S 2 and it is enough to consider S 2 ≥ 0 so it is sufficient to use the mixing matrix presented above.
To sum up the contents of the present Section, we have to point out several things. First, we formulated the algebraic structure of the Q-system by writing down the QQ-relations and 4th order Baxter equation originating from them. Second, the analytic structure of the Q-system was motivated from the solution of the classically integrable dual superstring σ-model and the QQrelations. The symmetries of the Q-system allowed us to introduce the gluing conditions for which we managed to impose several constraints. These constraints were partially solved for different values of the spins S 1 and S 2 . We examined the case when both spins are integer and non-integer. In the next Section we are going to appreciate an importance of the derived gluing conditions and see how they appear in the QSC numerical algorithm.

Numerical solution
The equations of QSC are especially well-suited for numerical analysis: simple analytical properties of the P-functions allow to parametrize them in terms of a truncated Laurent series and then constrain these coefficients by the gluing condition. Numerical algorithms for solving QSC equations were developed and applied in [1,32,35,36]. In a non-symmetric case, such as BFKL with S 2 = n = 0, the procedure has to be modified in a way which we will describe here. We attached a Mathematica notebook named code_for_arxiv.nb implementing the algorithm, which we used to obtain the results described in this Section.
Let us start by briefly reminding the main steps of the numerical algorithm. A comprehensive description of the algorithm for the left-right symmetric case can be found in [1,18]. Here we will point out the main features we have to take into account in the case without left-right symmetry. As in the left-right symmetric case, for the P-functions there is a sheet with only one cut in the complex plane where the following parametrisation is valid . The expansions (3.1) contain only the even powers of Zhukovsky variable x(u) because for the state in question the P-functions possess the certain parity determined by their asymptotics from (2.13) and (2.51). However, the coefficients c a,k and c a,k are not independent and are subject to the conditions following from In the left-right symmetric case the condition (3.2) was satisfied automatically.
Since Q a|i is analytic in the UHP and has a power-like behaviour at u → ∞, its asymptotic expansion for sufficiently large Im u in the UHP can be written as Plugging the ansatz (3.1) and the expansion of Q a|i into the equation we are able to fix the coefficients B a|i,k in terms of the operator charges and the coefficients c a,k and c a,k . The fact that the P-functions have the certain parity and they are given by (3.1) leads to the disappearance of the odd coefficients B a|i,2l+1 = 0 for l = 0, 1, 2, . . . in (3.3) and we obtain After doing this, using the same finite difference equation in the form we find the numerical value of Q a|i in the vicinity of the real axis. One remaining ingredient of the iterative numerical procedure is the loss function -a function which is zero for the exact solution and which should decrease as each iteration brings us closer to the exact solution. We have the following loss function which is zero when the gluing condition is satisfied. Here and {u i } is a set of points on the interval [−2g; 2g]. Every function F i (u) depends on the charges S, ∆, n, the coefficients c a,k and c a,k and the coefficients of the gluing matrix. As a starting point for the numerical algorithm one can use the weak coupling data from Appendix B.
In the present work we are interested in the case of non-integer spin S 1 = S. As it was already shown in the Section 2 in this situation we cannot keep all the gluing conditions (2.80). However, we found that only two gluing conditions F 2 and F 4 are sufficient to constrain all the coefficients c and are still valid even for non-integer S 1 providing thus a natural way to analytically continue to non-integer spins. In terms of (3.8) this means that the sum in that formula goes only over i = 2, 4. After the loss function and all the constraints are formulated, the algorithm searches for the parameters which minimize the loss function subject to the constraints using a numerical optimisation procedure (Levenberg-Marquardt algorithm). Then, using the obtained numerical values of the Q-functions, we are able to restore the ansatz for the gluing matrix, which tells us which elements of the gluing matrix contain the exponential terms and which of them are equal to zero. This allows to verify the modification proposed in the Section 2 for the gluing matrix (3.10) (in agreement with [23]). Thus for integer S 2 = n this leads to the gluing matrix for non-integer S 1 given by (2.96) with In the situation when S 2 = n is not integer the gluing matrix (3.10) needs further modification. To achieve this we use the fact that relaxing the conditions (3.10) is sufficient to make the numerical procedure convergent. Restoring again the gluing matrix, for general real value of the spin S 2 (or conformal spin n in high-energy scattering terminology) we are left with the gluing matrix which coincides with (2.96).
Using the proposed numerical algorithm, we managed to calculate several numerical quantities for the cases when n is non-zero and even non-integer. On the Figure 3 one can find the length-2 operator trajectory for n = 1.
It is also possible to numerically calculate the dependence of the spin S on the coupling constant g for the fixed dimension ∆. On the Figure 4 you can see the dependence S(g) for ∆ = 0.45 and n = 1 in comparison with the same result calculated perturbatively as the sum of LO and NLO BFKL eigenvalues.
Additionally, this numerical scheme allows us to compare the numerical values of BFKL kernel eigenvalues with the known perturbative eigenvalues at LO and NLO orders. In the Table 1     Section 4 calculated in the small coupling regime. The continuous lines correspond to the strong coupling expansion of the intercept function from Section 6, which was fitted from numerical data obtained in the present Section. In the next Section we are going to analyze the weak coupling expansion of the intercept function. To achieve this we apply the iterative method first applied in [23].

Weak coupling expansion
In this Section we explore the function S(∆, n) perturbatively at weak coupling for arbitrary integer conformal spin n. In particular, we are interested in the BFKL intercept j(n) = S(0, n) + 1. The calculation of this quantity consists of two steps.
First, we apply the QSC iterative procedure introduced in [23] to the calculation of the intercept function for some integer n 10 . To do this we adopt this procedure to the case without the left-right symmetry. We repeat the main points of the iterative algorithm introduced in [23] and describe the functions which are used in it for the case n = 0.
In the second part we formulate an ansatz for the weak coupling expansion of the intercept function for arbitrary value of conformal spin in terms of binomial harmonic sums and fix the coefficients of this ansatz using the values of the intercept at several integer n, which we calculate solving the QSC iteratively order by order. This approach appears to be successful in the NNLO order in the coupling constant allowing us to find the intercept function at this order, but at NNNLO order we were not able to fix the rational part of the result for arbitrary n (see the details in the Subsection 4.2). This is due to the lack of the generalized-"reciprocity" at the NNNLO order. It would be very interesting to understand why and how the reciprocity in n is violated, which would allow to obtain the rational part with a greatly reduced basis of functions. We postpone this quesiton to the future investigation.
Let us also mention that the method we explain here should be also applicable for non-zero ∆ when ∆ + n takes odd integer values. The case of n = 0 was considered in [23], and it was sufficient to take a few values of ∆ in order to fix the NNLO dimension. This would be also interesting to investigate in the future.

Description of the iterative procedure
The iterative procedure increasing the number of orders is essentially the same as described in [23], but modified to account for non left-right symmetric case. The procedure is based on applying a version of variation of parameters method applied to the equation which is a simple consequence of (2.1) and (2.5). Indeed, suppose the function Q (0) a|i solves the equation (4.1) up to a discrepancy dS a|i The exact solution can be represented as the zero-order solution plus a correction. We expand the correction in the basis over the components of the zero-order solution If the discrepancy is of the order of some small parameter to the power m, i.e. dS a|i = O( m ), then we can obtain the equation for the coefficients b k i with the doubled precision The discrepancy thus becomes two orders smaller with each iteration step.
For the problem in question we restrict ourselves to the situation when ∆ = 0, n is an integer number and we perform the expansion in the parameter = g 2 . To start solving the finite difference equation (4.4), we have to find the zero-order in g solution Q where η s (u) are examples of the so-called η-functions, whose definition is given below. As we checked by solving the 4th order Baxter equation (2.34) for different values of the conformal spin n, for odd values of n, as for the case n = 0 described in [23], the Q-functions in the LO in g at ∆ = 0 can be expressed as linear combinations of the η-functions with the coefficients being Laurent polynomials 11 in u plus a Laurent polynomial in u without η-function multiplying it as in (4.5). The η-functions were introduced in [37] and then used in [23,28,35] with their generalized version in [32] and for ABJM theory in [38] with an application in [39] as where s i ≥ 1, i = 1, . . . , k. Unfortunately, for non-zero even n we were not able to determine the class of functions to which Q (0) i belong, therefore from this moment we restrict ourselves to the odd values of n. Then, applying the equation

7)
11 Laurent polynomial is a polynomial with both positive and negative powers of the variable plus a constant term.
where P (0) a are given by (B.11) with ∆ = 0, we can find the functions Q a|i in the leading order in the coupling constant.
To describe the solution to (4.7) we have to explain some properties of η-functions. This class of functions is particulary convenient because it is closed under all relevant for us operations. First, a product of two η-functions can be expressed as a linear combination of η-functions using the so-called "stuffle" relations [40]. Second, the solution of the equation of the form can be expressed as a linear combination of η-functions with the coefficients being Laurent polynomials in u. These two properties make η-functions very useful when solving the QSC perturbatively at weak coupling [23,28,37]. The described properties of the η-functions lead us to the conclusion that at least Q We call the operation inverting the linear operator in the LHS of (4.8) "periodization" (for a more precise definition see 5.2.2). It is easy to see that the "periodization" operation solves the equation (4.4) if the zero-order approximation entering the RHS is expressed as a linear combination of η-functions with the coefficients being Laurent polynomials in u plus a Laurent polynomial in u. We were able to find such representation for odd values of n, but not for even ones. After the zeroorder solution is found, we iterate it as described above, applying the operation of periodization to the RHS of (4.4) in order to find the coefficients b k i . Because of the two properties of η-functions mentioned above, at each iteration the solution is again obtained in the form of a linear combination of η-functions with Laurent polynomial coefficients plus a Laurent polynomial.
After finding the corrected Q a|i (4.3) at the given iteration step we still have some unfixed coefficients in it including the quantity of our interest S(0, n). To find them we calculate the Qfunctions from (2.6) and (2.5) and apply the gluing conditions for the case of integer conformal spin n, i.e. (2.96) with (3.10) satisfied. As it was explained in the Section 2 in these gluing conditions the Q-functions has to possess the pure asymptotics (2.68). Therefore we find the combinations of the Q-functions with the pure asymptotics. It appears to be sufficient to use only 2 of 4 gluing conditions, namelyQ 2 =M 12 1Q1 , (4.9) taking the Q-functions on the cut on the real axis. This procedure allows to fix the remaining unknown coefficients including the function S(0, n). The described method allowed us to find the values of the intercept functions for odd conformal spins in the range from n = 1 to n = 91 up to NNNLO order in the coupling constant. These data will be used in the next Subsection, where we put forward the ansatz for the structure of the intercept function for arbitrary value of conformal spin.

Multiloop expansion of the intercept function for arbitrary conformal spin
Using the procedure described in the previous Subsection 4.1 we have calculated the expansion of the BFKL eigenvalue intercept for odd n up to n = 91 in the weak coupling limit up to the order g 8 (NNNLO). These data are valuable by themselves, as they can serve as a test for future higher-order or non-perturbative calculations. What is more important, however, is that it allowed us to find NNLO and partially NNNLO BFKL eigenvalue intercept as a function of the conformal spin n.
We start by noticing that LO and NLO BFKL intercept can be represented as a linear combination of nested harmonic sums of uniform transcendentality. Indeed, the LO and NLO BFKL Pomeron eigenvalues themselves can be expressed (see, for example, [3] and Appendix C, where this calculation is explained in details) through the nested harmonic sums described, for example, in [41] S a1,a2,...,an (x) = is negative, the formula (4.10) holds only for even integer x. In the literature [42][43][44][45][46] there was described the analytic continuation of the harmonic sums in question from the positive integer even x.
To work with such a continuation we utilize the Mathematica package Supppackage applied in [23]. One can take ∆ = 0 in these eigenvalues, which after some simple algebra gives Here and below the transcendentality is computed as follows: the transcendentality of a product is assumed to be equal to the sum of transcendentalities of the factors and transcendentality of a rational number is 0. Transcendentality of log 2 is 1 and transcendentality of ζ k is k. Since ζ k for even k is proportional to π 2 , it is easy to see that transcendentality of π is 1.
To conduct the calculations with harmonic sums one can use the HarmonicSums package for Mathematica [46][47][48][49][50][51][52] or the Supppackage utilized by the authors of [23]. It should be noted, that in the present work we utilize the same conventions for the analytic continuation of harmonic sums as in the latter work [23]. As we see, the argument of all the harmonic sums in (4.11) is (n−1)/2. This leads one to an idea of trying to find NNLO and NNNLO intercepts as analogous linear combinations of harmonic sums with transcendental coefficients of uniform transcendentality. The coefficients of the linear combination can be constrained using the data generated by the iterative procedure. But the number of harmonic sums of certain transcendentality grows fast as transcendentality increases. Fortunately, one can drastically reduce the number of harmonic sums in the ansatz by conjecturing a certain property of the result we call reciprocity.
The property in question [53][54][55][56] is parallel to the Gribov-Lipatov reciprocity [57,58] and was observed in the weak coupling expansion of the scaling dimensions of the twist operators. Let us remind the statement of the reciprocity: if one defines an auxiliary function P [53][54][55]  Thus the function P(M ) is much more convenient to work with than γ(M ) itself: the function γ(M ) can be expressed through the nested harmonic sums, while P, on the other hand, can be expressed through a much smaller class of functions which satisfy the property (4.13). Such functions were identified and used in [54,[59][60][61] as reciprocity-respecting harmonic sums. However, in our calculations we use another basis of the functions satisfying (4.13), which was applied in the works [29,[62][63][64][65][66][67]. These functions are called the binomial harmonic sums and for integer M they are defined as (see [52]) (4.14) Note that we consider only the positive indices i l , l = 1, . . . , k in the definition (4.14). Those are exactly the sums whose asymptotic expansion is even at infinity after the argument is shifted by 1/2. All this is directly applicable to our case and we are able to formulate an ansatz for the NNLO intercept function. From the LO and NLO expressions (4.11) we see that their asymptotic expansions at large n are even in n. Since we are using the harmonic sums of the argument M = (n − 1)/2, we need to keep only the harmonic sums invariant under the transformation M → −1 − M or n → −n in our notations. Those are exactly the binomial sums (4.14). The expressions (4.11) for LO and NLO intercepts can be easily expressed through them where the arguments of the sums are again (n − 1)/2. In order to find the NNLO intercept we make an ansatz in a form of a linear combination of binomial harmonic sums with transcendental coefficients. The maximal transcendentality principle, formulated by L.N. Lipatov and A.V. Kotikov [2,3], holds for the intercept as well: every term in the sum should be of the total transcenedentality 5. The terms of the sum can of course be multiplied by arbitrary rational coefficients which do not affect the transcendentality. Having constructed the ansatz in this way, we can now constrain its coefficients by the iterative data: we evaluate the ansatz (a linear combination of binomial nested harmonic sums) at several integer values of n and match the result to the data obtained from the numerical procedure for the corresponding n. Equating the coefficients in front of each unique product of transcendental constants in these two expressions, we get a linear system for the rational coefficients of the ansatz. Solving it we obtain a surprisingly simple expression The result (4.16) for the intercept function for arbitrary n can be compared with the other known quantities. First of them is the NNLO BFKL Pomeron eigenvalue for the conformal spin n = 0 calculated in [23]. Taking in this eigenvalue ∆ = 0 and comparing it with (4.16) for n = 0 we see perfect agreement. Second, for non-zero conformal spins the formulas for the Pomeron trajectories were found in [27], from which we can extract the intercept for given n. We also checked that the result of that work coincides with our result (4.16) for several first non-negative conformal spins n, thus representing an independent confirmation of the correctness of our calculation. The same procedure can be repeated in the NNNLO. The values of the NNNLO intercept for several first odd values of the conformal spin n are given in the Appendix D. Again, as for NNLO, an ansatz in a form of a linear combination of binomial harmonic sums with transcendental coefficients of uniform transcendentality 7 can be constructed and we attempted to fit it to the iterative data. However, we found that the basis of binomial harmonic sums is insufficient to fit the data. This signals that reciprocity understood as parity under n to −n seems to be broken down in this case. The reasons for this are unclear and will be the subject of further work. However, we managed to fit the certain part of the NNNLO data.
For each odd n we calculated (see Appendix D for several first conformal spins n) the value of the NNNLO intercept is a linear combination of the transcendental constants consisting of π, ζ 3 , ζ 5 and ζ 7 with rational coefficients and a rational number (see the file intercept_values_Nodd.mx with the data for the odd conformal spins from n = 1 to n = 91 in the arXiv submission of this paper). Let us restrict ourselves to the values of the conformal spin in our data equal to n = 4k + 1 with k = 0, 1, . . . , 22. In these points the values of the intercept functions are given by where all coefficients in front of the transcendental constants on the RHS of (4.17) are rational functions of k . Each coefficient in the RHS of (4.17) is conjectured to be be a linear combination with rational coefficients of the binomial harmonic sums with the transcendentality, supplementing the transcendentality of the corresponding coefficient to 7. We were able to fit all the contributions except for j ζ3 N N N LO and j rat. N N N LO , which, as other harmonic sums, take rational values at the points n = 4k + 1 for integer k ≥ 0. However, we found that the term j ζ3 N N N LO cannot be fitted with the ansatz consisting of the binomial harmonic sums (4.14). This motivated us to try to fit this contribution with the nested harmonic sums (4.10). This appeared to be really the case and we managed to fit this part with the ordinary harmonic sums, which means that the reciprocity, i.e. the symmetry n → −n in the asymptotic expansion, is violated. For the last, rational contribution j rat. N N N LO , we also found that it is not described by the binomial harmonic sums. Unfortunately, fitting this contribution with the ordinary harmonic sums did not lead us to completely fixing this contribution due to the lack of data. Therefore combining the obtained results we write down the non-rational part of the answer for the points n = 4k + 1, which is the sum of the terms in the RHS of (4.17) except for j rat.
One can find the full values of the NNNLO intercept function including the rational terms for the conformal spins n = 4k + 1 from 1 to 89 in the arXiv submission of this paper in the file named intercept_values_Nodd.mx. As the part of the harmonic sum consisting of the nested harmonic sums of the transcendentality 7, which constitutes the rational part at the points n = 4k + 1, was not fitted we are unable to write an expression for the NNNLO intercept working for all conformal spins leaving this task for future studies. Let us briefly summarize the results of the present Section. In the first Subsection applying the QSC iterative algorithm we found the values of the intercept up to NNNLO order in the coupling constant in the certain range of odd values of the conformal spin. In the second Subsection we saw that the LO and NLO intercept functions satisfy the reciprocity symmetry, which allowed us to rewrite them in terms of the binomial harmonic sums and using the ansatz in terms of these sums in the NNLO order, fix the answer for the NNLO intercept for arbitrary conformal spin. In the NNNLO order the reciprocity breaks down, nevertheless, for the conformal spins n = 4k + 1 we managed to describe the non-rational part of the NNNLO intercept function in terms of binomial and ordinary nested harmonic sums.

Near-BPS all loop expansion
In this Section we are going to analyze the QSC equations near the BPS point ∆ = 0, S 1 = −1 and S 2 = 1. It appears that it is possible to calculate two non-perturbative quatities in this point by the methods of QSC. Another BPS-point S 2 = 0, S 1 = 0 was analyzed in detail in [21,68,69]. In this Section we follow closely the Near-BPS expansion method by [21].

Slope of the intercept near the BPS point
A particularly important role in BFKL computations is played by the intercept function j(n) = S(0, n) + 1, where S = S 1 and n = S 2 . As we mentioned above, the point ∆ = 0, n = 1 is BPS, by which we mean that it is fixed for any 't Hooft coupling. The group-theoretical argument explaining this phenomenon should be based on the shortening condition. From the QSC perspective the BPS points are the points where A a A a = 0 simultaneously for all a = 1, . . . , 4. In this Subsection we study small deviations from this BPS point and calculate the slope of j(n) with respect to n in the point n = 1 in all orders in the coupling constant g.

LO solution
The fact that at the BPS point A a A a = 0 usually leads to more powerful condition P a = P a = 0, which is known to lead to considerable simplifications [21] (see also [70] for a similar simplification in the TBA equations). Based on that we also expect that in our situation the Q-system simplifies a lot near the BPS point. Let us show that, indeed, Q a|i takes a very simple form for S 2 = 1, ∆ = 0.
Recall that the functions Q a|i satisfy the equation Let us look at how the right and left hand sides of the equation (5.1) behave as S 2 approaches 1. Scaling of P a and Q i can be deduced from their leading coefficients A a and B i . For the convenience of the calculations we perform the rescaling of the P-functions. The H-symmetry (2.8) and its particular case rescaling symmetry (2.10) explained in [17] allow us to set some of the coefficients from (3.1) to some fixed values. To describe them we introduce a scaling parameter ν = √ S 2 − 1, which is real for S 2 ≥ 1. For the coefficients of the P-functions we obtain For the Q-functions in their turn As follows from the AA-, BB-relations (2.17), (5.2) and (5.3) in the small ν limit Notice that then from (2.17), (5.2) and (5.3) we derive that A a and B i all scale as ν and are given by the formulas is the slope-to-intercept function, which is the quantity of interest in the present Subsection. This means that the right-hand side in (5.50) is small and the functions Q a|i are i-periodic in the LO.
Since they are also analytic in the upper half plane, they are analytic everywhere. Recall that the asymptotics of Q a|i are given by  When ∆ = 0 and S 2 = n is arbitrary, the left-right symmetry is restored, which simplifies the solution a lot. Despite the normalizations (5.2) and (5.3) differ from (2.52) and (2.53) by a rescaling symmetry one can see that the left-right symmetry is restored (which is also confirmed by the weak coupling data for arbitrary ∆ and S 2 from Appendix B). At this point the symmetry (2.54) takes the form where χ ab and η ij are given by (2.55) and χ ab is the same matrix as for the left-right symmetric states as in [16,17]. As we consider the case when both spins S 1 and S 2 are not integer, let us recall the gluing matrix for this case (2.96) taking into account the hermiticity of the gluing matrix and keeping in mind that M 11 1 , M 33 1 and M 44 1 are real. As now we have the matrix η from (2.55), which relates the Q-functions with lower and upper indices, it is possible to obtain an additional constraint on the gluing matrix. Plugging this relation into the first gluing condition (2.29) we derivẽ which after comparison with the second gluing condition from (2.29) leads us to (M −t ) ij = η ik M kl η lj , and after multiplication of the both sides by (M t ) mj we obtain an additional constraint for the gluing matrix Substitution of (2.82) into (5.12) leads us to the following equations 14) It should be noted that the first equation (5.13) is satisfied for the gluing matrix (2.96) (the same as in (5.10)).
To start solving the constraint (5.12) order by order in ν we use the following expansion of the gluing matrix, which is motivated by the scaling of the Q-functions (5.4) and (5.5) are pure imaginary. Let us start solving these equations in the LO. We have already found Q a|i and thus As the scaling of the Pand Q-functions is determined by the scaling (5.4) and (5.5) of the leading coefficient at large u, we are left with the following expansions of these quantities in the small ν limit (5.24) Then, the correspondence between Qand P-functions (2.5) and (2.6) in the LO taking into account (5.8) looks as follows where In what follows we will usually omit the superscript with the variable u if the context does not imply the usage of cosh ± and sinh ± with different arguments and the expression 4πg in the argument of the generalized Bessel function for the sake of conciseness. Substituting P 3 from (5.27) into the first equation of (5.26), we havẽ We need to find the solution for P As the LHS of (5.35) does not contain the even powers of x for the RHS we have The equation (5.35) takes the form which still depends on the unknown coefficient M (0) 14 2 , which we will fix in the next Section, but going to the next order.
So far, starting from the gluing conditions for non-integer conformal spin n (5.10) we managed to solve the constraints on it in the LO getting (5.22) and then, using the connection between the Qand P-functions in the LO formulated the system of equations (5.26) for the P-functions in the LO. Then from this system we found all the P-functions in the LO except for P 2 for which we derived the equation (5.37), still containing one unknown constant. In the next Subsection using this equation we are going to show how to find the slope-to-intercept function.

Result for the slope-to-intercept function
From now on let us consider the equations (2.1) for the Q a|i functions in the NLO. We start with Q 3|3 should not contain log u in its large u asymptotic, which can only appear from the expansion of non-integer powers Au λ A + Aλ log u + . . . . Thus, in the large u expansion of the RHS of (5.38) the coefficient in front of 1/u (which would produce log u in Q (1) 3|3 ) has to be equal to 0, which is guaranteed by The substitution of (5.39) into (5.37) gives (5.40) We do not need to solve completely (5.40), because only the first coefficient of the expansion of P 2 in the powers of x is sufficient to find the slope-to-intercept function Remembering the leading coefficient of the large u asymptotic of P 2 from (5.4) we obtain the equation The previous equation (5.42) fixes θ, which is now equal to θ(g) = 1 + and constitutes our result for the slope-to-intercept function.
The weak coupling expansion of the obtained result (5.43) is given by Since the slope-to-intercept function by definition is the derivative of the intercept function with respect to the conformal spin n at n = 1 we can immediately compare the first few coefficients in the weak coupling expansion (5.43) with the derivative of (4.11) and (4.16). This will also show that these expressions provide the formulas compatible with our analytic continuation in n away from the integer values. To find the derivatives of the binomial harmonic sums with respect to the argument we apply the SuppPackage used in [23], which expresses the nested harmonic sums (4.10) in terms of the η-functions (4.6) and allows to find the derivatives of these sums. The derivatives of the intercept function (4.11) and (4.16) can be calculated and we find that that they are in full agreement with (5.44) confirming our result (5.43).
In the next Section we compute the strong coupling expansion of our result for the slope-tointercept function. As we will see the calculation is less straightforward than at weak coupling, even thought the result is still quite simple.

Strong coupling expansion of the slope-to-intercept function
To obtain the strong coupling expansion of the slope-to-intercept function (5.43) first we calculate the following expansion at strong coupling If we just sum the series (5.46) multiplying it by (−1) k , it appears to be divergent. However, simple ζ-regularization gives the right result, as we verified numerically with high precision. Namely, multiplying this expression by k δ and understanding the result as the limit δ → 0, we get the following answer Then, substituting (5.48) into the expression for the slope-to-intercept function (5.43), we obtain the strong coupling expansion for it This expansion (5.49) will be useful for us in the Section 6 when we are able to compare it with the derivative of our formula for the strong coupling expansion of the intercept function for arbitrary conformal spin n taken at the point n = 1, which is based on intensive numerical analysis.

Curvature function near the BPS point
As it was mentioned above in the Subsection 5.1 we are considering the expansion in the vicinity of the BPS point ∆ = 0, S 1 = −1 and S 2 = 1. In the previous Subsection we expanded in the powers of S 2 − 1, however, to find the curvature we keep S 2 = 1 and expand in the powers of ∆. The scheme of solving the QSC equations in this case is similar but with one difference. Since the function S(∆, 1) is an even function of ∆ we have to expand to the NLO order in ∆ to get a non-trivial result.
As in the previous Section we will utilize a simplification in the Q a|i function to solve the Q-system explicitly.

LO solution
Recall that the functions Q a|i satisfy the equation Let us look at how the right and left hand sides of the equation behave as ∆ approaches 0. Behaviour of P a and Q i can be deduced from their leading coefficients A a and B i . Analogously to the case of the slope-to-intercept function the H-symmetry (2.8) and its particular case rescaling symmetry (2.10) explained in [17] allow us to set some of the coefficients from (3.1) to some fixed values.
To describe them we introduce a scaling parameter = √ ∆, which is real for ∆ ≥ 0. For the coefficients of the P-functions we obtain For the Q-functions in their turn As follows from the AA-, BB-relations (2.17), (5.51) and (5.52) in the small limit Note that due to the parity symmetry ∆ → −∆ we expect α = 0, which will be the consistency check of our calculation. We only get a non-trivial result in the NLO order. Then from (2.17), (5.51) and (5.52) we derive that A a and B i all scale as and are given by the formulas This means that the right-hand side in (5.1) is small and the functions Q a|i are i-periodic in the LO in . Since they are also analytic in the upper half plane, they are analytic everywhere. Recall that the asymptotics of Q a|i are given by the result is the same as in the case of the slope-to-intercept function. As in the case of the slope-to-intercept function we can exploit the symmetry between the Pand Q-functions with lower and upper indices. Since ∆ is different from 0 we remember the symmetry (2.54). Despite the normalizations (5.51) and (5.52) differ from (2.52) and (2.53) by a rescaling one can see that the symmetry (2.54) takes the form (which is also confirmed by the weak coupling data for arbitrary ∆ and S 2 from Appendix B) where we set (−∆) = √ −∆ = i √ ∆ = i choosing the branch of the square root with the cut going from 0 to +∞ and (5.59) Usage of (5.58) and the substitution of (5.53) and (5.55) leads us to the equality as it should be. Substituting α from (5.60) into (5.53) and (5.55) we obtain A's and B's in the LO As we consider the case when the spin S 1 is not integer and the spin S 2 is integer, we can use the gluing conditions (3.10) from the Section 3, which were also mentioned in the Section 2, thus keeping in mind that M 11 1 and M 33 1 are real. As now we have the matrix η c from (2.55), which relates the Q-functions with lower and upper indices with ∆ replaced by −∆, it is possible to obtain an additional constraint on the gluing matrix. Plugging the relation (5.58) into the first gluing condition (2.29) we derivẽ Substitution of (2.82) into (5.64) leads us to the following equations It should be noted that the first equation (5.65) is satisfied for the gluing matrix (5.62) (the same as in (5.10)).
To start solving the constraint (5.64) order by order in analogously to the case of the slopeto-intercept function we use the following expansion of the gluing matrix, which is motivated by the scaling of the Q-functions (5.53) and (5.55) are pure imaginary. Let us start solving these equations in the LO. We have already found Q a|i and thus As the scaling of the Pand Q-functions is determined by the scaling (5.53) and (5.55) of the leading coefficient at large u, we are left with the following expansions of these quantities in the small limit Then, the correspondence between Qand P-functions (2.5) and (2.6) in the LO taking into account (5.57) looks as follows To summarize this part, we started from the gluing conditions for integer conformal spin n (5.62) and we managed to solve the constraints on it in the LO getting (5.74) and then, using the connection between the Qand P-functions in the LO formulated the system of equations (5.78) for the P-functions in the LO. Then from this system we found all the P-functions. In what follows using these functions we are going to show how to find the NLO solution.

NLO solution
Since the function S(∆, n) is even, we have to consider the next-to-leading order. Our strategy is to find the P-functions in the NLO order and extract the quantity of interest -the curvature function -from the leading coefficients of these functions. Indeed, from (5.51) and (2.17) we derive is the curvature function. We begin by finding the correction to Q (0) a|i . The equation (5.50), taking into account (2.5) and the scaling (5.76) and expanding it in the NLO, takes the form As the functions P (0) a were completely fixed in the previous calculations we are able to determine all Q (1) a|i up to a constant. In solving (5.94) we are going to act in a way similar to the one presented in [21]. We have to find an UHPA solution to the equation of the form where the RHS has one cut on the real axis and it can be represented as a series in the Zhukovsky variable x(u), whose powers are bounded from above. To build the solution of such an equation let us rewrite the RHS in the following form where h pol (u) is a polynomial in u and h − (u) is a series in the Zhukovsky variable x(u) starting from the power not greater than x −1 (u). Since h(u) is a series in x(u) bounded from above, we can always rewrite x a (u) with a > 0 as where x a + 1/x a is a polynomial in u. Thus, using (5.97) we are able to replace any positive power of x(u) in h(u) with the difference of polynomial in u and the negative power of x(u) which justifies the form (5.96). The UHPA solution of (5.95) is given by and c is a constant. The operator Γ U is determined by where the integration contour around the cut goes clockwise. Also we determine the operator Γ D which we will use to obtain the LHPA solution It is not hard to check that Rewriting the RHS of (5.94) in the form (5.96) we can find the solution of it  In what follows we will also need the complex conjugated functionQ (1) a|i , for which we have the equationQ where In finding the solution of (5.113) we follow the same method as in [21]. Let us briefly sketch its main points. One can notice that the equations for the P-functions in the NLO (5.113) have one of the two formsF (u) + F (u) = G(u) andF (u) − F (u) = G(u) , (5.115) where F (u) is a power series in Zhukovsky variable x(u) and the function G(u) can be represented as Laurent series in Zhukovsky variable x(u) in the vicinity of the point x(u) = 0. As it can be seen, the equations (5.115) are self-consistent only if, respectively, the conditions G(u) − G(u) = 0 andG(u) + G(u) = 0 (5.116) are satisfied. As it was explained in [17], the unique solutions to the equations (5.115) on the classes of functions non-growing and decaying at infinity are respectively where the integral kernels are given by It should be noted that if the asymptotic of F (u) is not specified to be non-growing of decaying respectively the solutions of (5.115) may include the zero-modes, i.e. the solutions of the homogeneous equations (5.115) with zero RHS.
With the usage of the integral kernels H(u, v) defined in (5.118) first the 3rd and 4th equations of (5.113) are solved, then the obtained P are satisfied exactly and do not fix any constants in the RHS of (5.113). But from the 1st and 2nd equations of (5.113) we obtaiñ . One can find the details of this calculation and the values of these constants in Appendix F. Therefore we checked that the RHS of (5.113) satisfy (5.116) and we can apply (5.117).
As P 3 has the constant asymptotic at infinity, the unique solution for P 4,1 and the curvature function γ(g) defined in (5.93), which we have to fix. With the usage of (5.51) and (5.92) we can express M 4,1 = where the integral kernel Γ is defined as and H k is the k-th coefficient in the large u expansion of the convolution of the integral kernel H(u, v) with the function on which it acts. The asymptotic of P  the powers x(u) and 1/x(u) and thus it is proportional to x(u) − 1/x(u). The coefficient in front of this zero-mode is determined by the coefficients (5.92), then P (1) 4 is equal to and contains the constants M (1) 34 1 and γ, which are not fixed yet. To proceed with the solution of the 3rd and 4th equations of (5.113) we take each equation and subtract the same equation conjugated As both the functions P 1 (u) and P 2 (u) have the decaying at infinity asymptotics, the solutions of (5.125) are uniquely determined by the formula (5.117) with the kernel K(u, v). Calculating R 3 with the usage of (5.121) and (5.122), we find the solution for P Doing the same for R 4 with the usage of (5.124), we find the solution for P where K k is the k-th coefficient in the large u expansion of the convolution of the kernel K(u, v) with the function on which it acts.
To sum up, we managed to solve the constraints for the gluing matrix in the NLO, which allowed us to write down the system of equations for the P-functions in the NLO (5.112). In the NLO some Q-functions contain infinite series of cuts, which led us to the usage of the integral kernels (5.102), (5.103) and (5.118). After finding the P-functions in the NLO by solving the system for them we are ready to fix the curvature function by utilizing the values of the coefficients found in the previous calculations.

Result for the curvature function
To obtain the curvature function, we have to remember that c 4 . Comparing this coefficient found from (5.124) with the usage of (5.129) and the one from (5.122) we find the answer (5.130) By expanding the integral kernels (5.118) at large u we can rewrite the curvature function in a more concise form After obtaining (5.131) we can compare it with the other known results. In the first two orders we know the BFKL Pomeron eigenvalues for arbitrary conformal spin including n = 1, therefore we are able to calculate the curvature from these eigenvalues in these two first orders. Comparing with the weak coupling expansion of the formula (5.131) for the curvature function −56π 2 ζ 7 − 6930ζ 9 g 8 + 136π 8 2835 ζ 3 + 668π 6 189 ζ 5 − 112π 4 3 ζ 7 + 508π 2 ζ 9 + 93720ζ 11 g 10 + + 754π 10 42525 we find that the first two terms of the expansion (5.133) coincide with the values obtained from the BFKL Pomeron eigenvalues, which represents a check of our result. Having computed the weak coupling expansion (5.133), in the next part we analyze the other interesting limit, i.e. the strong coupling expansion, which is a separate computational task in the case of the curvature function.

Strong coupling expansion of the curvature function
For the calculation of the strong coupling expansion of the curvature function (5.131) we utilize the same method as the one used in [21]. Let us briefly sketch the scheme of this calculation. The curvature function can be represented in the following form Integrating by parts (5.134) and changing the integration variables to x u and x v respectively, we obtain where the integrals go over the unit circle and G(x u , x v ) is a polynomial in the variables x u , x v , 1/x u , 1/x v , cosh u − , cosh v − and sinh u − . Then, expanding cosh − and sinh − in a series in the Zhukovsky variable x, we can represent (5.135) as an infinite series in the coefficients of the BES dressing phase [71][72][73][74], which admits the large g expansion in the form of asymptotic series [71].
Using the computed series we calculated the numerical value of the curvature function for a number of points in the range 5 ≤ g ≤ 40 and obtained the first coefficients of the strong coupling expansion of the curvature function with high precision by fitting the data with inverse powers of g. The first several coefficients appear to be simple rational numbers. In the 6th coefficient it was expected to have ζ-function, which we managed to fit utilizing the EZFace [75] webpage. As in the case of the slope-to-slope function [21], usage of the exact value fit in the order 1/g k increases accuracy of the fit in the next order 1/g k+1 , which represents a non-trivial check of the validity of the used method.
Application of the numerical method described above yields us the large g expansion of the curvature function where λ is given by (5.47). Together with the curvature function (5.131) and its weak coupling expansion (5.133) the formula (5.136), containing the strong coupling expansion of this function, concludes the list of the results of the present Section.

Intercept function at strong coupling
The other interesting limit of the intercept function besides the small coupling limit considered in the Section 4 is the strong coupling one. For n = 0 case the intercept function in the strong coupling limit was analyzed in [41,[76][77][78] and then extended to the next orders by the QSC method in [21]. As we have already the numerical data for the intercept for the different values of conformal spin n, then using the numerical algorithm described in the Section 3 and assuming that the coefficients are some simple rational numbers and extrapolating the high precision numerical data by the inverse powers of λ 1/2 we extract where λ is given by (5.47) and the result for n = 0 is taken from [21]. We see that the leading term is linear in n, sub-leading is quadratic and so on. This pattern is quite typical at strong coupling (see for example [31]). Assuming this polynomial pattern from the above data (6.1) we get S(0, n) = −n + (n − 1)(n + 2) which we can cross-check with our slope-to-intercept function (5.49) by differentiating it at n = 1! Comparison with (5.49) shows us complete agreement. We also verified our result numerically by taking n = 1.5 and fitting the data with inverse powers of λ 1/2 we reproduced precisely the coefficnets from (6.2).

Conclusions and future directions
BFKL regime is traditionally one of the "hard" problems of high energy theoretical physics. QSC method allowed to make progress in this direction when the traditional perturbative methods become too complex to make one's way through. In our paper we extended the area of applications of QSC even further, adding an additional parameter -conformal spin n. Importantly, we can deal with the situation when conformal spin takes an arbitrary real value. QSC can now handle the situation when all three global charges corresponding to AdS 5 -S 1 , S 2 and ∆ are non-integer. This required us to allow four components (or two independent if we take into account the hermiticity) of the gluing matrix to have exponential asymptotics. Gluing conditions are thus the main ingredient of the analytical continuation of QSC presented in this paper.
The knowledge of the gluing matrix allows us to implement an efficient numerical algorithm which solves the Pomeron eigenvalue problem numerically for any values of ∆, n and at any 't Hooft coupling g, limited only by the computing time 12 . Our Mathematica implementation of this method is attached to this arXiv submission.
As an illustration of our method we have computed two all-loop quantities -slope of BFKL intercept at ∆ = 0 with respect to n around n = 1 and curvature of the length-2 operator trajectory in the vicinity of the point ∆ = 0 and n = 1. We generated analytical perturbative and numerical data using the iterative procedures described in [23], which we had to modify to take into account exponential asymptotics for non-integer global charges. The iterative perturbative data allowed us to fix the NNLO intercept in terms of binomial harmonic sums and partially fix the NNNLO BFKL intercept in terms of binomial and ordinary nested harmonic sums.
Several further directions of work come to mind immediately. First, the basis of nested harmonic sums seems to describe well the perturbative expansion of BFKL eigenvalue and in particular its intercept. Since the iterative procedure can be used to arbitrary high order in g and for arbitrary odd n, it is just a question of time and computational power to fix BFKL eigenvalues at higher orders. The first task would be to find a closed analytic expression for the NNNLO BFKL Pomeron eigenvalue for arbitrary n and ∆. Since as we saw the reciprocity in this order is broken, it is an interesting open problem to understand the nature of this.
The second direction one can pursue is exploring the neighbourhood of the BPS point n = 1, ∆ = 0. After computing the slope of intercept with respect to n, one can compute the next order (n − 1) 2 . The computation should be similar in spirit of the computation of S 2 correction to the twist operator anomalous dimension [21].
Furthermore, one should study in a similar way the Odderon spectrum, which is expected to correspond to the length parameter to be taken L = 3. Most of the steps described in this paper should be applicable for arbitrary L and it would be very interesting to reproduce previously known perturbative results and extend them to finite coupling.
Finally, there are indications [79] that the structure constants can be also governed by the QSC Q-functions, which were evaluated in this paper in various regimes. So it would be interesting to compare and extend to finite coupling the results on the triple Pomeron vertex [80].

A Details of QSC construction
To show that the matrix Ω j i =Q − a|i C a b Q b|j− is i-periodic let us apply the conjugated equation (2.15) and the same equation for Q a|i which is valid now in the whole complex plane as we know the functionsQ i and Q i on the sheet with the short cuts. We have To show that the matrix Θ j i (u) = (−1) a+1 Q − a|i (−u)Q b|j− (u) is i-periodic let us apply the equation (2.15) with u replaced by −u and the same equation for Q a|i which is valid now in the whole complex plane as we know the functions Q i (−u) and Q i (u) on the sheet with the short cuts. We have B Weak coupling solution of the QSC for length-2 operators with nonzero S 2 = n In this Appendix we consider BFKL limit of the QSC with non-zero conformal spin. Let us first of all briefly remind what BFKL limit is. We are going to study the regime when at the same time the coupling constant g → 0 and one of the spins S 1 = S → −1 while keeping the ratio g 2 /(S + 1) finite. LO BFKL in this limit corresponds to resumming all the contributions of the form (g 2 /(S + 1)) k , NLO BFKL -to the contributions of the form (S + 1)(g 2 /(S + 1)) k and so on. However, in the present Subsection we take the second spin (conformal spin) S 2 = n = 0 and there appear some differences from the BFKL regime with zero conformal spin.
In order to find the BFKL kernel eigenvalue we are going to utilize the old fashioned method of Pµ-system. The Pµ-system consists of the functions P a (u), P a (u), which we introduced before and of an antisymmetric matrix µ ab (u) (see [16,17] for the detailed description). They satisfy the following equationsμ ab − µ ab = P aPb − P bPa ,P a = µ ab P b , (B.1) µ ab − µ ab = P aPb − P bPa ,P a = µ ab P b , P a P a = 0 , µ ab µ bc = δ c a ,μ ab = µ ++ ab ,μ ab = µ ab++ .
Before proceeding we are going to introduce a couple of new notations. Our notation for the BFKL scaling parameter is w = S + 1. It is also convenient to introduce the notation Λ = g 2 /w. To start solving the Pµ-system in the BFKL regime we have to determine the scaling of the Pand µ-functions in the limit w → 0. In what follows we are going to use the arguments from [22], where the left-right symmetric case with zero conformal spin was considered. First, we assume that the scaling of the P-functions coincides with the scaling of their leading coefficients in the large u asymptotics. Thus as from (2.17) for the length-2 state in question in the BFKL limit A a A a = O(w 0 ) for a = 1, . . . , 4 these functions can be chosen to scale as w 0 Second, because the asymptotic of the function P 1 (u) is the same as in the left-right symmetric case, for S 2 = n = 0 the argument from [22] about the scaling of the µ-functions is applicable and they scale as w −2 Additionally, as all the P-functions for the length-2 states being considered possess the certain parity from the Pµ-system equations (B.1) we can conclude that the functions µ + ab (u) have the certain parity, which will be specified below.
Finally, we are to determine the asymptotics of the µ-functions for the case of non-integer S 1 = S for non-zero S 2 = n. Let us restrict ourselves from now on in this Section to the case of integer conformal spin S 2 = n. The µ-functions with the lower indices are given by the formula from [17] combined with (2.44) from the Section 2 Substituting the matrices (3.10) and (2.72), which are applicable to the case of integer conformal spin, into (B.4), we are able to determine the asymptotics of the ω-functions. Taking into account also the antisymmetry of ω ij = M ik Ω k j from (2.50) we find that they are Using (B.5) we see that the leading asymptotics of the µ-functions with the lower indices at u → ±∞ are while the µ-functions with the upper indices have the same asymptotics but in the reverse order.
In addition, we know that the P-functions have only one cut on one of the sheets, therefore they can be written as a series in the Zhukovsky variable. The parametrization of the P-functions for the case J 1 = 2 and J 2 = J 3 = 0 is already known to us and can be taken from the formula (3.1). Utilizing the P-functions rescaling (2.10) we can set the coefficients A 1 = 1, A 2 = 1, A 3 = −1 and A 4 = 1. Then, we are also allowed to apply the certain H-transformation of the P-functions, which do not alter their asymptotics and parity As (B.8) has two parameters then, applying this transformation, we can set the coefficients c 3,1 and c 2,1 to zero. Thus, we arrive to the formulas In what follows for the sake of convenience we change the numeration of the µ-functions to µ 12 = µ 1 , µ 13 = µ 2 , µ 14 = µ 3 , µ 23 = µ 4 , µ 24 = µ 5 and µ 34 = µ 6 and the same for µ with the upper indices. Having set the scaling and the expansions of Pand µ-functions in the scaling parameter together with their asymptotics we are able to proceed with the solution of the Pµ-system order by order in the scaling parameter.

B.1 LO solution
In the present part we will find the LO solution of the Pµ-system in the BFKL regime. Taking into account thatP scale as 1/w 2 , in the LO in w we obtain for the P-functions from (B.9) and (B.10) Λu , (B.11) where we implied c 4,1 = Let us now examine the µ-functions in the vicinity of the point S = −1. In the LO, when S = −1, µ-functions with the lower indices at u → ∞ have the asymptotics (µ 1 , µ 2 , µ 3 , µ 4 , µ 5 , µ 6 ) ∼ (u 0 , u 1 , u 2 , u 2 , u 3 , u 4 )e 2π|u| . (B.13) To find the µ-functions with the lower indices in the LO analogously to the left-right symmetric case we notice that P where P(u) is an i-periodic function. Let us first determine only the polynomial part of the solution.
Plugging the ansatz (B.14) into the equations µ ++ ab = µ ab + P a µ bcP c − P b µ acP c we fix all the coefficients besides b 1,1 and c 1, 4,1 . The latter equality automatically leads to the validity of the requirement in the LO. Now let us return to the general solution. Because in the leading order the equations for µ's are homogeneous, the solution is determined up to some i-periodic function. This periodic function grows not faster than e 2π|u| , then the most general ansatz for it is Also we have to remember about the requirements of analyticity, i.e. that the following expressions have no cuts The constraints (B.17) lead us to the results Having obtained the LO solution given by (B.11), (B.19) and (B.12) we are able to proceed with finding the NLO solution of the Pµ-system together with the coefficient B 1 , which will be done in the next Subsection.

B.2 NLO solution
Now it remains to fix the coefficient B 1 and find the P-functions in the next order in w. To do this we need to calculateP. First, we have to understand, which coefficients it is necessary to find in order to calculate P in the NLO order w 1 using the scaling or the P-(B.2) andP-functions (w −2 ) Notice, that we already know c Thus, we have to find c 4,1 , c 4,2 , c 1,1(2) , c 1,2(0) and c 3,2(0) . Comparing theP-functions calculated from the LO result and the ones found with the equationP a = µ ab P b we can fix some unknown coefficients. On the one hand we get On the other hand The result is the following Thus, the final answer for the LO µ-functions is Using the equation we are able to determine the µ-functions with upper indices in the LO up to a multiplication constant. Then, inserting them into Pµ-system equations, we fix this constant obtaining and some other coefficients in the P-functions with upper indices c 1,2(0) = i 96 ∆ 2 − n 2 2 − 2 ∆ 2 + n 2 + 1 , c 3,2(0) = −1 . (B.28) We are able even to go further and calculate completely the µ-functions in the NLO and partly fix the expansion coefficients in the P-functions in the NNLO order. Let us now fix the other coefficients in the NLO order in w. To find the unknown coefficient c (2) 4,1 , we have to build an ansatz for µ in the NLO. The asymptotics of the NLO µ are {u 0 , u 3 , u 2 , u 1 , u 5 } log u. First we introduce a new function This function has the infinite series of poles in the points i 2 + iZ. As the case of general integer conformal spin n has to be consistent with the known left-right symmetric case n = 0, it is natural to assume that the ansatz for the NLO µ-function is essentially the same as for n = 0, but relaxing the requirement that µ 3 is equal to µ 4 . So, we use the following ansatz Substituting this ansatz (B.14) into the Pµ-system equations, first, we get the unknown coefficient and obtain µ with the lower indices in the NLO order The µ-functions with the upper indices in the NLO are given by the formula where the matrix χ is given by (2.55). Usage of the equationP a = µ ab P b allows to partially fix the P-functions in the NNLO order. They are written in the formulas below 4 = π 2 Λ 2 3u 4 + 5Λ 2 u 6 . Also, from the obtained results one can notice that the Pand µ-functions possess the following symmetry P a (n, u) = χ ab P b (−n, u) , µ ab (n, u) = χ ac µ cd (−n, u)χ db . (B. 35) where the matrix χ is given by (2.55). In what follows we will assume that this symmetry is present in all orders of the perturbative expansion.

B.3 Passing to Qω-system
To proceed we need to use the system which is dual to Pµ-system -Qω-system. The equations of this system look as follows [16,17] We assume the scaling of the ω-functions is the same as in the case n = 0 considered in [22], which is motivated by the fact that the case of general n has to be consistent with the known case of n = 0. Thus, the function ω 13 scales as w −2 , ω 12 , ω 14 , ω 23 and ω 34 scale as w 0 and ω 24 scales as w 2 . For the lower indices we have ω 24 as w −2 , ω 12 , ω 14 , ω 23 and ω 34 as w 0 and ω 13 as w 2 . Let us remind the connection between µ-and ω-functions µ + ab = 1 2 Q ab|ij ω ij+ , µ ab+ = 1 2 Q ab|ij ω + ij . (B.37) With the obtained Pand µ-functions in the LO and NLO we can extract the functions Q ab|13 and Q ab|24 in the LO and NLO as well. Let us proceed with further calculations taking Q ab|13 only as the actions with Q ab|24 are completely the same but with the exchange n → −n.
One of the QQ-relations says and it allows to express Q 3|1 and Q 4|1 in terms of Q 1|1 , Q 2|1 and Q ab|13 . Using this fact and the equation we are able to derive the following second order Baxter equation for Q 1|1 in the LO We found two solutions q 1 and q 2 with the pure asymptotics u

D Values of the NNNLO intercept
Using the iterative procedure we managed to calculate the values of NNNLO intercept for a set of values of the conformal spin n. They are written below The parts of the solution are given by where we used the fact that (Γ U · const) = 0. Using the large u expansion of the kernel Γ U , it is possible to fix the coefficients c a|i where we used the fact that Γ D (const) = 0. Let us rewrite the equation (5.120) once more introducing a new designation for the LHS of (5.114) Z 1,2 ≡R 1,2 + R 1,2 + 2i where R a , a = 1, . . . , 4 are given by (5.114). After the substitution of the P-functions, Q-functions with two indices and the gluing matrix we see that the integral kernels Γ U and Γ D in (F.1) combine into the difference