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 N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 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 N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 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 tr(ϕ1J) where ϕ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 N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 SYM.


JHEP01(2018)095
Contents 1 Introduction: bi-scalar theory and results for "wheel" graphs 1 2 SU(2, 2) picture for fishnet graphs 7 2.1 "Graph-building" transfer matrix and its integrable SO (1,5)  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 JHEP01(2018)095 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 multiloop 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 biscalar 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 4-point 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 JHEP01(2018)095 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 figure 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 JHEP01(2018)095 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 JHEP01(2018)095 finite-difference equations of the form This equation depends on two integrals of motion, ∆ and m, and coincides with the Baxter equation for the SL(2) spin chain of length J = 3 and spin 0 representation at each site. The quantization condition which fixes the dependence of ∆ and m on the coupling constant is 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 figure 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 figure 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 figure 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 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 JHEP01(2018)095 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 figure 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 figure 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.

SU(2, 2) picture for fishnet graphs
To prepare the formalism for computing the wheel-graphs shown in figure 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 figure 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 figure 4. The contribution of such graphs to (2.1) can be written as where T J,M (x 1 , . . . , x J |y 1 , . . . , y J ) describes a cylinder-like fishnet graph with M wheels and J external points. For M = 0 we have so that the corresponding contribution to (2.2) is given by the product of J scalar propagators. For M ≥ 1, due to the structure of wheel graphs, T J,M admits the following recursive representation

JHEP01(2018)095
where the kernel H J represents one wheel of the diagram in figure 1. It is given by the product of scalar propagators with periodic boundary condition y 1 = y J+1 . Here each horizontal link is a free scalar propagator (2π) −2 (y j − y j+1 ) −2 . Each vertical line produces similar propagator connecting points x i and y i . Replacing the scalar propagator with inverse d'Alembert operator, we find that (2.5) can be identified as a kernel of the following "graph building" operator defined in [1] 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 JHEP01(2018)095 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).
2.1 "Graph-building" transfer matrix and its integrable SO (1,5) 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 x = (x 1 , . . . , x J ) 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 5 Except the case J = 2 which has overlapping singularities and leads to the leading asymptotics ln D(x) ∼ −1/ 2 [1,3]. 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.

JHEP01(2018)095
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] R u (x 1 , x 2 |y 1 , y 2 ) = c(u) 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 figure 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

JHEP01(2018)095
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 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 figure 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 (2.20) It acts on the tensor product of J copies of the conformal group representation (2.10) and the trace is taken over the same (infinite-dimensional) auxiliary space representation.
Replacing each R−operator in (2.20) by its integral representation (2.13) we can realize T J (u) as an integral operator where Φ is a test function and the kernel T J,u (x|y) is given by a J−fold integral over the auxiliary space

(2.22)
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 figure 6. In the standard manner, in virtue of Yang-Baxter equation, the transfer matrices (2.20) commute among themselves as well as with the local integrals of motion of the spin chain. Let us examine (2.21) and (2.22) for two special values of the spectral parameter, u → 0 and u → −1. In both cases, some of the propagators in the integral representation of the transfer matrix have index 2 which allows us to apply the identity (2.18). For u = − and → 0 we find Figure 6. Diagrammatic representation of the transfer matrix (2.22). It is obtained by gluing together J rectangles, each representing the R−operator. The leftmost and rightmost vertices are identified. The values of indices α ± and β are the same as in figure 5. 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 JHEP01(2018)095 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 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 (3.8) The values of the parameters m and α = (∆ − 2) 2 as functions of the coupling ξ will be fixed in section 5 from an additional quantization condition.

Baxter equation for
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.

JHEP01(2018)095
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. Figure 7. Left: analytic structure of P a (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 Q i (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 toQ(u), which has an infinite ladder of cuts in the upper half-plane (right).

JHEP01(2018)095
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 figure 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 figure 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 [−2g, 2g] [32] 6 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.
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 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 Qi 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 Qi as νi for different i could differ by an integer, allowing for Qi 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 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.

JHEP01(2018)095
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 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 statê 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] .

JHEP01(2018)095
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 "Fourthorder 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.1.

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)

JHEP01(2018)095
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

JHEP01(2018)095
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.

JHEP01(2018)095
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.2) 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.3) 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].

JHEP01(2018)095
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.3) 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.3, 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.3) and (4.4) into (5.6) we find after some algebra 10 where s = √ κκ and the notation was introduced for We conclude from (5.7) 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.

JHEP01(2018)095
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.4), 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.10) that the residue at this pole have to satisfy Using the Baxter equation (5.2) we get for u → 0 Replacing u → −u + i in (5.10) and matching it into the last relation we obtain The relations (5.11) and (5.13) 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.10) and (5.3) to establish linear relations between the functions Q i (−u) to Q j (u). Using the definition (5.4) we obtain from (5.16) the leading JHEP01(2018)095 asymptotic behavior of Ω i j (u) at the origin (5.17) In the next section we compare this relation with the expression for the discontinuity of Ω i j given by (5.7) 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.17) 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.18) for g → 0 According to (5.7), ∆Ω 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.18) and taking into account (5.19) JHEP01 (2018)095 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.22) should take the expected form (5.17). In particular, it should scale as 1/u 3 for u → 0. Note that l.h.s. of (5.22) scales as g 0 . We have to deduce the value of n. If we assume that n = 4 the second term in (5.22) 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 where in the second relation we took into account (5.20) and (5.21). This equation provides us with a set of nontrivial relations between the entries of matrices (5.7) and (5.17). In particular, we notice that the matrix element ∆Ω 2 1 is zero implying that Ω 2 1 should vanish.
This leads to the quantization condition 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.24) this gives (5.25) Finally imposing the relation (5.23) 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.23) are satisfied automatically! The quantization condition (5.24) fixes the dependence of ∆ on m. Together with (5.26) 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. 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 (∆, S1, S2|J1, J2, J3). 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 (∆, S1, S2|J1, J2). For our particular BMN-state and its analogs considered further we take (∆, 0, 0|J, 0).

JHEP01(2018)095 6 Numerical solution
In this section, we describe the numerical solution to Baxter equation (4.20) supplemented with the quantization conditions (5.24) and (5.26). 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.2) and find the coefficients of the expansion by 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 figure 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.
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.24) 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.26). 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.15). It is a meromorphic function of u with the second-order poles located in the lower half-plane q I = P (L−1)/2 (u) , q II = P (L−1)/2 (u)η 2 (u) + Q (L−3)/2 (u) , (7.1)

JHEP01(2018)095
where P n and Q n are polynomials of degree n and the notation was introduced for the special function η s 1 ,...,s k (u) with an appropriate pole structure 12 [34] η s 1 ,...,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.2) at infinity. As follows from (5.2), 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)

JHEP01(2018)095
at large u and matching the coefficients we find that solutions to the Baxter equation satisfying (5.2) are given by To analyze the quantization condition (5.24), we have to evaluate these expressions at the origin. Making use of the identity η s 1 ,s 2 ,...,s k (i) = (−i) s 1 +s 2 +···+s k ζ s 1 ,s 2 ,...,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.24) 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.24), 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:

JHEP01(2018)095
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 State with L = 9. We found two states and their properties are similar to those of L = 5 states 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 bi-scalar χ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 V ij describe the mixing Figure 9. Mixing of O 3 = tr(φ † 2 φ 1 φ 2 φ 2 1 ) with the operators O 4 and O 2 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.
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 figure 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 figure 9 γ 1 = −2 , γ 2 = 1 . 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.

JHEP01(2018)095
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 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 figure 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].

JHEP01(2018)095
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 figure 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 figure 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 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.24) 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.24) and (5.26) 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

JHEP01(2018)095
where the notation is introduced for Here we used (5.26) 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.24) 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 section 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 figure 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

JHEP01(2018)095
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 figure 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 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 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 The relations (8.7) and (8.13) fix the dependence d = d(ξ) for N = 1. As can be seen from figure 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  get a simple 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

JHEP01(2018)095
For N even, large d asymptotics of this function should contain odd powers of d only and, therefore, F − (v) should change the sign under d → −d. Notice that the same transformation exchanges the two solutions to the Baxter equation q(v, ξ) andq(v, ξ). The function (8.10) is obviously invariant under q(v, ξ) ↔q(v, ξ) and, as a consequence, its coefficients (8.11) are even functions of d. To get a parity odd function, it is sufficient to flip the sign in (8.10) 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 Evaluating the function (8.19) for the solutions (8.21) we obtain

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 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 figure 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 section 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.

JHEP01(2018)095
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 13 we use the term "exponentially" in the sense faster-than-any-power grows/decay. 14 Alternatively one can set it to zero at any other v ∼ 1, this will change the definition of q− by adding to it q+ with an exponentially small coefficient.

JHEP01(2018)095
eq. (5.15)). 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 combinatioṅ q 0q0 + q 0q0 that enters the exact quantization condition (5.24). In other words, for the function d(ξ) satisfying the quantization condition (5.24) 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.24). 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 getq

JHEP01(2018)095
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 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 p 0 (x) is an analytic function on a complex x−plane with two cuts running the between the branch points at which e p 0 (x i ) = ±1. Introducing the function we find that the branch points satisfy y(x i ) = 0 leading to

JHEP01(2018)095
The position of the branch points depends on the value of d defined in (8.24). As follows from the last relation, two of the branch points coincide for d → ∞ and d 3 = 27. At large d, the branch points are located at x 1 = −(d/4) 1/2 , x 2 = 1/d, x 3 = 1/d + O(1/d 4 ) and x 4 = (d/4) 1/2 . It is convenient to choose the cuts of p 0 (x) to run along intervals (−∞, x 1 ], [x 2 , x 3 ] 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 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, −ξ) anḋ q(v) =q(v, −ξ). In the semiclassical approximation, they can be obtained from q(xξ) and q(xξ) through the transformation x → −x and ξ → −ξ. The resulting expressions foṙ q(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 exp 2ξ Here the integration contour encircles the segment [x 2 , x 3 ] and N is an arbitrary integer. For the cut [−x 3 , −x 2 ] we can get analogous relation by replacing y → −y and ξ → −ξ.
We show in the next subsection that, upon replacing the quasimomentum in (8.46) by its explicit expression (8.39) and (8.40), the relation (8.46) allows us to obtain the dependence of d on the coupling ξ.
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 righthand 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 exp ξ 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
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) 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 figure 12). 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)

JHEP01(2018)095
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 figure 12). . 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 figure 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 figure 13). 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 x 3 x 2 x 1 x 4 x 2 x 3 x 1 x 4 x 3 x 2 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 tions ∆(ξ) (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 figure 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.

JHEP01(2018)095 9 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 figure 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 transfer-matrices has been used, but not yet formally proven, in derivation from Lax operators. It would be good to check these assumptions algebraically, directly 15 It was protected BMN vacuum operator in untiwisted N = 4 SYM theory but it gets corrections after γ-twisted N = 4 SYM and, consequently, in the bi-scalar χFT4.

JHEP01(2018)095
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 SU(2, 2|4) Heisenberg spin chain, which could led to the prove of integrability of N = 4 planar SYM from first principles.

JHEP01(2018)095
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 finitedimensional 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 σ 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 .

JHEP01(2018)095
A.2 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

JHEP01(2018)095
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: The higher Lax operators L 2,ν and L 3,ν are linear operators on, respectively, V ν ⊗ 6 [sl (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 with φ k ∈ V . As it will be useful in a short while we also introduce the quantum determinant of the Lax operator L 1,ν : In order to fix the normalisations κ k and X k,ν in (A.26), we apply the R-matrix on the highest-weight state Φ 0 = 1 ∈ V ν ⊗ V σ , using the known expression for the eigenvalues of the general R-matrix [25] 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

JHEP01(2018)095
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 (4)] .
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].

JHEP01(2018)095
A. 4 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".

A.5 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 = 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 and One easily see that and, since 19 q 1 = 0, With some basic algebra we obtain the following expressions 19 This is because q1 = J m only depend on which of these summands the state of the spin chain belongs to, not on the specific state considered.

JHEP01(2018)095
Then, substituting (A.79) into the Baxter equation (A.77), we find that this equation remains invariant under u → −u and q [α] → (−1) |α| q [α] . This means that for any solution q(u, ∆, q [α] ) to (A.77) there should exist another solution q(−u, ∆, (−1) |α| q [α] ) describing the state with the same scaling dimension ∆. Thus, the states with nonzero q [α] s.t. |α| ∈ 2Z + 1 have a two-fold degeneracy with respect to ∆. To avoid the degeneracy, we have to require that all charges with odd indices should vanish [36] q [α] = 0 , ∀α s.t. |α| ∈ 2Z + 1 . (A.82) 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 As it has been showed in section 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.

A.7.1 Baxter equation for J = 2
In this case, the general expressions for the transfer matrices (A.78) are 3 v + υ

B.2 Formulas for ∆ through ansatz coefficients
This appendix contains expressions for ∆ through the coefficients of the anatz (C.4).
We assume that the state is parity invariant which implies that at the level of the P a functions we have P a (−u) = λ a b P b (u) , P a (−u) = λ b a P b (u) (B.3) for some constant coefficients λ a b and λ b a , which obey λ a b λ a c = δ b c . Let us show that Indeed Ω i j (u)Q j (u) = −Q a|i (−u + i/2)λ a b P b (u) = −Q a|i (−u + i/2)P a (−u) = Q i (−u) , (B.5) using identities like P a P a = 0 and Q i = −Q a|i (u ± i/2)P a and Q i = +Q a|i (u ± i/2)P a it is easy to check that Ω i j (u + i) = Ω i j (u).
Furthermore, we can easily find the discontinuity of Ω i j (u):

C Details of the derivation of the Baxter equation from QSC
Here we give details of the derivation of the Baxter equation from QSC approach for J = 2, 3, 4.

C.1 Left-right symmetry
Most problems solved using QSC method before possessed so-called left-right symmetry, which in particular means that the indices can be raiser or lowered with a constant matrix:

JHEP01(2018)095
States of twisted N = 4 SYM are not left-right symmetric for general twists, however for the particular case of operator trφ J 1 and twists (4.7) indices can still be raised by χ if one also exchanges κ withκ: C.2 Baxter equation from double scaling limit of QSC In this section we will derive a finite-difference equation for Q i (u) in the double scaling limit. Our starting point will be the Baxter equation (4.13). The computation is similar to that in [23]: we need to construct an ansatz for P a , expand it in the double scaling limit and plug into (4.13). The ansatz will necessarily contain many unknown coefficients which we will fix by solving (4.13) expanded at large u. In the process we will also get a relation between ∆ and the coefficients of the ansatz for P a . Let us start with an ansatz for P a : along the lines of [13], we pull out of P a the exponential and power-like prefactors, leaving the part which scales like 1 at infinity: Here we choose A a = 1 and A a according to (4.9). Remember that P a and P a have only one cut -Zhukovsky cut on the real axis. This cut can be resolved by considering P a as a function of Zhukovsky variable x(u). In other words, we can represent f i , g i as series in x(u): This ansatz was constructed so as to satisfy the condition that P a stays finite as g → 0 andP a grows as O(g −J ). One can see that this is the case if we assume c a,k ∼ 1. Indeed, the operation tilde (monodromy around a branch point of the Zhukovsky cut) transforms x(u) to 1/x(u). Since we want to keep u = g(x + 1/x) finite in the weak coupling regime, when x ∼ 1/g → ∞, each power of x should be compensated by at least one power of g. where Poly(u) is some polynomial in u. The scaling condition for P a -the finiteness in the weak coupling limit g → 0 -is obviously satisfied. Now considerP a : The scaling condition forP a is satisifed as well. To summarize: the scaling of the terms in f 1 , g 2 is constrained by the scaling ofP. The scaling of the polynomial terms in g 1 , f 2 is constrained by finiteness of P as g → 0 and the scaling of the singular terms in g 1 , f 2 -byP. We also want to have a consistency with the weak coupling result [13], that is why the potentially singular parts in f 1 , g 2 ,g 2 ,f 1 have a g 2L prefactor for the infinite sums.
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. where f 2 , f 3 are rational functions of their arguments given in appendix B.2. 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.1), so we expect it to drop out from final result, although it may be present in the intermediate computations.
Each coefficient c a,k is a regular function of g, κ,κ and so it needs to be expanded in the powers of g: c a,k = ∞ m=0 c a,k,m (ξ)(g/ξ) 2m (C.9)