Integrability of Conformal Fishnet Theory

We study integrability of fishnet-type Feynman graphs arising in planar four-dimensional bi-scalar chiral theory recently proposed in arXiv:1512.06704 as a special double scaling limit of gamma-deformed $\mathcal{N}=4$ SYM theory. We show that the transfer matrix"building"the fishnet graphs emerges from the $R-$matrix of non-compact conformal $SU(2,2)$ Heisenberg spin chain with spins belonging to principal series representations of the four-dimensional conformal group. We demonstrate explicitly a relationship between this integrable spin chain and the Quantum Spectral Curve (QSC) of $\mathcal{N}=4$ SYM. Using QSC and spin chain methods, we construct Baxter equation for $Q-$functions of the conformal spin chain needed for computation of the anomalous dimensions of operators of the type $\text{tr}(\phi_1^J)$ where $\phi_1$ is one of the two scalars of the theory. For $J=3$ we derive from QSC a quantization condition that fixes the relevant solution of Baxter equation. The scaling dimensions of the operators only receive contributions from wheel-like graphs. We develop integrability techniques to compute the divergent part of these graphs and use it to present the weak coupling expansion of dimensions to very high orders. Then we apply our exact equations to calculate the anomalous dimensions with $J=3$ to practically unlimited precision at any coupling. These equations also describe an infinite tower of local conformal operators all carrying the same charge $J=3$. The method should be applicable for any $J$ and, in principle, to any local operators of bi-scalar theory. We show that at strong coupling the scaling dimensions can be derived from semiclassical quantization of finite gap solutions describing an integrable system of noncompact $SU(2,2)$ spins. This bears similarities with the classical strings arising in the strongly coupled limit of $\mathcal{N}=4$ SYM.

1 Introduction: bi-scalar theory and results for "wheel" graphs The γ-deformed planar N = 4 SYM theory has an interesting double scaling limit [1] combining a vanishing 't Hooft coupling constant g 2 = N c g 2 Y M → 0 and the increasing twist parameters q j = e −iγj /2 → ∞ 1 in such a way that ξ j = gq i are kept fixed 2 . In the particular case ξ 1 = ξ 2 = 0 and ξ ≡ ξ 3 = 0, the limiting theory describes two complex scalar fields and is dubbed the bi-scalar χFT 4 [1]. The Lagrangian of this theory is given by where φ i = φ a i T a are complex scalar fields and T a are the generators of the SU (N c ) gauge group in the fundamental representation, normalized as tr(T a T b ) = δ ab . Notice that the quartic scalar interaction term in (1.1) is complex making the theory nonunitary. As we show below, this leads to a number of unusual properties of the bi-scalar χFT 4 .
The γ-deformed version of N = 4 SYM breaks the global P SU (2, 2|4) superconformal symmetry down to SU (2, 2) × U (1) 3 . In the bi-scalar theory (1.1) this symmetry is further reduced to SU (2, 2) × U (1) 2 . For the large majority of physical quantities, such as the multi-point correlators 3 or amplitudes, this theory shows a CFT behavior in the planar limit. This bi-scalar χFT 4 is shown to have a very limited set of planar Feynman graphs in the perturbative expansion of any physical correlator. The number of planar graphs at each loop order is not growing with the order. Moreover, at sufficiently large loop orders, these graphs have in the bulk the "fishnet" structure of a regular square lattice pointing on the explicit integrability of the model (1.1) in the planar limit [1]. This integrability is due to the fact noticed in [4] that these fishnet graphs define an integrable lattice model.
In this paper, we show that the bi-scalar model (1.1) represents a four-dimensional field-theoretical realization of the planar fishnet lattice model and identify the underlying integrable model as being a noncompact Heisenberg spin chain with spins belonging to infinite dimensional representation of the four-dimensional conformal group SU (2,2). This fact sheds a certain light on the origins of, still mysterious, integrability of the planar N = 4 SYM [5]. Similar noncompact Heisenberg spin chains have been previously encountered in the study of high-energy asymptotics in QCD [6,7]. We shall apply the technique developed in [8,9] to clarify the basic features of integrability of the fishnet graphs using the method of the Baxter Q−operator. In particular, we describe the Lax operators, transfer-matrices and construct the Baxter T Q equation that we later use to compute the anomalous dimensions of certain operators in the bi-scalar model (1.1). However, in order to extract the anomalous dimensions from the T Q relations we had to use an additional insight coming from the side of N = 4 SYM, where the problem of finding the spectrum is solved by the Quantum Spectral Curve (QSC) method [10][11][12][13].
The integrability of the bi-scalar χFT 4 , together with a specific, limited Feynman graph content of this theory provides us with a powerful method of computing new multi-loop massless four-dimensional Feynman integrals. In the paper [14] a large variety of such graphs, relevant for the computations of anomalous dimensions, is described and some of these graphs are computed by the help of integrability. Similar double scaling limit has been observed in [14] for the ABJ(M) model, for which the QSC is also known [15,16]. In the ABJ(M) theory the Feynman graphs are dominated in the bulk by the regular rectangular "fishnet" lattice structure. In the paper [17], a similar "fishnet" structure for the bi-scalar amplitudes and single-trace correlators in the theory (1.1) has been observed and their integrable structure in terms of an explicit Yangian symmetry has been established. The authors of [18] introduce a similar tri-scalar χFT 6 with φ 3 −type chiral interactions, realizing the hexagonal finsnet graphs of [4] and compute the simplest all-loop 2-point correlation function. This theory seems to define a genuine CFT 6 in planar limit for all local operators. In the paper [19] a single trace 4point correlation function given by a rectangular fishnet graph has been explicitly computed using the bootstrap methods for AdS/CFT integrability.
The operators in the theory (1.1) can be classified with respect to the irreducible representations of its global symmetry SU (2, 2)×U (1) 2 , characterized by the values of Cartan generators (∆, S,Ṡ|J 1 , J 2 ), where ∆ is the scaling dimension of the operator, the pair (S,Ṡ) defines its Lorentz spin and the U (1) charges J i count the difference between the total number of φ i and φ † i . In this paper, we study the scaling dimensions of a family of scalar single-trace operators of the following schematic form O J,n, = P 2 (∂)tr[φ J 1 φ n 2 (φ † 2 ) n ] + . . . , (1.2) where P 2 (∂) denotes 2 derivatives acting on scalar fields inside the trace with all Lorentz indices contracted. The dots stand for similar operators with the scalar fields φ 1 , φ 2 and φ † 2 exchanged inside the trace. The operators (1.2) belong to the representation (∆, 0, 0|J, 0) with ∆ = J + 2(n + ) for zero coupling ξ = 0. For nonzero coupling, the operators (1.2) with the same ∆ and J mix with each other and their scaling dimensions can be found by diagonalizing the corresponding mixing matrix. In the special case n = = 0 the operator (1.2) takes the form Similar operator can be also defined in N = 4 SYM in which case it is protected from quantum corrections and is known as the BMN vacuum operator. In the γ-deformed N = 4 SYM, the operator (1.3) is not protected and its scaling dimension starts to depend on the coupling. Furthermore, in the bi-scalar theory (1.1) the coupling dependent corrections to the scaling dimension ∆ only come from the wheel-type Feynman graphs [1] shown in Fig. 1. Each wheel contains J interaction vertices and the contribution to the scaling dimension of the wheel graph with M frames scales as ξ 2JM . For M = 2 and arbitrary J, the contribution of the double-wheel (with two frames) graph to the scaling dimension of the operator (1.3) in γ-deformed N = 4 SYM was found in [20] using the TBA-type computations and represented in [1] in terms of explicit multiple zeta values (MZV).
In this paper, we use integrability technique to compute the scaling dimensions of the operators (1.2) with J = 3 for an arbitrary coupling ξ. At zero coupling, the scaling dimensions of such operators take odd values ∆(0) = 3, 5, . . . . For ∆(0) = 3 there is only one operator tr(φ 3 1 ). For ∆(0) = 5 we can separate all operators into those involving derivatives, tr(φ 2 1 φ 1 ) and tr(φ 1 ∂ µ φ 1 ∂ µ φ 1 ), and those built from five scalar fields. Making use of the equations of motion in the theory (1.1) and discarding operators with total derivatives, we can express the former operators in terms of the latter. The relevant operators are of the following form tr(φ 3 1 φ 2 φ † 2 ), tr(φ 2 1 φ 2 φ 1 φ † 2 ), tr(φ 1 φ 2 φ 2 1 φ † 2 ) and tr(φ 2 φ 3 1 φ † 2 ). Examining their mixing matrix, we find that, due to a particular form of the interaction term in (1.1), these operators form the so called logarithmic multiplet, typical for a non-unitary theory [21,22]. As a result, for ∆(0) = 5 there are only two unprotected conformal operators. We determine their scaling dimension for an arbitrary coupling ξ. The two other operators form as logarithmic multiplet leading to a logarithmic factor in their correlation functions. For ∆(0) ≥ 7, the basis of the operators (1.2) contains operators with derivatives and the spectrum of their scaling dimensions has a more complicated structure. In particular, in virtue of integrability, the spectrum contains pairs of operators with coinciding scaling dimensions. We argue that for each odd ∆(0) ≥ 7 there are two operators that are not degenerate. We determine the scaling dimensions of such operators for an arbitrary coupling.
To compute the scaling dimensions of the operators (1.2), we employ the method of quantum spectral curve (QSC). This method has been previously used to compute the anomalous dimensions of local single-trace operators in planar N = 4 SYM theory [10,11] and it can be generalized to twisted version of this theory [12,13]. The twisted QSC method (tQSC) should be adopted for the twisted N = 4 SYM theory in the double scaling limit, which is the principal technical problem we are confronted to. A similar problem has been recently solved in [23] for the case of the cusped Wilson line in a very similar limit to the one leading to the bi-scalar model (1.1). We will adopt the most important elements of the method of [23] for our current, more complicated case.
The tQSC method consists of two ingredients: i) derivation of the Baxter equation for particular set of the operators (1.2), ii) derivation of the quantization condition which fixes the dependence of the parameters of this equation on the coupling constant ξ. This information is sufficient to determine the scaling dimension of the operators ∆ = ∆(ξ). We solve the problem i) by applying the double expansion with respect to the parameters g and 1/q 3 , of the coefficients of the most general Baxter equation derived within the QSC method in [24]. To check this result, we use the fact that the same Baxter equation describes the conformal, SO(1, 5) ∼ SU (2, 2) spin chain in principal series representation, for a given set of Cartan charges (∆, 0, 0) corresponding to the scalar operators (1.2). It is known to have a universal form of a fourth-order finite difference equation with the coefficients depending on the integrals of motion including ∆. For the sake of consistency, we also re-derive the Baxter equation directly from the equivalent conformal spin chain describing the dynamics of "fishnet" diagrams. Moreover, we show that the same equation can be defived by imposing a certain symmetry among its coefficients, inspired by algebraic properties of the conformal Lax operator [25,26] along with the u → −u symmetry and large u asymptotics of its 4 solutions. The problem ii), allowing us to fix the dependence of ∆ and the remaining coefficients of Baxter equation on the coupling constant, is solved by exploiting the explicit analyticity properties of the QSC equations [10][11][12][13], including the Riemann-Hilbert relations completing the general QSC Baxter equation, in the form proposed in [27]. The direct derivation of the same quantization conditions from the conformal SU (2, 2) spin chain is not yet available.
Let us summarize the main results of this paper. We show that for nondegenerate operators mentioned above the Baxter equation possesses the additional symmetry under the exchange of the spectral parameter u → −u. This symmetry fixes the values of all but two integrals of motion of the spin chain and leads to a remarkable factorization of the 4th order Baxter equation into two 2nd order finite-difference equations of the form where q 2 (u, m) and q 4 (u, m) denote special solutions to (1.4) which have the following large u asymptotic expansions (the reason to label them as 2 and 4 will be clear below) It is important here that the expansion on the right-hand side runs in powers of 1/u only and there is no admixture of O(u −∆/2+3/2 ) and O(u ∆/2−1/2 ) terms in the first and in the second relations, respectively. That is the reason why we refer to (1.6) as "pure" solutions.
Having constructed the functions q 2 (u, m) and q 4 (u, m), satisfying (1.4) and (1.6), we can use (1.5) to compute the scaling dimension ∆ for any value of the coupling constant ξ. This procedure can be carried out numerically, with practically unlimited precision following the general method [23,27].
As an example, we show in the Fig. 2 our results for the scaling dimension ∆ 3 (ξ) of the simplest operator (1.3) for J = 3. At weak coupling, ∆ 3 (ξ) receives corrections from the wheel diagrams shown in Fig. 3. As a function of the coupling constant, ∆ 3 − 2 decreases with ξ and turns from real to purely imaginary values at ξ = ξ with ξ 3 0.2, indicating that the weak coupling expansion has a finite radius of convergency. We find that (∆ 3 (ξ) − 2) 2 is a smooth real function of the coupling and, therefore, ∆ 3 (ξ) has a square-root singularity at ξ = ξ . For ξ > ξ , the scaling dimension takes the form ∆ 3 = 2 + id(ξ) where real valued function d(ξ) monotonically increases with the coupling and scales as d ∼ ξ 3/2 for ξ → ∞.
The relations (1.4) -(1.6) are invariant under ∆ → 4 − ∆. Therefore, for any solution ∆(ξ) to the quantizaton condition there should exist another one 4 − ∆(ξ). Then, the appearance of the imaginary part of ∆(ξ) can be interpreted as a result of the collision of two 'energy levels' ∆(ξ) and 4 − ∆(ξ) at ξ = ξ (see Fig. 2). Although the level crossing cannot happen in a unitary conformal field theory, it occurs in the theory (1.1) since it is not unitary.
At weak coupling the quantization conditions (1.5) can be solved analytically order-by-order in ξ 6 . The approach is algorithmic and is only limited by the computer power available. The first four Numerical results for the scaling dimension of the operator tr(φ 3 1 ) as a function of the coupling ξ 3 . We observe a "phase transition" at ξ 3 0.21 where the scaling dimension takes the value ∆ = 2 and becomes imaginary. This point defines the radius of convergency of the weak coupling expansion. The second branch, starting from ∆(0) = 1, arises due to the symmetry of the Baxter equation (1.4) terms of the weak coupling expansion of ∆ 3 − 3 are 4 where ζ i1,...,i k = n1>···>n k >0 1/(n i1 1 . . . n i k k ) are multiple zeta functions. Here the coefficients in front of ξ 6M give the residues at simple pole in dimensionally regularized Feynman integrals corresponding to the J = 3 wheel graphs with M = 1, 2, 3, 4 frames (see Fig. 3). The first two terms of (1.7) (up to double wrapping) coincide with the known results [28,29].
At strong coupling, for ξ 1, the scaling dimension ∆ 3 − 2 takes pure imaginary values that scale as ξ 3/2 . In this limit, the quantization conditions (1.5) can be solved using semiclassical methods. We find that the strong coupling expansion runs in powers of 1/ξ 3 and the first few terms are given by Notice that the expansion coefficients grow factorially indicating that the series is not Borel summable.
In a close analogy with the strong coupling expansion of the cusp anomalous dimension in planar N = 4 SYM [30], this suggests that ∆ 3 receives at strong coupling exponentially small, nonperturbative corrections of the form e −cξ 3 . The same technique can be applied to compute the scaling dimensions of the operators (1.2) with the same R−charge J = 3 and ∆(0) = 5, 7, . . . in a free theory. As explained above, for ∆(0) = 5 we encounter two unprotected operators. We present the numerical and analytic study of quantum corrections to their scaling dimensions both at weak and at strong coupling. We show that the resulting expressions have a different behaviour with the coupling constant as compared with ∆ 3 (see Fig. 2). Namely, they acquire an imaginary part at weak coupling and increase monotonically with the coupling. For ∆(0) ≥ 7 we compute the scaling dimensions of a pair of nondegenerate operators, one for each ∆(0). We demonstrate that for ∆(0) = 7, 11, 15, . . . and ∆(0) = 9, 13, 17, . . . the dependence of the scaling dimensions of the coupling follows the same pattern as ∆ 3 (ξ) and ∆ 5 (ξ), respectively (see Fig. 8 below). The limit of large length ξ for such operators will be also treated and the explicit formula for the dimension will be found in the form of a Bohr-Sommerfeld type quantisation condition.
2 SU (2, 2) picture for fishnet graphs To prepare the formalism for computing the wheel-graphs shown in Fig. 1, we demonstrate in this section that these graphs can be identified as transfer matrices of an integrable conformal SU (2, 2) spin chain.
To begin with, we consider the two-point correlation function of the operators (1.3) where the normalization constant d J (ξ) and the anomalous dimension γ J (ξ) depend on the coupling constant and the length L = J of the operator. As was shown [1], in the planar limit this correlation function receives contribution from globe-like graphs shown in Fig. 4. The same graphs can be viewed as fishnet graphs with periodic boundary conditions along horizontal (longitudinal) direction and all lines in the vertical (latitudinal) direction joined at the points x and 0, respectively, as depicted in 0 x Figure 4. A "wheel" Feynman graph with M frames and J spokes whose external legs have been joined at the point x (M = 5 and J = 9 in this picture).
Here the operator H J adds an extra row to fishnet graph and '•' denotes the convolution of the kernel (2.5). As we show in the next subsection, H J can be identified as a transfer matrix of the non-compact Heisenberg spin chain of length J with spins being the generators of the conformal group SU (2, 2). We expect that the anomalous dimension in (2.1) is different from zero. This means that the correlation function (2.2) has to develop ultraviolet divergences and, therefore, requires a regularization. If we introduced dimensional regularization with D = 4 − 2 , these divergences would appear as poles of (2.2) in 1/ . Notice that for arbitrary x i and y i , the function T J,M (x 1 , . . . , x J |y 1 , . . . , y J ) is well-defined in D = 4 dimensions and does not require any regularization. The divergences appear in (2.2) when we identify the points y 1 = · · · = y J and/or integrate over x i → x. In both cases they come from integration in (2.4) in the vicinity of the two points z i = 0 and z i = x and have a clear UV origin. To lowest order O(ξ 2J ), the dimensionally regularized integral in (2.2) has simple pole 1/ whose residue gives the anomalous dimension. At high orders, the integrals entering the two-point correlation function (2.2) have overlapping divergences and, as a consequence, they produce higher power of 1/ . The divergent part of (2.1) has the following form The question arises how to use a four-dimensional representation (2.2) to extract correctly the divergent part of the correlation function (2.8). To this end we notice that, in distinction from D(x), its logarithm has at a most simple pole ln D(x) ∼ −γ J (ξ)/ and, therefore, does not contain overlapping divergences 5 . The simple pole of ln D(x) originate from the two integration regions in which all internal vertices approach one of the external points, z i = x or z i = 0. Then, applying the dilatation operator to ln D(x), we can remove its divergent part and evaluated the anomalous dimension γ J (ξ) = −[(x∂ x )/2 + J] ln D(x) using the four-dimensional representation (2.2).
where G = J i=1 G i are given by the sum of differential operators acting on points We recall that general expressions for the generators of the conformal group depend on the parameters (∆, S,Ṡ) defining the scaling dimension and the Lorentz spin. The relations (2.10) correspond to a free scalar field representation of the conformal group with ∆ = 1 and S =Ṡ = 0.
To verify (2.9) we apply its both sides to a test function Φ(x) depending on the set of coordinates where we used a shorthand notation for dy = dy 1 . . . dy J and introduced subscript to indicate that G x acts on x variables. Here in the second relation we integrated by parts and introduced notation for the conjugated generators G † y . For instance, P † µ = −P µ and D † = −(x∂) − 3. As follows from (2.11), the conformal symmetry implies that the kernel of the operator H J has to satisfy the following relation It is then straightforward to check that the kernel (2.5) verifies this relation indeed.
Having identified the representation of the conformal group (2.10), we can now construct the Heisenberg spin chain with the spin generators given by (2.10) and establish the relation between H J , Eq. (2.6), and transfer matrices of this model. As a first step, we define the R−operator acting on the tensor product of two representations (2.10). Similar to (2.5) it can be realized as an integral operator [8] with Φ(x 1 , x 2 ) being a test function. The requirement for the operator R 12 (u) to satisfy the Yang-Baxter equation leads to a differential equation for the kernel R u (x 1 , x 2 |y 1 , y 2 ). Its solution is given by [31] where c(u) is a normalization factor Its value was fixed by imposing the so-called T −inversion relation R 12 (u)R 12 (−u) = 1, or equivalently It proves convenient to use a diagrammatic representation for the R−operator (2.13) and (2.14).
Representing each factor of the form 1/[(x − y) 2 ] α as a solid line connecting points x and y with index α attached to it, we can depict (2.14) as a rectangular graph shown in Fig. 5. Notice that the indices of all four lines depend on the spectral parameter. Let us examine (2.13) for u → 0. Putting u = − in (2.14) and (2.13) we find for → 0 where we replaced c( ) by its leading asymptotic behaviour. For the expression on the right-hand side to be different from zero, the integral should produce a double pole 1/ 2 . Indeed, making use of the identity Figure 5. Diagrammatic representation of the R−operator (2.14). Solid line with index α stands for 1/(x 2 ) α . The values of indices depend on the spectral parameter and are given by α+ = −1 − u, α− = 1 − u and β = 2 + u.
we find that integration in (2.17) yields R 12 (0)Φ(x 1 , x 2 ) = Φ(x 2 , x 1 ), so that R 12 (0) coincides with the permutation operator, The same property can be easily found from the diagrammatic representation of the R−operator (see Fig. 5). The operator R 12 (0) is described by the diagram in which two lines carry index β = 2. According to (2.18), such lines collapse into two points, x 1 = y 2 and x 2 = y 1 , leading to (2.19).
We can now use the R−operator (2.13) to construct the so-called fundamental transfer matrix where Φ is a test function and the kernel T J,u (x|y) is given by a J−fold integral over the auxiliary space Using diagrammatic form of the R−operator, we can represent this integral in the form of J rectangles glued together through common vertices in a pairwise manner as shown in have index 2 which allows us to apply the identity (2.18). For u = − and → 0 we find 24) or equivalently T J (0) coincides with the operator of the cyclic shift. The same relation can be also obtained using (2.19). For u = −1 + and → 0 we find in the similar manner where y J+1 = y 1 and the additional power of 1/ appears due to the fact that c(−1 + ) defined in (2.15) is regular for → 0. Substituting this relation into (2.21) we obtain the following result for the transfer matrix (2.26) Comparing this relation with (2.6) we conclude that the operator H J defines the leading asymptotic behaviour of the transfer matrix T J (−1 + ) for → 0. Since the transfer matrix commutes with the local integrals of motion for an arbitrary spectral parameter, the same is true for the operator H J .
Having established integrability of the operator H J , we can now turn to computing the correlation function (2.2) and (2.7). A direct calculation of (2.2) would require diagonalizing the operator H J with subsequent summation over the whole spectrum of its eigenstates. Since the underlying Heisenberg spin chain is defined for an infinite-dimensional representation of the conformal group, this proves to be an extremely nontrivial task. In the next section, we present another approach based on the Baxter T Q equation, combined with the QSC formalism, that allows us to avoid these difficulties and obtain a compact representation for the anomalous dimension γ J (ξ).

The Baxter T Q equation
We have shown in the previous section that the corrections to the scaling dimension of the scalar operators (1.3) can be described in terms of the Heisenberg spin chain with the spin operators being the generators of the conformal group. In this section, we derive the Baxter T Q equation in this integrable model for the states corresponding to a special class of scalar operators defined in (1.2).
The Baxter equation is a fourth-order finite difference equation for the function q(u) with the coefficients given by transfer matrices of SU (2, 2) spin chain in specific, antisymmetric representations. We argue below that, for some states with the charges (∆, 0, 0|J, 0), the general form of this equation can be determined, up to a few constants (to be fixed by additional quantization conditions, see section 5), from the known asymptotic behaviour of four solutions for the Baxter function q(u) at infinity combined with the symmetry under parity transformation u → −u. We partially justify these assumptions in Appendix A by making use of the Lax operator formalism.
It was shown in the paper [26] (see eq.(5.5) there) that the Baxter TQ-equation for SL(4) takes a standard form of a linear 4th order finite difference equation (A.24) on Q-functions, with coefficients being the fundamental transfer-matrices of the SL(4) spin chain with spins in principal series representations. In Appendix A we use this formula and u ↔ −u symmetry arguments for the states with the charges (∆, 0, 0|J, 0) to bring the Baxter equation to the following symmetric form 1) where A(u) = u J is completely fixed by the choice of the representation (2.10), whereas B(u) and u J C(u) are polynomials in u of degree J and 2J, respectively. The coefficients A(u) B(u) and C(u) are explicitly calculated in Appendix A in terms of global charges of the state and they contain a certain number of state-dependent constants, to be defined later by quantization conditions proper to our problem. To fix the form of polynomials A(u), B(u), C(u) it appears to be enough to impose two additional conditions: if q(u) is a real solution to (3.1) then q(−u) should also be a solution. This gives the parity condition on the coefficient functions: Secondly, for large u, the solution to (3.1) should have asymptotic behavior q ∼ u δ with the exponent δ taking the following values (see next section or the twisted QSC formalism of [13]) These conditions fix the functions B(u) and C(u) to be Baxter equation for J = 2 where all constants are fixed in terms of α = (∆ − 2) 2 . We notice that this equation factorizes as where D = e i∂u is the shift operator.
Baxter equation for J = 3 where α = (∆ − 2) 2 and m is arbitrary. Remarkably, this equation also factorizes It depends on 3 extra constants, b, c 1 , c 2 , apart from α = (∆ − 2) 2 , to be fixed from quantization conditions, yet to be derived. The Baxter equations (3.5), (3.7) and (3.9) match perfectly the equations derived in detail for particular cases J = 2, 3, 4 in Appendix A from the Lax operator formalism and they will be confirmed in section 4 from the QSC formalism.
4 Baxter equations for bi-scalar χFT 4 from the double scaling limit of QSC In this section we use an alternative integrability based approach which originates from N = 4 SYM. The integrability structure for the spectrum in this theory is very well studied and is given in terms of so-called Quantum Spectral Curve (QSC) [10,11]. We will use that fact that χFT 4 can be obtained as a limit of a twisted version of N = 4 SYM in order to get further clues about the spectrum of anomalous dimensions and integrability structure of χFT 4 . In particular in this section we demonstrate how the Baxter equations obtained from the integrability of the Feynman graphs can be obtained directly from QSC establishing an important link between these two seemingly different integrability based approaches. We will see that QSC in addition to the Baxter equation also provides an essential missing ingredient -the quantization condition for the spectrum, which we were not able to obtain from the Feynman graphs. The quantization condition is discussed in Section 5. In this section we describe how the Baxter equation arises from QSC approach in the double scaling limit where the twisted SYM reduces to the scalar χFT 4 .

QSC generalities
The description of QSC and the notations closely follow [11][12][13]. In formal terms, the QSC is represented by a Grassmannian consisting of 2 8 so-called Q-functions of spectral parameter u related to each other by Plücker QQ-relations, in such a way that 8 of them are enough to parameterize all other. General Q-functions can have quite complicated analytical structure as functions of the spectral parameter. In practice it is sufficient to restrict ourselves to a subset of functions with simplest analytic properties, namely They are related to each other by the following set of equations: and similar relations with all upper/lower indices lowered/raised. Here and below we use a common notation Apart from these algebraic relations, there are also analytic constraints on the structure of cuts of Q-functions in the complex plane, which, in particular, express their monodromies around the branch points through other Q-functions.
The basic analytic structure is the following: P a and P a have one cut on the real axis, from −2g to 2g, as shown in Fig. 7. This cut can be resolved by a transfromation to Zhukovsky variable x(u), defined as x(u) + 1/x(u) = u/g, |x(u)| > 1. Thus it will be convenient to parametrize P in x instead of u.
Functions Q i and Q i have infinite ladder of Zhukovsky cuts: the original one, from −2g to 2g and its copies, shifted by an integer number of i into the lower half plane, as shown in Fig. 7. For any function of spectral parameter Q(u) we denote byQ(u) the analytic continuation of Q(u) under the Zhukovsky cut on the real axis. It is sufficient to impose the gluing condition on Q i at the cut which schematically can be written asQ i (−u) = H i j Q j (u) 7 , to close the system of equations. This means that imposing the above relation will give us a discrete set of isolated solutions each corresponding to a certain state of the theory.
In order to identify a particular state we have first to know its quantum numbers, which are hidden in the large u asymptotic of Q i and P a as we discuss below. Figure 7. Left: Analytic structure of Pa(u) and P a (u): P(u) has one cut on the first sheet (left). Analytical continuation through this cut leads to the second sheet corresponding toP(u) (right), which has an infinite ladder of cuts Right: Analytic structure of Qi(u) and Q i (u): Q(u) has an infinite ladder of cuts in the lower half-plane (left). Analytical continuation through the cut on the real axis leads to a sheet corresponding tõ Q(u), which has an infinite ladder of cuts in the upper half-plane (right).
Asymptotics and quantum numbers. Asymptotic behavoir of Q-functions at large u is determined by the twists and by the quantum numbers of the particular state of N = 4 SYM we are studying. Here we will consider a particular kind of state: BMN vacuum tr(φ J 1 ), where φ 1 is one of three complex scalars of N = 4 SYM theory 8 . In the notation of [13] its quantum numbers are {J, 0, 0|∆, 0, 0} . (4.5) In the untwisted theory tr(φ J 1 ) would have been protected BPS state, but in the presence of twists it has a nontrivial scaling dimenstion ∆ J (g) which we will be computing. This scaling dimension was computed at weak coupling up to g 4J−2 (two wrappings) in [20] in terms of infinite double sums 6 We used here the left-right symmetry relations (C.2) and already put at this stage κ =κ. 7 In the unitary theory like N = 4 SYM one can use the complex conjugation instead, which is applicable for all operators and not only for the parity symmetric ones. Also, for this gluing condition to be of this simple form that is crucial to choose Q i to be of a "pure" form, meaning that their large u asymptotic are of the form u −ν i ∞ i=0 a i u i to all orders in 1/u. In general this condition is not sufficient to uniquely determine Q i as ν i for different i could differ by an integer, allowing for Q i to mix. There are several ways to deal with this problem, one possibility is to introduce a twist in AdS 5 . Another possibility, which we use in this paper, is to make use of the u → −u symmetry, applicable for some states, and keep only even powers u 2i which also leads to the correct non-ambiguous gluing. Finally, it was noticed in [32] that the conditions (4.4) are over-defined and it is usually sufficient to impose only one of them. 8 We use here notation φ i for complex scalars, instead of the familiar {X, Y, Z} since two of them appear with a different notation in bi-scalar model (1.1), say Z = φ 1 , Y = φ 2 . and integrals and it was brought to the standard explicit MZV form in [1]. For the case of our current interest, the bi-scalar model, the single wrapped graph at any J was computed in [28] and the double-wrapped graph at J = 3 was computed in [29]. We use these results to verify our computation.
In the general case of twisted QSC, the asymptotics of one-indexed functions are given by [13] In a particular case of BMN vacuum with charges (4.5) the powersλ a ,λ * a ,ν i ,ν * i are determined by the quantum numbers of the stateλ and the twists are chosen to be x a = κ J , κ −J ,κ J ,κ −J . They are related to the twist parameters introduced at the beginning of introduction as whereas the dependence on q 1 is absent for BMN vacuum state. We can plug the asymptotics (4.6) into the equations (4.2) and obtain a set of constraints for combinations A a A a and B i B i (no summation). After the fixing the rescaling symmetry the solution for A a can be chosen to be [13] .
In a similar way for B i : Since, as we discussed, the P a and P a functions has only one cut, they admit the following Laurent expansion The constants b a,n and b a,n contain all the information about the state. In order to constrain them we have to find Q i and Q i for the given P a and P a and impose the analyticity condition (4.4) on the cut [−2g, 2g]. The Fourth-order Baxter equation is an efficient way of doing this.
Fourth-order Baxter equation. Instead of using the chain of QQ-relations (4.2) one can equivalently deduce the Q i using the given P a and P a functions from the "Fourth-order Baxter equation". More precisely, four one-indexed functions Q j are the 4 linear independent solutions of the following 4th order finite-difference equation [24]: where D i are determinants of matrices with elements of form P [k] a given in the appendix B.

Double-scaling limit
The limit we have to take to obtain the χFT 4 is g → 0 and s ≡ √ κκ → ∞ with ξ = gs fixed. First obvious thing which happens in this limit the branch points of the Zhukovsky cuts collapse into a point producing poles. Furthermore, the coefficients b a,n and b a,n in (4.12) typically scale as g 2n and thus the infinite series truncates. In Appendix C we argue that for the BMN state the following scaling in g has to be imposed: (4.14) where f i , g i as series in x(u): This follows the argument similar to [33]. We see that the terms with heigher dergree of 1/x n get more and more suppressed. However, we still have to keep a few first term simply because the Baxter equation (4.13) contains P a with shifts u → u + in, n = −2, . . . , 2. Due to the twists tending to infinity along the complex axes the factors x iu a in P a we could get an enhancement of these suppressed with g terms. At the end we have to plug the truncated expressions for P a , which become rational functions with poles at u = 0 in the g → 0 limit, into the equation (4.13) we obtain a finite difference equation with rational coefficients. We worked out explicitly the J = 2, 3, 4 cases in the Appendix C. Below we present the results.
Baxter equation for twist 2. For J = 2 the only unfixed coefficient left is ∆ and the Baxter equation (4.13) coincides in the double scaling limit with (3.6). The expression on the left-hand side of (3.6) is given by the product of two finite difference operators separated by the factor of u 2 . It is easy to verify that (3.6) stays invariant under the exchange of these two operators. As a consequence, the four solutions to (3.6) have to satisfy 2nd order finite difference equations, e.g. 16) and the second equation is obtained by replacing ∆ → 4 − ∆. By matching the asymptotics (4.6), we find that Q 2 and Q 3 satisfy (4.16) whereas Q 1 and Q 4 satisfy the second equation. The 2nd order Baxter equation (4.16) has been previously studied in [9,24] and its explicit solutions have been found in terms of hypergeometric 3 F 2 −functions 17) and the second solution is given by q 0 (−u). Using these results, we construct their linear combinations These expressions are solutions to the Baxter equation (4.13) with the correct pole structure. Moreover, they are "pure" solutions in the sense that their expansion at infinity has a form u α 1 + c 1 /u 2 + . . . . The two remaining solutions Q 4 and Q 1 are obtained from Q 2 and Q 3 by replacing ∆ with 4 − ∆ (and changing B i accordingly).
Baxter equation for twist 3. For J = 3 the equation is a bit more complicated. After using the identification the equation (4.13) reduces in the double scaling limit to the Baxter equation (3.7) where α = (∆−2) 2 and m are two yet unfixed parameters. Its factorized form is given by (3.8). In other words again we reproduced precisely the equation obtained from the integrability of the Feynman graphs! In the next section where we will derive the quantization conditions for fixing m and, finally, ∆(ξ). Thus we will extract the all-loop anomalous dimension of Trφ 3 1 operator, and of some related operators. It can be also easily checked that the equation (3.7) is invariant w.r.t. the change m → −m. That is why we can find all four solutions from a much simple 2nd order Baxter equations as another pair of solutions can be obtained by replacing m → −m.
Baxter equation for twist 4. For J = 4 with the identification (4.19) we obtained again the same Baxter equation as in (3.9), for which we have not been able to show the factorization property. We postpone further investigation of this state for a future publication.
In conclusion, the derivation in this section of Baxter equations from QSC formalism in the double scaling for J = 2, 3, 4 confirms the Baxter equation for arbitrary J given at the end of the previous section, obtained from the group-theoretical considerations, in usual assumption of universality of Baxter equations (independence on auxiliary state representation). Now we have to impose additional quantization conditions on possible solutions of Baxter equation, which we do in the next section, so far only for J = 3 case, from the double scaling limit of QSC formalism.

Quantization condition for the Baxter equation
In this section we focus on the operators (1.2) with the R−charge J = 3 whose scaling dimensions ∆ are described by the Baxter equation (4.20). The method presented here should also apply for general J and even to general operators in the bi-scalar theory (1.1).
The Baxter equation (4.20) does not single out a particular value of the parameters ∆ and m. Furthermore, it does not depend on the coupling constant ξ. The goal of this section is to use the underlying QSC description of the full N = 4 theory to derive the quantization conditions which fix the dependence ∆(ξ) and m(ξ). A priori, these conditions can be derived from the first principles using integrability of the fish-net diagrams [1,4]. However, this is beyond the scope of the current paper.
Let us first discuss analytic properties of the solutions of the Baxter equation (4.20). To build a solution we start from u with large positive imaginary part, for which the finite difference equation (4.20) can be replaced by an ordinary differential equation. For arbitrary ∆ and m, it has two linear independent power-like solutions that we can choose to be "pure solutions" For arbitrary non-integer ∆ these solutions are well defined from the requirement that q 2 (u, m) does not contain O(u −∆/2+3/2 ) terms and q 4 (u, m) does not contain O(u ∆/2−1/2 ) terms. The expansion coefficients a i and b i can be fixed from (4.20), e.g.
Using (5.1) we can apply the finite difference equation (4.20) to recursively decrease the imaginary part of u in integer steps. Since the coefficients of the Baxter equation (4.20) are analytic for Im u > 0 and have a 3rd order pole at u = 0, the solutions constructed in this way are also analytical in the upper half-plane and have generically a 3rd order pole at the points u = −in with n = 1, 2, 3, . . . . We will discuss in more details how to build these solutions in the section 6.
Having constructed solutions to the Baxter equation (4.20), we can now identify four Q−functions defined in (4.1). In the double scaling limit we have It is straightforward to check that these expressions verify the defining relations (4.6) and their normalization is fixed by (4.10) (we recall that s = √ κκ and κ,κ → ∞ in the double scaling limit). To simplify the calculation, we have assumed Q−functions to be even functions of m. We can relax this condition and carry out the calculation in the general case. This will not affect the final result for ∆(ξ) but will significantly complicate some intermediate expressions. At the same time that is clear that the m → −m symmetry of the system has to be reflected in the Q i so instead of deriving the expressions for Q 2 and Q 4 we fix the proportion of q 4 (u, m) and q 4 (u, −m) from the symmetry requirement from the beginning.
The quantization condition for ∆(ξ) arises from the requirement for the functions (5.2) to have correct analytic properties in twisted N = 4 SYM presented at the beginning of section 4. We recall that for a finite value of the twist the functions Q i have the cut [−2g, 2g]. Going under the cut, these functions have to satisfy the gluing conditions [32]. In the double-scaling limit, the cut [−2g, 2g] shrinks into a point and the relations (4.4) lead to nontrivial constraints for the functions Q i (u) near the origin [32].
Another condition comes from the symmetry of the Baxter equation (4.20) under parity transformation, u → −u and m → −m. Indeed, the functions q i (u, m) and q i (−u, −m) satisfy (4.20), implying that the functions Q i (−u) can be expanded over the basis of the solutions (5.2) with some periodic coefficients: where Ω i j (u) has the following i-periodicity property 9 Similar relation should also hold under the cutQ i (−u) =Ω j i (u)Q j (u). As we show in the Appendix B, the discontinuity of Ω i j (u) across the cut, ∆Ω j i ≡Ω j i − Ω j i , satisfies the following exact Riemann-Hilbert equation (valid for any value of the twist parameter) In the double-scaling limit, this relation allows us to compute ∆Ω j i (u) for u → 0 in terms of solutions to the Baxter equation, q 2 (u, ±m) and q 4 (u, ±m), evaluated at the origin. Indeed, for u → 0, substituting (5.2) and (4.4) into (5.5) we find after some algebra 10 where s = √ κκ and the notation was introduced for with compact notations We conclude from (5.6) that ∆Ω j i vanishes near the origin in the double scaling limit as ∆Ω j i ∼ u 3 . In the next section, we will find Ω j i (u) independently for an arbitrary coupling from the Baxter equation. Comparing the results of these two calculations we will be able to determine the quantization condition for ∆ and m.

5.1
Extracting Ω j i from the Baxter equation As was mentioned in the beginning of this section, the Baxter equation (4.20) is invariant under simultaneous change u → −u and m → −m. Following the same logic that led (5.3), its solutions have to satisfy the relation where σ l k (u) are some nontrivial i-periodic functions. We can use the Baxter equation to find the leading behavior of σ j k (u) for u → 0. We recall that solutions to the Baxter equatin q j (−u, −m) have 3rd order poles at u = in (with n = 1, 2, . . . ) whereas q k (u, m) are regular at these points. Together with periodicity condition σ l k (u+ i) = σ l k (u), this implies that σ l k (u) must have 3rd order pole at the origin. Since the functions q j (−u, −m) and q k (u, m) are regular for u → 0, it follows from (5.9) that the residue at this pole have to satisfy lim Using the Baxter equation (5.1) we get for u → 0 Replacing u → −u + i in (5.9) and matching it into the last relation we obtain The relations (5.10) and (5.12) viewed as a linear system of equations on u 3 σ l k (u) for u → 0 lead to We notice that the expression in the denominator coincides with the Wronskian of the finite difference equation (4.20). As such, it should not depend on u and can be computed using asymptotic behavior of functions at infinity (5.1) In this way, we finally obtain We can now combine the relations (5.9) and (5.2) to establish linear relations between the functions Q i (−u) to Q j (u). Using the definition (5.3) we obtain from (5.15) the leading asymptotic behavior of Ω i j (u) at the origin In the next section we compare this relation with the expression for the discontinuity of Ω i j given by (5.6) to obtain the constraints on the parameters of the Baxter equation m and ∆.

Quantization condition from gluing
We demonstrated in the previous subsection that the matrix Ω i j (u) has 3 rd order pole at u = 0. We would like to stress that the relation (5.16) holds in the double scaling limit, for ξ = gs fixed with s → ∞ and g → 0. Due to our conventions different Q i scale differently with g and as a consequence different components of Ω i j scale differently in our limit. This is totally due to the choice of our normalization. We introduce γ ij so that Ω i j ∼ s γij . As explained in section 4, at finite coupling g the only singularities Ω i j (u) could have at finite u are due to the branch cuts of the Zhukovsky variables x(u) locate at u = ±2g. As a consequence, it can be represented in the form where 4 × 4 matrices f (g, u) and h(g, u) are regular around the origin and each term in their small g expansion should be regular as well.
The poles of Ω i j (u) at the origin can only appear as an effect of expansion of the Zhukovsky cut At the same time, computing discontinuity of Ω i j (u) across the cut [−2g, 2g] we find from (5.17) for g → 0 According to (5.6), ∆Ω i j ∼ u 3 for u → 0 in the double scaling limit. Then, it follows from the last relation that the series expansion of f (g, u) in g and u must be of the form where we assume that f (g, u) scales as g −n . The expansion coefficients depend on ξ and ∆ (but not of g). Substituting the last relation into (5.17) and taking into account (5.18) we get where we use a shorthand notation for a matrix s −γ Ω ≡ s −γij Ω j i (u) and dots denote terms regular for u → 0 and/or suppressed by powers of g 2 . Notice that regular function h(g, u) does not contribute here to the leading, singular terms.
We expect that in the double scaling limit the matrix (5.21) should take the expected form (5.16). In particular, it should scale as 1/u 3 for u → 0. Note that l.h.s. of (5.21) scales as g 0 . We have to deduce the value of n. If we assume that n = 4 the second term in (5.21) become negligible and we will not be able to reproduce the 1/u 3 singularity. If we take n = 6 we get a contradiction as the l.h.s. scales as g 0 whereas in the r.h.s. we get a term 1/g 2 . The only way to avoid this problem is to set f 0,2 + f 1,0 = 0. Using this relation and setting n = 6 we get for the residue of Ω(u) at 3rd order pole at u = 0 lim where in the second relation we took into account (5.19) and (5.20). This equation provides us with a set of nontrivial relations between the entries of matrices (5.6) and (5.16). In particular, we notice that the matrix element ∆Ω 2 1 is zero implying that Ω 2 1 should vanish. This leads to the quantization In the similar manner, the vanishing of Ω 3 1 implies ∆Ω 3 1 = 0 leading to q 2 4 − q 2 4 β 2 1 + q 2 2 −q 2 2 = 0. Together with (5.23) this gives Finally imposing the relation (5.22) for Ω 1 1 and ∆Ω 1 1 we find that the parameter m 2 is related to the coupling constant in the bi-scalar theory (1.1) It is straightforward to verify that all remaining conditions (5.22) are satisfied automatically! The quantization condition (5.23) fixes the dependence of ∆ on m. Together with (5.25) this allows us to find the spectrum of anomalous dimensions of the states/operators (1.2) with charges (∆(ξ), 0, 0|J, 0), obeying the parity invariance. 11 In the next sections we will explore these quantization conditions along with the Baxter equation (4.20) to study the dimensions of the underlying operators with charge J = 3 numerically, as well as in strong and weak coupling approximations.

Numerical solution
In this section, we describe the numerical solution to Baxter equation (4.20) supplemented with the quantization conditions (5.23) and (5.25). The numerical method for solving the whole QSC was developed in [27] and it can be adopted for our much simpler case (see also [23]). One first constructs analytically the solution to the Baxter equation (4.20) for large values of the spectral parameter u where one can simply take the asymptotic series (5.1) and find the coefficients of the expansion by 11 We remind that the full superconformal symmetry of N = 4 SYM is broken by γ-deformation to P SU (2, 2|4) → SU (2, 2) × U (1) 3 and an arbitrary state is still characterized by Cartan charges (∆, S 1 , S 2 |J 1 , J 2 , J 3 ). For the bi-scalar theory under consideration the remaining symmetry is SU (2, 2) × U (1) 2 (only two complex scalars out of three left) we label the states/operators with Cartan charges (∆, S 1 , S 2 |J 1 , J 2 ). For our particular BMN-state and its analogs considered further we take (∆, 0, 0|J, 0). plugging it into the equation and expanding for large u. We kept around 12 first orders which allowed us to gave extremely accurate approximation for q 2 (u, m) and q 4 (u, m) for Im u > 100. After that we used the exact Baxter equation (4.20) recursively to decrease the imaginary part of u until we reach u = 0. After that we define the function which must be zero for the physical values of ∆ and m 2 = −ξ 6 . One can simply use the FindRoot function in Mathematica to find ∆(ξ) determined by the condition F (∆(ξ), ξ) = 0.
In this way we obtained the results for the scaling dimension of the operator tr(φ 3 1 ) shown in Fig. 2. We see that the dimension approaches ∆ = 2, where it collides at ξ 3 0.21 with the "mirror" solution obtained by ∆ → 4 − ∆. After that the two dimensions stay in the plane Re ∆ = 2, while their imaginary parts increase at large ξ as ∆ ∼ ξ 3/2 . In the next section we describe this strong coupling asymptotics analytically. We also found that the quantization condition (5.23) also describes other states with Cartan charges (∆, 0, 0|L, 0). As we first found numerically (and then confirmed analytically, see the next section) the quantization condition F (∆, ξ) = 0 is satisfied for several values of ∆! At zero coupling, these extra solutions all start at odd integer values of ∆. For ∆(0) = 5+4n , n = 0, 1, 2, . . . the scaling dimension ∆(ξ) become complex for arbitrary small ξ > 0 whereas solutions with ∆(0) = 3 + 4n , n = 1, 2, . . . are real for small ξ. Our numerical analysis suggests that, similar to the state with ∆(0) = 3, all solutions with ∆(0) = 3 + 4n reach the value ∆ = 2 and then acquire an imaginary part.
Finally, using the high precision numerics we extracted the expansion coefficients of the weak coupling expansion of ∆(ξ) for the two lowest lying states We use these results in the next section to verify our analytic expressions.

Weak coupling solution
In this section, we describe the method for finding analytical solutions to the Baxter equation (4.20) with quantization condition (5.23) at weak coupling. We keep the presentation short since it goes along the same lines as in [23,32].

Perturbative solution of the Baxter equation
At the first step we have to find two linearly independent solutions to the Baxter equation (4.20) with the parameter m satisfying (5.25). At weak coupling, the solutions to (4.20) can be constructed perturbatively in powers of m. To lowest order, for m = 0 and L ≡ ∆(0) an odd integer, the Baxter equation (4.20) reduces to that for the SL(2) spin chain of length 2. As such, it has a solution (4.17) which for odd L reduces to a polynomial. The second solution then can be deduced from the Wronskian relation (5.14). It is a meromorphic function of u with the second-order poles located in the lower half-plane where P n and Q n are polynomials of degree n and the notation was introduced for the special function η s1,...,s k (u) with an appropriate pole structure 12 [34] η s1,...,s k (u) = In this way we find that the solutions to (4.20) satisfying (5.1) are given to the leading order in m = iξ 3 by L = 3 : q I = u , q II = 1 , To incorporate corrections in m, we use the ansatz ∆ = L+ k m k ∆ (k) and look for solutions to (4.20) in the form q I (u)+δq I (u) with δq I (u) = c 1 (u)q I (u)+c 2 (u)q II (u) and similar for q II (u)+δq II (u). This leads to the system of first-order finite difference equations for the coefficient functions c 1 (u) and c 2 (u) which can be solved order-by-order in m in terms of the functions (7.2) with polynomial coefficients. Let us consider the state with L = 3. Numerical solution (6.2) suggests that corrections to ∆ run in powers of m 2 = −ξ 6 . Using the ansatz we compute corrections to the solutions (7.3) At the next stage we have to find particular linear combinations of these functions, q 2 (u) and q 4 (u), which have the correct asymptotic behavior (5.1) at infinity. As follows from (5.1), the leading correction to q 2 (u) and q 4 (u) at large u should scale as δq 2 = − 1 2 (m 2 δ)u ln u and δq 4 = 1 2 (m 2 δ) ln u, respectively. Expanding q I (u) and q II (u) at large u and matching the coefficients we find that solutions to the Baxter equation satisfying (5.1) are given by To analyze the quantization condition (5.23), we have to evaluate these expressions at the origin. Making use of the identity η s1,s2,...,s k (i) = (−i) s1+s2+···+s k ζ s1,s2,...,s k , as well as relations between the multiple zeta values, e.g. ζ 2,2,2 = π 6 /5040, we arrive at where in our conventionsζ 1 = iη 1 (i) = γ is the Euler constant andζ 1,2 = i 3 η 1,2 (i) = γπ 2 6 − 2ζ 3 . We use these relations to find The quantization condition (5.23) yields δ = −12ζ 3 and, together with (7.4), fixes the dependence of the scaling dimension on the coupling constant.
With a help of Mathematica we pushed the calculation of (7.4) up to the order m 12 and arrived at the weak coupling expansion of ∆(ξ) given by (1.7) in the Introduction.
The above analysis can be repeated for the states with L = 5, 7, 9. We present below the results for weak coupling expansion of scaling dimensions of these states.
States with L = 5: Solving the quantization conditions (5.23), we found two states with ∆ = 5 at zero coupling. In distinction from (7.4), their scaling dimensions at weak coupling run in powers of m, or equivalently in powers of iξ 3 , and are complex conjugated to each other: States with L = 7: As in the previous case, we found two states. Similar to the state with L = 3, their scaling dimensions are real at weak coupling and have an expansion in powers of ξ 6 We verified that these relations are in perfect agreement with the numerical consideration at weak coupling. The strong coupling expansion of the scaling dimensions is discussed in section 8.

Logarithmic multiplet
As was mentioned above, the scaling dimensions of operators with bare dimension ∆(0) = 5, 9, . . . have expansion in powers of ξ 3 and not in powers of ξ 6 as it happens for operators with bare dimension ∆(0) = 3, 7, . . . . In this section we show that this property reflects a very unusual feature of biscalar χFT 4 theory first noticed by Joao Caetano [22]: the operators with ∆(0) = 5, 9, . . . behave as conformal primary operators in a logarithmic four-dimensional conformal field theory.
To simplify the consideration, we examine the simplest case of operators with bare dimension ∆(0) = 5 and the R−charge J = 3. As was mentioned in the Introduction, such operators can be obtained from tr(φ 3 1 ) operator by inserting a pair of scalar fields φ 2 and φ † 2 inside the trace and by dressing scalar fields with derivatives, tr(φ 2 1 φ 1 ) and tr(φ 1 ∂ µ φ 1 ∂ µ φ 1 ). In the latter case, we can use the equations of motion in the theory (1.1) to show that the operators with derivatives can be expressed in terms of the former operators as well as the conformal descendant operator tr(φ 3 1 ). This allows us to define the basis of dimension−5 operators At quantum level, these operators mix with each other and their anomalous dimensions can be found by diagonalizing the corresponding mixing matrix V ij where µ is an ultraviolet cut-off and the matrix elements To the lowest order in the coupling, the quartic scalar interaction vertex in (1.1) induces the following transitions The corresponding Feynman diagrams are shown in Fig. 9. Applying these rules we find that the operators (7.14) mix as follows where the factor of 2 is due to the fact that the first two transitions in (7.16) yield the same operator. The last transition in (7.16) produces the additional mixing between the operators O 3 and O 2 . Notice that the operator O 4 is not affected by the transitions (7.16). As follows from (7.17), the mixing matrix for the operators (7.14) takes the following form to leading order in ξ 2 Here ξ 2 γ 1 and ξ 4 γ 2 describe the first two and the last transitions in (7.16), respectively. The terms appearing as O(ξ 4 ) and O(ξ 6 ) in the matrix do not contribute to anomalous dimensions so we ignore them in what follows. Using the dimensional regularization with D = 4 − 2 , we can find them as the residue at a simple pole 1/ of diagrams shown in Fig. 9 γ 1 = −2 , γ 2 = 1 . Figure 9. Mixing of O3 = tr(φ † 2 φ1φ2φ 2 1 ) with the operators O4 and O2 to the lowest order in the coupling. Three diagrams correspond to three transitions defined in (7.16). Solid line denotes field φ1, dashed line represents φ2 and φ † 2 depending on the direction of the arrow.
Notice that γ 1 and γ 2 have an opposite sign. The mixing matrix (7.18) has a number of unusual properties. In unitarity conformal field theory this matrix has to be Hermitian. Since χFT 4 theory is not unitary, the matrix (7.18) does not have this property. This implies that its eigenvalues are, in general, complex valued functions of the coupling ξ. Indeed, the scaling dimensions (7.11) develop imaginary part at weak coupling.
Secondly, the matrix (7.18) has rank 3 and it can be brought to the Jordan canonical form by a similarity transformation If the matrix J were diagonalizable, its eigenspectrum would define four different conformal operators. Since J contains 2 × 2 Jordan block, the situation is more complex. Namely, we can use the lower diagonal 2 × 2 block of the matrix J to define two conformal operators ξ 2 O 2 ∓ 2iξO 3 − 4O 4 with the anomalous dimension γ ± = ±2i ξ 3 + O(ξ 6 ) . (7.21) Notice that this expression scales as O(ξ 3 ), in agreement with (7.11). The Jordan block of J describes the mixing matrix for the pair of the operators 16O 4 and O 1 + 4O 3 /ξ 2 that we denote as A and B. The form of this block is fixed by the interaction term in (1.1) and is protected from quantum corrections. The pair of the operators A and B belongs to the same conformal multiplet with the conformal weight ∆ = 5, a phenomenon typical for logarithmic conformal field theories [35]. To show this, we define the set of operators conjugated to (7.14) . . , 4. The reason why we label these operators in such a way is that they satisfy the same evolution equation (7.15) as operators O i with the mixing matrix given by (7.18). As a consequence, two of the operators have the anomalous dimension (7.21) and the two remaining onesĀ = 16Ō 4 andB =Ō 1 + 4Ō 3 /ξ 2 have a mixing matrix given by 2 × 2 Jordan cell. Computing the correlation functions of the operators A(x), B(x) andĀ(0),B(0) we find where c is the normalization factor. Here is given by the product of five scalar propagators. At the same time, the correlation function A(x)Ā(0) ∼ O 4 (x)Ō 4 (0) vanishes since, due to different ordering of scalar fields in the operators O 4 andŌ 4 , the same product of scalar propagators is accompanied by a nonplanar color factor. The correlation function B(x)B(0) receives a logarithmically enhanced correction coming from the transition O 3 → O 4 (see the first two diagrams in Fig. 9) and from similar transitionŌ 3 →Ō 4 . It is easy to verify that the relations (7.23) are in agreement with the evolution equation (7.15). They coincide with analogous expressions for two-point correlation functions in a logarithmic conformal field theory [35].
It is straightforward to extend the above analysis to higher orders in the coupling. We can use the transitions (7.16) to generate higher order Feynman diagrams. Starting from order O(ξ 6 ) a new transition appears. To see this we notice that the rightmost diagram in Fig. 9 has an intermediate state of three scalar fields φ 1 . In a close analogy with tr(φ 3 1 ) operator, we can dress this state by an arbitrary number of wheel graphs. Such graphs provide higher order contribution to the transition O 3 → O 2 and modify the eigenvalues of the mixing matrix.

Strong coupling expansion
In this section, we study the properties of scaling dimensions of the operators in bi-scalar χFT 4 theory (1.1) at strong coupling ξ 1. The numerical results shown in Fig. 8 suggest that the scaling dimensions exhibit remarkable regularity at strong coupling and fall into two different groups. The first group consists of functions ∆(ξ) that start at zero coupling at ∆(0) = 3, 7, 11, . . . and behave at strong coupling as ∆(ξ) = 2 + id(ξ) , with d(ξ) ∼ ξ 3/2 .
The second group consists of functions that start at ∆(0) = 5, 9, 13, . . . , take complex values for ξ = 0 and scale at strong coupling as ∆(ξ) ∼ ξ. To explain these properties, we solve the Baxter equation (4.20) with quantization conditions (5.23) at strong coupling using semiclassical methods. We remind that in planar N = 4 SYM theory the scaling dimensions of operators are identified through the AdS/CFT correspondence with energies of classical strings on the AdS 5 × S 5 background. The analysis in this section suggests that asymptotic behavior of the scaling dimensions in strongly coupled bi-scalar χFT 4 theory is described by a classical integrable conformal spin chain with a finite number of non-compact spins. We postpone its detailed exploration to a future publication.

Baxter equation at strong coupling
Our strategy in this section is to solve the Baxter equation (4.20) and, then, use the quantization conditions (5.23) and (5.25) to find ∆(ξ) at large ξ.
It is convenient to change variables u = iv and introduce notation for By definition, these functions have third-order poles at negative integer v and for large positive v they behave as In addition, they satisfy the Baxter equation where the notation is introduced for Here we used (5.25) to replace m = iξ 3 . As we show below, the scaling dimensions are real functions of m and, therefore, the second solution m = −iξ 3 leads to a complex conjugated expression for ∆. To find the function ∆(ξ) for ξ 1 we have to solve the quantization condition (5.23) that takes the form

Shortcut to the solution
In this subsection, we present a shortcut to finding the exact solution to (8.6). For this purpose, we concentrate on the states that start at zero coupling at ∆(0) = 3, 7, 11, . . . (the remaining states will be discussed in Sect. 8.6). As mentioned above, at strong coupling their scaling dimensions scale as (8.1) with a real valued function d(ξ) depending on the state. We introduce integer positive N to count different functions d(ξ) in the order in which the corresponding functions ∆(ξ) approach the plane Re ∆(ξ) = 2 (see Fig. 8), i.e. N = 1 for ∆(0) = 3 state, N = 2, 3 for ∆(0) = 7 states, N = 4, 5 for ∆(0) = 11 states and so on.
It is convenient to invert the dependence d = d(ξ) and introduce the function We expect from (8.1) that ϕ(d) should approach a constant value for d → ∞ and look for its general expression in the form where expansion runs in even powers of 1/d.
At large v the solutions to the Baxter equation (8.4) can be constructed as a formal series in 1/v. The leading term of the expansion is given by (8.3), e.g.
and similar forq(v). The coefficients q 1 , q 2 , . . . can be found by plugging the expansion into (8.4) and matching the coefficients in front of different powers of 1/v. They are given by rational functions in d and have a power-like behavior at large d.
Let us now examine the same combination of the q andq functions that enters into the quantization condition (8.6) where dots denote terms suppressed by powers of 1/d 2 . Notice that F 1 , F 2 , . . . scale at large d as F k ∼ d k+1 and the coefficients in front of powers of d are invariant under ϕ i → −ϕ i . We would like to emphasize that the relation (8.10) holds at large v whereas in order to solve the quantization condition (8.6) we need to know the function F (v) for small v. One may try to resum the series in (8.10) and analytically continue it to small v. As we show below, this proves to be a nontrivial task since the function F (v) has complicated analytical properties along the positive v−axis. Instead of following this route, we use (8.10) and impose the following additional condition: for large but fixed v and d → ∞ the function (8.10) should scale as with N = 1, 3, 5 . . . . As we will see in a moment, this condition fixes unambiguously the coefficients of the expansion (8.8) and yields a prediction for the function d(ξ) in (8.1), which is in a perfect agreement with the numerical results shown in Fig. 10. However it is not obvious a priori why the relation (8.12) is equivalent to the exact quantization condition (8.6). We clarify this issue in section 8.3. Let us examine the relation (8.12) for a few values of N .

States with odd N
For N = 1 we find from (8.11) that the condition This fixes the values of the coefficients ϕ 0 = ±1/8 and ϕ 1 = ∓1/8. It is quite nontrivial that the same procedure can be applied to higher F k since the number of terms to cancel grows with k and the system of equations for the ϕ−coefficients become overdetermined. We checked explicitly up to F 23 term that this system has two solutions for (8.8) and the second one that differs by the sign. Substituting these expressions into (8.11) we verify that the function (8.10) has the expected asymptotic behavior (8.12). Moreover, its expansion simplifies significantly and can be easily resummed (8.14) The relations (8.7) and (8.13) fix the dependence d = d(ξ) for N = 1. As can be seen from Fig. 11, the resulting function d N =1 (ξ) is in agreement with the numerical results for the state with ∆(0) = 3. For N = 3 the relation (8.12) does not lead to any condition for F 1 in (8.11) but requires the vanishing of the O(d 4 ) term in F 3 . This gives a new solution with ϕ 0 = ±3/8. As in the previous case, the condition F k = O(d 2 ) for k = 5, 7, . . . allows us to determine high ϕ−coefficients leading to  expression for the function (8.10) and (8.11) We verified that the function d N =3 (ξ) defined by (8.7) and (8.15) is in agreement with the numerical value for |∆ − 2| for one of the two states with ∆(0) = 7. For N = 5 going along the same lines as before we find from (8.12) These expressions correctly describe numerical values of |∆−2| for one of the two states with ∆(0) = 11. The same pattern persists to higher odd N and ∆(0).

States with even N
So far we demonstrated that the relation (8.12) describes half of the states with odd ∆(0) ≥ 5. The question arises of how to get the remaining states corresponding to even N = 2, 4, . . . . In an analogy with (8.10) and (8.12), we expect that they should satisfy the condition Its expansion at large v looks as where dots denote terms suppressed by powers of 1/v 2 and 1/d 2 . We find that for N = 2 and N = 4 the conditions , respectively, lead to the system of equations for the ϕ−coefficients whose solutions are These expressions have a simple form suggesting a generalization to arbitrary N (see Eq. (8.26) below).

Summary of the Results
Examining the expressions for the functions ϕ N = ξ 3 /d 2 with odd and even N found above we notice that they are described by the following universal formula x 3 x 2 x 1 x 4 x 2 x 3 x 1 x 4 x 3 x 2 where only odd powers of N appear. Inverting this series we can find the function d(ξ) = |∆(ξ) − 2|. We depicted this function for various states on Fig. 10. It is convenient however to examine the following combination that enters into the Baxter equation (8.5) We obtain that at strong coupling it is given by We expect that this relation holds for all states with ∆(0) = 3 + 4n for n = 0, 1, 2, . . . . For each given n > 0 it describes two states with N = 2n and N = 2n + 1. In the next section we will see an independent confirmation for this statement. Finally, the function F N (v) is given for general N by We would like to emphasize that this relation was derived at v 1 and it cannot hold for small v. Indeed, according to its definition, Eqs. (8.10) and (8.19) for odd and even N , respectively, the function F N (v) is built out of solutions to the Baxter equation which are analytic for positive v and have poles at negative integer v. As a consequence, F N (v) cannot have poles at positive v.

Asymptotic regime
In this and in the following subsections we justify the quantization conditions (8.12) and (8.18). Following the outline described in Sect. 8.1, we shall examine the Baxter equation (8.4) in two regimes corresponding to different range of the parameters, 1 v ξ and v = O(ξ). In what follows we refer to them as to asymptotic and quasiclassical regimes, correspondingly.
The asymptotic regime corresponds to the limit in which transfer matrix t(v) defined in (8.5) is large, Here in the second relation we introduced the notations for the roots of t(v). At large d we find that two of the roots are large, λ ± ∼ ±d/(2 √ 2), and the remaining one λ 0 = 4ϕ(d) stays finite in the large d limit. As a consequence, the condition |t(v)| 1 is verified for In this regime, we can distinguish two solutions to (8.4), q + (v) and q − (v), such that q + (v) q + (v + 1) and q − (v) q − (v + 1). In the leading large d limit they satisfy with v satisfying (8.28).
Having constructed the solutions q + (v) and q − (v) at large v, we use (8.29) to continue them to finite v. Since t(v) ∼ d 2 for v = O(d 0 ) the solution q + (v) becomes exponentially large when approaching the lower bound in (8.28), q + (v) ∼ d +αd with some α ∼ 1. In the similar manner, q − (v) becomes exponentially small q − (v) ∼ d −αd−1 in the same limit 13 . To build the solutions satisfying (8.29) explicitly we write them in terms of the well-defined and uniquely fixed q(v) andq(v) in the following way where we used notations for q 0 = q(0),q 0 =q(0) and q +,0 = q + (0). Indeed, generically each of the two independent solutions q(v) andq(v) will contain both q + and q − . As q + is defined up to an arbitrary addition of q − we use this freedom to define q + as in the first line above. Next, to extract q − from, say, q(v) we have to project out the growing part q + (v), i.e. we have to find some coefficient α such that q(v) − αq + (v) is not exponentially large at v ∼ 1. This can be done simply by requiring the difference to vanish at v = 0 14 . Inverting the above relations can now write q(v) andq(v) in terms of the solutions to (8.29) By construction, these expressions are dominated for v = O(d 0 ) by the first term. We recall that q(v) andq(v) satisfy Wronskian relation q(v + 1)q(v) −q(v + 1)q(v) = i d, (see Eq. (5.14)). Applying (8.31) and taking into account (8.29) we find that it leads to the following relation for q − (v) and q + (v) . (8.32) For v = O(d 0 ) the expression on the right-hand side scales as O(1/d), in agreement with the expected asymptotic behavior q Let us now examine the function F (v) defined in (8.10). Replacing the functions q(v) andq(v) by their expressions (8.31) we find . By construction, this relation holds for v inside the region (8.28). For v = O(d 0 ) the four terms on the right-hand side of (8.33) have different scaling behavior at large d: the first term grows exponentially, the last term decreases exponentially whereas the two terms in the middle have a power-like behavior. We notice however that the first term in (8.33) involves the same combinationq 0q0 + q 0q0 that enters the exact quantization condition (5.23). In other words, for the function d(ξ) satisfying the quantization condition (5.23) the first term in (8.33) vanishes leading to a power-like scaling of the function F (v) at large d. Turning the logic around, we find that the condition for the function F (v) do not have an exponential growth at large d is equivalent to the quantization condition (5.23). As we demonstrated in the previous subsection, the same condition fixes the function ϕ(d).
We can also make an additional consistency test of this observation. After canceling the leading term in (8.33) we are left with the two finite terms which we can evaluate explicitly. We can use (8.29) to get Here in the second relation we replaced t(v) by its expression (8.27) and took into account thatṫ(v) is given by the same expression (8.27) with the sign of all roots flipped. It is then easy to see that the contribution of large roots λ ± cancels in the ratio of the transfer matrices and it only depends on the root λ 0 = 4ϕ(d). Subsequenty applying (8.34) we obtain where the normalization constant is related to the value of the same product at the lower boundary of the region (8.28), c + = q − (4ϕ)q + (4ϕ)/(4ϕ) 3 . Repeating the same calculation forq − (v)q + (v) we geṫ The normalization constants c + and c − are not independent. Taking the product of the last two expressions and making use of (8.32) we obtain To find the dependence of c ± on d, we can examine the relation (8.35) and (8.36) for v close to the upper bound in (8.28) and match them into analogous expressions obtained in the quasiclassical regime v = O(d). One can show that this leads to c + /c − = d 16ϕ and, as a consequence, for positive ϕ = O(d 0 ) the contribution of (8.36) to F (v) is suppressed by a power of 1/d as compared to that of (8.35). We observe that for ϕ = N/8 expression on the right-hand side of (8.35) becomes a rational function of v. Its contribution to (8.33) matches precisely the expression (8.26)! This provides the additional support to the results presented in the previous subsection.

Quasiclassical regime
The Baxter equation (8.4) has an interesting scaling behavior for v = O(ξ). In this regime, it is convenient to introduce a new variable x = v/ξ and look for solution to (8.4) in the form of WKB expansion where 1/ξ plays the role of Planck constant and the quasimomentum p(x) admits an expansion in powers of 1/ξ The normalization of q(v) depends on the choice of x 0 . Plugging the expansion (8.39) into the Baxter equation (8.4) and expanding both sides in 1/ξ we obtain [36,37] Solving the first equation we find that e p0(x) is an analytic function on a complex x−plane with two cuts running the between the branch points at which e p0(xi) = ±1. Introducing the function we find that the branch points satisfy y(x i ) = 0 leading to and [x 4 , ∞) on the real axis, so that the middle cut shrinks into a point for d → ∞. At large positive x we can define two branches of the quasimomentum p 0,± (x) = p 0 (x ± i0) corresponding to its value above and below the cut, correspondently, where we replaced (−d) 1/2 = id/(2ξ) + O(1/d). Substituting this expression into (8.38) we find that the corresponding solution to the Baxter equation scales at large v and d 1 as q(v) ∼ v id/2 and q(v) ∼ v −id/2 , in agreement with (8.3).
Applying the above relations we can define the semiclassical solutions to the Baxter equation q(xξ) andq(xξ) on the complex x−plane with the two cuts. We recall that the quantization condition (8.6) involves another pair of solutionsq(v) = q(v, −ξ) andq(v) =q(v, −ξ). In the semiclassical approximation, they can be obtained from q(xξ) andq(xξ) through the transformation x → −x and ξ → −ξ. The resulting expressions forq(xξ) andq(xξ) have analytical properties similar to those of q(xξ) andq(xξ) with the only difference that their cuts run along intervals (−∞, −x 4 ], [−x 3 , −x 2 ] and [−x 1 , ∞) on the real axis. Since x 1 = −x 4 = (d/4) 1/2 in the large d limit, the two sets of functions share the same semi-infinite cuts whereas the two 'short cuts' are symmetric with respect to the origin.
Let us consider the following ratio The exact quantization condition (8.6) requires this function to vanish for v → 0. In the previous subsection we demonstrated that this condition is equivalent to the requirement for F (v) to have a power-like behavior (8.12) and (8.18) for v satisfying (8.28). As we will see in a moment, in the quasiclassical regime, for v = O(ξ), the same quantization condition follows from the requirement for (8.44) to be a single-valued function of v = xξ on a complex x−plane with the cuts to be specified below. The advantage of considering the ratio (8.44) is that it is independent on the choice of normalization of the solutions of the Baxter equation, or equivalently, on the choice of the point x 0 in (8.38) provided that x 0 is located away on the cuts. In the quasiclassical regime, for v = xξ, the expression on the right-hand side of (8.44) has two short cuts [−x 3 , −x 2 ] and [x 2 , x 3 ]. For the ratio (8.44) to be a single-valued function on the x−plane, it should acquire the same value after going around each of these cuts. Sinceq(v)/q(v) is analytical on [x 2 , x 3 ], the monodromy only comes from q(v)/q(v). Using (8.38) we find that the above condition translates to the Bohr-Sommerfeld quantization condition We can use the semiclassical analysis to clarify two issues that were mentioned in the previous subsections. We remind that in order to reproduce half of the states in section 8.2.2, we had to flip the sign between the two terms in the exact quantization condition (8.6) and impose the condition (8.19) instead. This can be understood as follows. The relation (8.19) holds at large v and d and it should be applicable in the quasiclassical regime for x = O(v/d) away from the cuts, that is for x 3 < x < x 4 . Since the exact quantization condition holds at v = 0 one may wonder whether the first term on the right-hand side of (8.44) acquires a monodromy as v moves from v = 0 to large v across the cut [x 2 , x 3 ]. Indeed, q(v) andq(v) are given in this case by the same expression (8.38) in which the integration goes slightly above and below the cut, respectively, so that the ratio q(v)/q(v) generates the additional factor which flips the sign for even N .

WKB expansion
In this subsection we solve the Bohr-Sommerfeld quantization (8.46) and show that the resulting expression for d(ξ) coincides with (8.25). Examining (8.25) we find that the function d(ξ) has an interesting scaling behavior for large ξ and N with their ratio N = N/(2ξ) fixed Notice that the expansion runs in powers of 1/ξ 2 with the coefficients being nontrivial functions of N . We shall determine these functions exactly using (8.46).
Integrating by parts, we can rewrite (8.46) in the following equivalent form 49) where N = N/(2ξ) and the integration goes along the contour encircling the cut [x 2 , x 3 ]. Replacing the quasimomentum with its expansion (8.39) in powers of 1/ξ we obtain analogous expansion for a a = a 0 (q) + a 1 (q) ξ + a 2 (q) ξ 2 + . . . . (8.50) The first term of the expansion is given by The integral on the right-hand side has a simple interpretation in terms of a Riemann surface defined by a complex (elliptic) curve (8.41). Namely, a 0 is given by the period of the 'action' differential over the α−cycle. For d → ∞, the α−cycle shrinks into a point and the integrand develops a pole at x = 1/d. Expanding the integrand around this point we can easily compute the integral by residues where the second relation was obtained from summing the series in 1/d and it holds for d 3 > 27.
For the second term in (8.50) we find where the integral takes a universal form and its value does not depend on d. We observe that the contribution of a 1 to the left-hand side of (8.49) is given by 1/(2ξ) and it matches an analogous term on the right-hand side. For the third term in (8.50), a 2 = − 1 2πi α dx x p 2 (x), we find after some algebra a compact representation in terms of hypergeometric functions As before, the second relation holds for d 3 > 27.
It is straightforward to continue this procedure and compute subleading terms in (8.50). In this way, we found that all terms with odd powers of 1/ξ vanish a 3 (d) = a 5 (d) = · · · = 0 (8.55) In the Appendix D we give expressions for the subleading terms a 4 , a 6 and a 8 .
Finally, we substitute the obtained expressions for a 0 , a 1 and a 2 into (8.49) and obtain the Bohr-Sommerfeld quantization condition (8.56) To find the dependence of d on the coupling ξ and N = N/(2ξ), we have to invert this relation. We verified that at large d this yields (8.48). We also checked that (8.56) is in a perfect agreement with the numerical results for d 3 > 27 (see Fig. 13). This gives a strong support to our conjecture that the Bohr-Sommerfeld quantization condition (8.56) gives the correct result for the scaling dimensions to all orders in 1/ξ expansion.
We would like to emphasize that the hypergeometric series on the right-hand side of (8.56) were obtained by summing up 1/d expansion. Curiously, at d = 3 these series develop a logarithmic branch cut. The reason why this singularity appears is that the two branch points of the curve (8.41) collide at this value of d. Indeed, for d = 3 the branching points are given by x 1 = −1, x 2 = 1/3 and x 3 = x 4 = 1/2. Going from d 1 to d = 3 we find that the α−cycle encircling the cut [x 2 , x 3 ] expands whereas the β−cycle encircling the interval [x 3 , x 4 ] shrinks into a point. For d < 3 the branch points x 3 and x 4 develop imaginary part, x 3 = x * 4 , and move away from the real axis. The definition of the α−cycle becomes ambiguous in this case since there are two possible choices of the cuts, [x 2 , x 3 ] and [x 2 , x 4 ]. The integral (8.49) evaluated along the closed contour encircling these cuts takes different complex conjugated values but its real part is the same for the two cuts. This suggests that for an arbitrary positive d the Bohr-Sommerfeld quantization (8.49) condition should look as For d > 3 this relation coincides with (8.49). Evaluation of the integral on the left-hand side of (8.57) leads to (8.56) with hypergeometric series replaced by their real part. We verified that this relation correctly describes numerical results for d < 3 (see Fig. 13). . We have shown that it can be found by solving the Bohr-Sommerfeld quantization condition (8.57) for a particular choice of the cycle α on the Riemann surface (8.41). This cycle is uniquely defined by the requirement that the branch points should merge for d → ∞. In this section we argue that the remaining states with ∆(0) = 5, 9, . . . also admit an analogous semiclassical description. The only difference compared to (8.57) is that these states correspond to another choice of the integration contour on the Riemann surface (8.41).
We remind that the branch points satisfy the relation (8.42). There are four different values of d for which two of these points collide. Two of them, d → ∞ and d = 3, we encountered in the previous subsection. As we will see in a moment, the two remaining ones, d = 3 e ±2πi/3 , describe the strong coupling limit of the states ∆(0) = 5, 9, 13, . . . . As follows from the definition (8.24), the value d = 3 e ±2πi/3 corresponds to the following scaling behavior of ∆(ξ) at strong coupling It should be compared with analogous relation (8.1) for the states ∆(0) = 3, 7, 11, . . . . We verified that (8.58) correctly describes numerical results for the trajectories in Fig. 8 that start at ∆(0) = 5, 9, . . . . To find subleading corrections to (8.58) we shall employ the Bohr-Sommerfeld quantization condition analogous to (8.57). Following the logic of the previous subsection, the integration contour should be chosen to encircle the pair of branch points that collide at d = 3 e ±2πi/3 . It is easy to see from (8.41) that for this value of d the branch points are aligned along the same ray in a complex plane, x 1 = −e ∓2iπ/3 , x 2 = e ∓2iπ/3 /3 and x 3 = x 4 = e ∓2iπ/3 /2. Thus, in the vicinity of d = 3 e ∓2πi/3 the integration contour should encircle the segment [x 3 , x 4 ]. This corresponds to the choice of the β−cycle on the Riemann surface (8.41) (see Fig. 12). The resulting Bohr-Sommerfeld quantization condition reads where N = N/(2ξ) and the quasimomentum is given by (8.39) and (8.40). The relation (8.59) has been previously encountered in the study of semiclassical limit of the SL(2) spin chain [37]. Using the results of [37] we find that in the vicinity of d = 3 e ±2πi/3 the first two terms of the expansion of a D in powers of 1/ξ are given by where d = 3 e ±2πi/3 (1 + x) and x is small. Here O(1/ξ) term describes the contribution of p 1 (x) correction to the quasimomentum. As in the previous case, it cancels against the same term on the right-hand side of (8.59). The series in the first term in (8.60) can be summed up and expressed as discontinuity of a hypergeometric 3 F 2 −series More precisely, the hypergeometric series develops a logarithmic cut that starts at d 3 = 27. The discontinuity across this cut is a rational function of d = 3 e ±2πi/3 (1 + x) whose expansion at small x matches the first term on the right-hand side of (8.60). Notice that the same hypergeometric series enters (8.52). This is not accidental of course since the relations (8.61) and (8.52) define the period of the same 'action' differential over the two cycles on the Riemann surface (8.41).
Solving the Bohr-Sommerfeld quantization condition (8.59) we can obtain the large ξ expansion of d To find the leading term we substitute (8.60) into (8.59) and invert the series where the notation was introduced for complexN = N e ±iπ/6 = N e ±iπ/6 /(2ξ). In agreement with our expectations, this relation defines two complex-valued functions ∆(ξ) (see (8.24)). We verify that at large ξ these functions have a correct asymptotic behavior (8.58).
To determine the subleading correction to (8.62), we exploit the relation between the periods a and a D mentioned above. It allows us to find O(1/ξ 2 ) correction to a D by taking a discontinuity of (8.54). In this way, we find  where expansion runs in powers of e ±iπ/6 ξ. The relations (8.65) and (8.66) are valid at large ξ with N fixed. We checked them against numerical data on Fig. 14 for different N . We found that the identification of integers N corresponding to different states is not trivial. Namely, for odd N the relation (8.65) correctly describes the states with ∆(0) = 5, 9, 13, . . . whereas for even N it predicts some states which are not present in the spectrum. To some extend that is expected as these states come in the complex conjugate pairs and in order to have the same average number of states in some interval of the quantum numbers N as for the states with ∆(0) = 3, 7, 11, . . . we should miss exactly half of all N 's.

Problems and perspectives
In this paper, we used the quantum integrability methods to compute the exact scaling dimensions of various operators in four-dimensional bi-scalar theory defined by the Lagrangian (1.1) (also called χFT 4 ) in planar limit. This theory has the symmetry SU (2, 2) × U (1) 2 -a subgroup of the original P SU (2, 2|4) conformal symmetry of the full (untwisted) N = 4 SYM theory from which it was obtained in a special double scaling limit combining large imaginary γ-twist and small coupling. Consequently, the operators in this theory are classified by the set of Cartan charges (∆, S,Ṡ|J 1 , J 2 ). The basic example of such operators is tr(φ J 1 ) -the "BMN vacuum" operator 15 with charges (∆, 0, 0|J, 0). The perturbative corrections for this operator are given by so called "wheel" graphs of Fig.1 (a single graph at each non-zero order). Our method makes it possible to compute these conformal graphs at least to O(1/ ) order in dimensional regularization, in an algorithmic way. We computed the J = 3 wheel graphs up to 12 loops, which can be easily pushed further with more computer time or on a more powerful computer. Our numerical computation of anomalous dimension have been done for this operator at virtually unlimited precision for all reasonable couplings ξ. Moreover, it is also applicable to the other operators, of the type (1.2), with the same charges but a bigger length L corresponding to the insertion of any number of couples of fields φ 1 , φ † 1 and φ 2 , φ † 2 inside tr(φ J 1 ), as well as of the derivatives, where we performed similar analytic and numerical computations. Furthermore, we investigated the strong coupling limit of the anomalous dimensions for this family of operators. In this limit we found that the spectrum is described by a classical algebraic curve and the quantization condition reduces to a modified Bohr-Sommerfeld quantization condition with half-integer filling fractions. We managed to develop a systematic strong coupling expansion to a high order in 1/ξ.
To compute the spectrum of scaling dimensions ∆, we needed to provide two basic ingredients: i) formulate a 4th order finite difference Baxter equation for a quartet of Q−functions of spectral parameter; ii) derive quantization condition fixing the unique solution with prescribed analytic structure, and at the same time all the constants in Baxter equation and the spectrum of dimensions ∆(ξ) with given charges.
We managed to solve the problem i) -to find the Baxter equation at any J for the operators in the sector (∆, 0, 0|J, 0) . We did it in two ways: from AdS 5 /CFT 4 QSC formalism (so far only for the charges J = 2, 3, 4), where the full Baxter equation is known, in the double scaling limit leading to the bi-scalar theory; and also from explicit Lax operator of SU (2, 2) conformal spin chain describing the underlying fishnet graphs of the bi-scalar model [1]. We used the fact that the Baxter equation is a universal object in quantum integrability, independent of the auxiliary representation of Lax operator. The resulting Baxter relation, even though cross-checked in many ways, is not derived in a rigorous way in none of two formalisms. Certain very natural assumptions of symmetry of fundamental transfermatrices has been used, but not yet formally proven, in derivation from Lax operators. It would be good to check these assumptions algebraically, directly from the from the spin chain formalism. The QSC derivation serves as a good cross-check but it also involves certain natural, but unproved assumptions about the scaling behavior of QSC Q-functions in the double scaling limit. One of the ways to verify these assumptions is to use numerics for the full twisted N = 4 SYM and then approach the limit with infinite coupling numerically.
Concerning the ingredient ii), we managed to derive the quantization condition from the QSC formalism, so far only in the case J = 3, where the 4th order Baxter equation nicely factorizes into two 2nd order equations. Its derivation for an arbitrary J can be certainly obtained in a similar way but we leave this, likely more involved, calculation to future publications. In particular for J = 4 one should reproduce [1] γ J=4 vac = −40ζ 5 ξ 8 + 8 309ζ 11 +16ζ 3,8 +20ζ 5,6 −4ζ 6,5 +40ζ 8,3 −8ζ 3,3,5 +40(ζ 3,5,3 +ζ 5,3,3 )−200 ζ 2 5 ξ 16 +O(ξ 24 ) . (9.1) It would be also interesting to derive the quantization condition in a rigorous way, from Lax formalism of the underlying conformal spin chain. It would involve the computation of certain spin chain eigenfunctions and matrix elements relevant to the "wheel" graphs. Probably, the most adequate formalism for it should be the Sklyanin's separated variables (SOV). Such a formalism would not only allow to compute in a rigorous way the wheel graphs but also get hand on a much bigger variety of physical quantities: more general than in [4] fishnet-type graphs at arbitrary loop orders, such as described in [14], multi-point correlators and amplitudes [17] in χFT 4 , etc. Recently observed simplifications in the SOV approach to the wave functions [38] could pay an important role here.
Our methods are potentially applicable to the χFT 3 theory obtained in the double scaling limit of twisted ABJM [14], where certain classes of Φ 6 -type graphs can be computed at any loop order. One can try to apply for it both the Lax [31] and the QSC formalisms [15,16]. One could also try to approach the same program for tri-scalar χFT 6 with Φ 3 -type chiral interactions, which seem to define the genuine CFT in planar limit [18]. The 6D twisted SYM "mother"-theory, from which the latter model could be obtained in a double scaling limit, is unfortunately not known.
It would be interesting to include into our formalism more general χFT 4 containing fermions and more of scalars, described in [1,14]. This could be done by applying more sophisticated limits to AdS/CFT QSC equations or, alternatively, establishing the Lax formalism for conformal spin chain with fermionic degrees of freedom. The latter would open a way to multi-loop computations of more sophisticated than Φ 4 type graphs, containing also the Yukawa couplings.
It would be also interesting to compute next-to-leading corrections to the double scaling limit of γ-twisted N = 4 planar SYM which might give a better understanding of the origins of integrability in the full theory. Another, more ambitious direction, would be to generalize the results of the current paper to the full P SU (2, 2|4) Heisenberg spin chain, which could led to the prove of integrability of N = 4 planar SYM from first principles.
representation is parameterized by a set of 4 complex numbers ν := (ν 1 , ν 2 , ν 3 , ν 4 ), such that Complex z ij variables can be used to define the lower-triangular matrix where (ê ij ) l k := δ ik δ l j are the unit 4 × 4 matrices. Define further the matrix differential operator The algebra sl (4) can then be represented in Aut [V (z)] by defining its generators as whereν = diag(ν 1 , ν 2 , ν 3 , ν 4 ). The corresponding sl(4)-module V ν is irreducible [42,43] iff ν ij := ν i − ν j / ∈ N , ∀i < j, otherwise V ν contains an invariant subspace. Let us consider a set of complex numbers σ = {σ ii+1 } such that Then it can be shown [40] that the associated principal series sl (4)-module V (k) σ decomposes as a tensor product: One of the factors is the infinite-dimensional space V 3−k of polynomials in the variables appearing in the first 3 − k columns and rows of Z. The other factor v σ k+1 is a finite-dimensional sl (k + 1) module with highest weight Λ σ = (λ k , · · · , λ 3 ) and λ j = σ j j+1 − 1. Consequently, choosing k = 3, we have that the principal series module V (3) σ collapses to a finite dimensional sl (4) module: This module is associated to the Young tableau defined by the partition = { 1 , 2 , 3 } with j = 3 k=j λ k .

The general T − Q relation from sl(4) invariant R-matrix
Let us consider the following invariant R-matrix: where both V ν and V σ are sl (4) principal series modules, defined above. As was shown in [25,41], this operators adimt a factorised form: where P ij is the permutation operator of spaces i and j. The associated T -operator, defined in the usual way as the trace over the auxiliary space V σ of the product of J R-matrices 17 will then enjoy a factorised expression where the Q-operators are defined as transfer matrices associated to the factorising operators R (j) : This factorisation is a general property, irrespective of the fact that V σ is irreducible or not. If we were to choose σ such that V σ = V (k) σ , then the R-matrix will assume an upper-triangular form It is possible to prove, by using Bernstein-Gel'fand-Gel'fand resolution for finite dimensional sl (n) modules [26,42], that Finally, by choosing k = 3, we obtain the following nice determinant expression which can be rewritten in terms of Young tableaux indices as where l j = 4 − j + j , f = 1 4 and denoting the T -operators associated to the k-th fundamental representation as and understanding t 4 (v) ≡ t 0 (v), we easily obtain the general form of Baxter T -Q relation for a π ν [sl(4)]-spin chain [26] 4 k=0 The fundamental R-matrices and the associated Lax operators In order to obtain an explicit expression for the transfer matrices t k , we need to compute the traces where R (k) is the R matrix with physical space V ν and auxiliary space π Λ k and (Λ k ) i = δ i k is the k-th fundamental weight. These R-matrices for the fundamental representations are related to the Lax operators L k by a normalization: where κ k are some rational numbers and X π k ,ν (v) are functions. The Lax operator for the basic representation is a linear operator on the space V ν ⊗ 4 [sl (4)] and is defined as follows: In order to consistently define the T -and Q-operators it is necessary to introduce a regulator in the R-matrix. We choose to avoid this subtlety as it is of no relevance for our goal. Further informations can be found in [26]. 18 Note that in our conventions, the defining representation 4 is associated to the Young tableau 1 = (1, 1, 1) with three boxes, while the conjugate one 4 to the Young tableau 1 = (1, 0, 0) with a single box. This is all a matter of preference, since the Z 2 symmetry of the Dynkin diagram makes the two representations isomorphic. However the isomorphism is realised in a complicated fashion at the level of R-matrices which results in the fact that the associated transfer matrices are, in general, radically different, as we see later. Our choice was dictated by the request that t 1 should be the simplest among the transfer matrices.
The higher Lax operators L 2,ν and L 3,ν are linear operators on, respectively, V ν ⊗ 6 [sl (4)] and V ν ⊗ 4 [sl (4)]. Although they have not a simple form such as L 1,ν , from classical representation theory it is known that (4) which means that L 2,ν and L 3,ν can be defined as follows With the notation ∧ k A k (v), where A is a linear operator on a space V , we denote a linear operator on the space k V , acting as follows As it will be useful in a short while we also introduce the quantum determinant of the Lax operator L 1,ν : (A. 31) In order to fix the normalisations κ k and X k,ν in (A.26), we apply the R-matrix on the highestweight state Φ 0 = 1 ∈ V ν ⊗ V σ , using the known expression for the eigenvalues of the general R-matrix with p (v) being a periodic function whose exact expression is irrelevant for our needs. Plugging in the values of the labels σ corresponding to the fundamental representations (see table 1), we obtain where we denoted as r k (v) the eigenvalue of R (k) (v) on the highest weight state Φ 0 and we introduced the notation Now we will compare these eigenvalues with those obtained by acting with the Lax operators L k,ν on the highest weight states (A. 35) These are easily computed by using the explicit representation of the operators E ji (A.7) and the definition (A.30): We also easily obtain By direct comparison we then make the following identifications Note that we might have identified Lax operators with fundamental R-matrices also as giving us the relations which can be rewritten as These relations are the generalisations to the quantum case of the equalities valid for any 4 × 4 matrix [44].

The monodromy matrix as a Manin matrix and Talalaev's formula
From now on we will drop the explicit dependence on the labels ν.
Thanks to the identifications (A.38) and defining the Monodromy matrix M (v) of the spin chain as follows we are able to write the general Baxter equation (A.24) in the following form Now we notice that the following relations hold These are the defining relations for a Manin matrix [45], a particular class of matrices with not necessarily commuting entries. They are a natural extension of the matrices over commutative rings, in the sense that most of the standard theorem of linear algebra continue to hold true for them [46]. Of particular interest for us is the fact that the spectral determinant of M (v) e ∂v has the following expansion where "cdet" stands for the column ordered determinant As a consequence, the Baxter T − Q equation (A.44) takes the following suggestive expression which is the quantum version of a classical spectral curve equation. The above expression has first appeared in the works of Talalaev [47] and we thus refer to it as "Talalaev's formula".
The explicit computation of the transfer matrices More important on the practical level is the fact that the higher traces satisfy the Newton's identities [48] kt k (v) e k∂v = k j=1 Explicitly we have Thanks to these relations we are now in the position to derive the expressions for the coefficients of Baxter T − Q equation (A. 44). To this end, we introduce the following objects One easily see that and, since 19 q 1 = 0, With some basic algebra we obtain the following expressions is global, in the sense that it only depends on the labels of the representation V θ of the full spin chain 20 : it is thus a Casimir operator and we call it "total Casimir". On the other hand, not all of the charges q [i1,...,im] are independent. There exist relations amongst them which also involve the usual Casimir operator The following functions possess an interesting symmetry property. In fact let us explicitly denote their dependence on the local conserved charges: , (A.65) 20 The Hilbert space of the spin chain decomposes as a direct integral of representations V θ [49]: The eigenvalues of the charges C (J) m only depend on which of these summands the state of the spin chain belongs to, not on the specific state considered.
Up to now, we have used solely the representation theory which allowed us to fix 9 coefficients amongst the 6J ones. We need further constraints. As it has been shown in Sect. 4, considering the equation (A.77) as a double scaling limit of the Baxter equation in γ-twisted N = 4 SYM (see Eq. (4.13)) means we must be able to bring it to the symmetric form (3.1). Imposing these symmetries on the functions (A.78), will allow us to fix some of the unknown coefficients of the functions τ k . First of all, we see that, in order to have the coefficients in front of q(v + 1) and q(v − 1) in (A.77) to be, respectively, B(v + 1 2 ) and B(v − 1 2 ) for some function B(u), the transfer matrices have to satisfy the following relation From this requirement it directly follows that the transfer matrices have then a definite parity in this case, τ 1 (v) = (−1) J τ 1 (−v), τ 2 (v) = τ 2 (−v) and τ 3 (v) = (−1) J τ 3 (−v), and are given by Finally, asking the symmetry (A.79), the Baxter equation takes the form (3.1) and (3.4) advocated in section 3. For J = 2, J = 3 and J = 4 the expressions (A.84) contain, respectively, 0, 1 and 3 unfixed coefficients.
As it has been showed in Sect. 5, the rest of coefficients have to satisfy the additional quantisation conditions. The relation (A.82) imposes a tight selection rule on the states of the spin chain. It is these states that play a distinguished role in our analysis as they allow us to find the scaling dimensions of the operators (1.2) for an arbitrary coupling.
Let us examine relations (A.78) and (A.84) for few lowest values of the length J of the spin chain.

B QSC supplementary relations
Coefficients of 4th order Baxter equation The determinants D i used the equation (4.13) are defined as follows:

Formulas for ∆ through ansatz coefficients
This appendix contains expressions for ∆ through the coefficients of the anatz (C.4).
In order to proceed we need to constrain the coefficients c a,k of the ansatz above as much as we can. We also want to be able to express ∆ through c a,n . To this end, we plug P a given by the ansatz into (4.13) and expand this equation at u → ∞ assuming that the asymptotics of Q i are given by (4.6). This yields a system of equations for coefficients c a,n and ∆. The structure of this system of equations depends on L. Typically we have to expand the Baxter equation to the fourth subleading order in 1/u to get a non-trivial relation for ∆ and a closed system for c a,n . We performed the calculations for J = 2, 3, 4.
For J = 2 and J = 3 the resulting relations for ∆ have a form where f 2 , f 3 are rational functions of their arguments given in appendix B. Now we need to take the double scaling limit of the relations just derived and of the equation (4.13). Let us define (see (4.8)) s = √ κκ , r = (κ/κ) J , ξ = gs . (C.8) The double scaling limit consists of taking g → 0 and s → ∞ when keeping ξ constant. The parameter r is not entering the Lagrangian (1.