Exact Correlation Functions in Conformal Fishnet Theory

We compute exactly various 4-point correlation functions of shortest scalar operators in bi-scalar planar four-dimensional"fishnet"CFT. We apply the OPE to extract from these functions the exact expressions for the scaling dimensions and the structure constants of all exchanged operators with an arbitrary Lorentz spin. In particular, we determine the conformal data of the simplest unprotected two-magnon operator analogous to the Konishi operator, as well as of the one-magnon operator. We show that at weak coupling 4-point correlation functions can be systematically expanded in terms of harmonic polylogarithm functions and verify our results by explicit calculation of Feynman graphs at a few orders in the coupling. At strong coupling we obtain that the correlation functions exhibit the scaling behaviour typical for semiclassical description hinting at the existence of the holographic dual.

. Three types of the 4-points functions topologies. These cases differ by the number of "particles" (red dashed lines) transfered from the bottom to the top. These particles can be associated in terms of the ABA with magnons in the intermediate states on the OPE. We refer to them as zero, one and two particle cases correspondingly. The black lines correspond to the X-particles.

Introduction
Conformal quantum field theories (CFT) have demonstrated their importance for very diverse fundamental problems in physics, its applications ranging from the physics of phase transitions (see [1] and citations therein) to various problems of fundamental interactions and cosmological scenarios (see [2] and citations therein) and QCD [3]. Whereas in d = 2 dimensions the CFTs are well studied and classified [4], for d > 2 both the classification and the tools for study of CFTs are notoriously incomplete. The supersymmetric QFTs, and in particular the supersymmetric Yang-Mills theories include a rather large class of CFTs which are relatively well classified and, at least qualitatively, understood [5][6][7][8][9][10][11], especially due to the discovery of the AdS/CFT correspondence. In rare cases, such as 4-dimensional N = 4 SYM theory or 3-dimensional ABJM theory, the integrability allows us to study in-depth, at least in the 't Hooft limit, the basic quantities of operator product expansion: all-loop anomalous dimensions [12][13][14][15], where the comprehensive and efficient solution is given by the quantum spectral curve (QSC) approach [16,17] (see also recent reviews [18,19]), OPE coefficients (structure constants) can be studied in various limits [20][21][22][23] and even obtain some non-perturbative information on multi-point correlators of local operators [24,25], cusped Wilson loops [26,27] and 1/N -corrections [28].
Much less is known about non-supersymmetric four-dimensional CFTs. Mostly, they are known to be IR or UV fixed points of various renormalization group flows. Usually these CFTs are strongly coupled at these fixed points and, apart from a few rather exotic cases, such as the Banks-Zaks CFT [29], their Lagrangian description is unknown. Such theories have been recently quite efficiently studied by conformal bootstrap methods [30,31] (using basic assumptions for CFTs, such as unitarity, crossing, symmetry etc.) which involve heavy numerical methods and can give very accurate predictions of OPE data. However, as for any numerical approaches, the physics behind these computations often remains obscure.
In this respect, various integrable deformations of N = 4 SYM, breaking partially or even entirely the supersymmetry [6,7,32,33], open a unique window into the dynamics of four-dimensional CFTs.
In particular, the γ-deformed N = 4 SYM, where breaking of the R-symmetry leads to the complete loss of supersymmetry, is a CFT with the well-defined classical action. In order for the theory to be consistent at the quantum level, one has to add to the action a finite number of particular scalar double-trace terms, for which the couplings have to be fine tuned to special values corresponding to the fixed points of the underlying beta-functions, or rather functions of the Yang-Mills coupling [34,35]. The only drawback of such a CFT is the loss of unitarity, since the double-trace couplings induced by the renormalization [36][37][38][39] take complex values at these fixed points. This poses an interesting challenge of construction of nontrivial unitary non-supersymmetric CFTs. On the positive side, the γ-deformed N = 4 SYM at the fixed point seems to be a genuine CFT, well defined by its explicit action, including the double-trace terms [35]. Last but not the least, quantum integrability property of planar γ-deformed N = 4 SYM [33] described in terms of γ-deformed quantum spectral curve formalism in [8], occurs precisely at these fixed points, as was conjectured and argued in [35]. The powerful machinery of quantum integrability allows us to study in great detail its complicated nonperturbative dynamics [40,41]. Moreover, in a specific double scaling limit proposed in [42], combining weak coupling with strong imaginary γ-twists, the γ-deformed N = 4 SYM drastically simplifies and gives rise to a family of chiral non-unitary CFTs with 3 effective couplings describing the scalar and Yukawa interactions of three complex scalars and three fermions. Its spectrum of anomalous dimensions, scalar and fermion amplitudes have been studied in a series of papers [35,[42][43][44][45][46]. In the particular, single coupling version of these models, the bi-scalar "fishnet" CFT, studied in this paper, the four-point correlation function of certain protected operators was computed in [35], providing a rich non-perturbative OPE data for the exchange operators with an arbitrary spin. These results have been generalized to any dimension d in [47], where the d-dimensional version of the bi-scalar model was proposed. In this paper we extend these results to more general correlators. In addition to the wheel graphs we also consider single and double spiral graphs as shown on Fig.1. We also analyse the results at weak and strong coupling.
An important feature of such models is the drastic simplification of their weak coupling expansion, where in many particular cases (when we turn on a single coupling) it is dominated by various kinds of "fishnet" Feynman graphs [43]. These graphs represent integrable two-dimensional statisticalmechanical systems by themselves [48] and can be efficiently studied by the quantum spin chain methods and the double-scaled version QSC [44,49]. In particular, the individual, so called "wheel" multi-loop Feynman graphs can be computed exactly in terms of multiple zeta values (MZV) [44].
Importantly, the explicit graph-by-graph integrability property in such models sheds some light on the origins of the planar integrability of their "mother"-theory -the N = 4 SYM, where the perturbation theory is much more complicated and the reasons for integrability are still obscure. In particular, in the bi-scalar CFT discussed in this paper the integrability is manifest due to the explicit integrability of fishnet graphs dominating the perturbation theory in this model [35].

The conformal "fishnet" theory
We will focus in this paper on a particular example of strongly γ-deformed N = 4 SYM -the bi-scalar theory [42]. At the classical level, the Lagrangian of the bi-scalar theory is given by 1 L = N c tr ∂ µX ∂ µ X + ∂ µZ ∂ µ Z + (4π) 2 ξ 2XZ XZ , where X, Z are complex N c × N c matrix fields andX ≡ X † ,Z ≡ Z † are their hermitian conjugates. The model retains the SU (N c ) global symmetry which is a remnant of the gauge symmetry of the original N = 4 SYM theory. The effective coupling constant ξ 2 = g 2 YM N c e −iγ3 /(4π) 2 is given by the product of the Yang-Mills coupling and the complex deformation parameter. The general γ−deformed N = 4 SYM theory depends on three deformation parameters γ 1 , γ 2 , γ 3 . The Lagrangian (1.1) arises in the double scaling limit, g 2 → 0 and Im γ 3 → ∞ with ξ 2 and γ 1,2 fixed. In this limit, all fields except two scalars get decoupled leading to (1.1).
The relative coefficients between the operators in both lines of (1.2) follow from the invariance of (1.1) under the transformations of fields with the conjugate fieldsX,Z transforming accordingly. As we show below, these transformations can be used to establish relations between different correlation functions. The theory with the Lagrangian L + L dt is renormalizable. The coupling constants depend on the renormalization scale and the corresponding beta functions have been computed perturbatively in [35] in the planar limit, for N c → ∞ and ξ 2 , α 2 1 , α 2 2 = fixed. Examining zeros of the beta functions, we find that the theory has two fixed points (α 2 1 = α 2 + , α 2 2 = ξ 2 ) and (α 2 where α 2 ± is given at weak coupling by the following expression . (1.5) Notice that the expansion of α 2 ± runs in powers of iξ 2 with real coefficients. The planar bi-scalar theory (1.1) possesses a conformal symmetry at the fixed points (1.4) [34,35] and is integrable [42,44]. Viewed as a function of ξ 2 , the relation (1.4) defines two lines of the fixed points. It has been argued in [35] that the "mother" theory of the bi-scalar model -the γ−deformed N = 4 SYM -is a nonunitary CFT on a line of (complex) fixed points of double couplings as functions of the 't Hooft coupling, also integrable in planar limit. It is also natural to expect the existence of such a complex conformal trajectory even at finite N .
The bi-scalar "fishnet" CFT (1.1) is the most studied case of the abovementioned chiral CFTs proposed in [41]. The spectrum of long local operators of the type tr(Z m X n )+permutations, can be efficiently investigated by the asymptotic Bethe ansatz (ABA) equation [43] -the double scaled version of the Beisert-Staudacher ABA equations for N = 4 SYM [12]. The short operators of this and other types (also with insertions of derivatives andZ,X fields) can be studied by QSC methods [44] and by the quantum non-compact spin chain methods [42,44], when the spins take values on conformal group SU (2, 2). The spin-chain approach to this theory is very promising since it would allow us to study there non-perturbative physics starting from the first principals, without any assumptions. Unfortunately, efficient methods of study of non-compact, Heisenberg spin chains are not very well developed, especially for principal series representation in physical space and for higher ranks symmetries, such as SU (2,2), though an important progress has been made in the study of spectral problem for SL(2, C) spin chain, in relation with high-energy (Regge) limit of QCD [50,51]. Another remarkable observation in bi-scalar theory concerns the planar scalar amplitudes: they are dominated by a single multiloop fishnet graph with open boundary and obey a Yangian symmetry, potentially allowing for their computation [45,46]. A particular, single-trace four-point correlation function given by such fishnet graph was explicitly computed in [52] We refer to the single trace operators tr(X n Z m ) + perm. for n ≥ m as m-magnon states, in accordance with the ABA description, where the asymptotic anomalous dimension is described by the Bethe state with m Bethe roots and conformal spin chain of length n. The related Feynman graphs have been described in [43]. They have a shape of multi-spiral when m radial lines of the field Z coming out of the point where the operator is placed, are "braided" by m parallel spirals, as shown on the Fig.1 (in the middle for a single spiral, and on the right for the double spiral).
Correspondingly, the simplest set of non-trivial single trace operators has length n = 2 and numbers of magnons m = 0, 1, 2. In addition, one can also introduce Lorentz spin S by inserting lightcone derivatives in the following way: tr(X n (∂ − ) S Z m ) + perm. The most efficient way of studying these operators is to extract their conformal data from the OPE of 4-point functions. Accordingly, we will analyze the 4-point functions of 3 different topologies, corresponding to the number of the magnons (see Fig.5). The simplest, zero-magnon four-point correlation function was computed to all orders of weak-coupling expansion in [35]. It is dominated by the wheel graphs containing only two "spikes". It was shown in [35] that this quantity has correct conformal properties in the perturbation theory only if one takes into account the double-trace interactions.
In the current paper, we will review the findings of [35] and continue to study the properties of the four-point correlation functions. In addition, we compute a few related four point-functions of short local operators, corresponding to one-and two-magnon cases. These three types of four-point functions are distinguished by the relative simplicity: for their computation one does not need to appeal to the integrability -the conformal symmetry is enough for this purpose. We then use the obtained results for the four point correlation functions to extract explicit expressions for the anomalous dimensions and structure constants.
All these results represent a unique opportunity of study properties of bi-scalar CFT, in the hope to better understand the non-perturbative structure of non-supersymmetric CFTs in d > 2 dimensions. They provide rich data for the future integrability based calculations of the correlation functions.

Correlation functions and their perturbative structure
In this paper, we exploit conformal symmetry to find exact expressions for correlation functions of local protected dimension-two operators as well as of bi-local operators of a "one-magnon" type The reason for the choice (1.6) is that, in the planar limit, the two-point correlation functions of operators (1.6) are protected at the fixed points (1.4) where the normalization factor c = 1/(4π 2 ) comes from free scalar propagator. The pair correlation function of bi-local operators of the type (1.7) defined below as type C will represent another type of four-point functions, having one-magnon exchange states in OPE, dominated by single-spiral graphs of the type shown in Fig. 1.
In what follows we consider the simplest unprotected four-point correlation functions of the local operators (1.6), of the following two types: • Type B The remaining four-point correlation functions of the operators (1.6) vanish due to nonzero total U (1) charge. Notice that type B correlation function expansion in small x 2 12 ≡ (x 1 − x 2 ) 2 limit is saturated by the two-magnon operators. Two such spin-zero operators, tr(XZXZ) and tr(XXZZ), are not protected and mix with each other, in such a way that their dimensions are related by the change ξ 2 → −ξ 2 . The analogous operator in N = 4 SYM theory is Konishi operator tr[X, Z] 2 , where as the second operator with the same R-charge, tr(2XXZZ +XZXZ), is protected 3 .
• Type C We will also define the type-C four-point functions containing one-magnon exchange states: (1.11) We can also define a similar pair of correlation functions 12) which, due to the relations (1.3), coincide with G C and G C , respectively.
• Type D In addition to (1.9) and (1.10), we also consider a four-point correlation function of scalar fields computed in [35] G D = tr(X(x 1 )X(x 2 )) tr(X(x 3 )X(x 4 )) . (1.13) As will become clear in a moment, its calculation is closely related to that of G A given by (1.9) (see (2.9) below). The function G D is obtained from G A by Wick contracting two pairs of Z andZ fields. Figure 5. Three basic types of the skeleton graphs corresponding to zero, one and two magnons. They have the structure of wheels, single and double spirals. These basics graphs are the building block of all correlators we discuss in this paper. Another way to represent the same graphs given on Fig.1. (1.14) For the correlation function (1.9), the relations (1.3) imply that G A should be invariant under the exchange of points A remarkable feature of all considered correlation functions is their iterative structure: the (nonzero) contribution at each successive order of perturbation theory can be obtained from the previous one by action on it by some graph generating integral operators. Thus the relevant graphs have a chain structure and they can be studied using the Bethe-Salpeter equation. In addition, the emerging graph generating operators commute with the generators of the conformal group and their eigenspectrum can be easily found with a help of the conformal symmetry. Thanks to these features the above mentioned correlation functions of the bi-scalar model are explicitly computable in a relatively simple way.

Relation to skeleton scalar graphs
In this section, we will describe the structure of Feynman graphs for all types of the studied 4-point correlation functions. We express these correlators in terms of the basic generating functions of the wheel graphs (G 0 ), single (G 1 ) and double spiral graphs (G 2 ) as on Fig.5. In the section 4 we evaluate these 3 types of the graphs by the Bethe-Salpeter method, by diagonalizing their graph generating kernels. Consequently, we compute all the related structure constants defining the full explicit OPE representation of each of these 4-point functions (1.11)-(1.13).
We will first discuss the connected and disconnected parts of all these 4-point functions and then study the connected planar part of each of them. For the sake of explicitness, we will restrict our discussion to d = 4 dimensions, but the final formulas will be readily generalized in section 9 to any d.

Connected part of the correlation functions
The four-point correlation functions of the types A and B (1.9) -(1.10) receive both connected G (c) and disconnected (in the sense of factorisation of the coordinate dependence) contributions G (d) . The former are suppressed with respect to the latter by a power of 1/N 2 but all the most interesting physics resides of course in the connected part. The disconnected contribution is given by the product of two-point correlation functions (1.8) leading to For the correlation function of type C and D the disconnected contribution is Due to the color structure of bi-local operators in the definition of G C and G D , the disconnected part in both cases is of the same order in 1/N c as the connected part. Finally, the disconnected part of the correlation function of type C is suppressed by the factor of 1/N 2 c as compared with its connected part.

Relation of 4-point functions of all types to basic fishnet graphs
The computation of all above-mentioned types of 4-point functions can be reduced to the evaluating Feynman graphs having 3 basic structures of fishnet graphs. They are depicted on the Fig.1 and they are distinguished by the the magnon numbers -0, 1, 2 -the number of propagating "particles" (dotted spirals lines on the Fig.1) in the exchange channel. We will denote the related generating functions as G 0 , G 1 and G 2 and call them the n-magnon functions. As we will see shortly, the Feynman graphs for these functions are computable by the Bethe-Salpeter approach, due to their periodic fishnet structure and conformal properties.
To be more precise, we will define n-magnon functions as perturbative expansion w.r.t. the coupling ξ 2 as follows.
For zero-magnon case (the "wheels") we have where G For one-magnon case ("spiral", or "spiderweb" graphs) the structure function looks as where G (n) 1 (x 1 , x 2 |x 3 , x 4 ) denotes the Feynman graph depicted in the middle of Figs.1 and 5 and having n interaction vertices (black dots). Each graph takes a shape of a spiral consisting of propagators of type X wounding around two lines of propagators of type Z. Note that we have two distinguished structures depending on the parity of expansion term w.r.t. ξ 2 : for even powers of ξ 2 the spiral on Fig.1 starts and ends on the same (black) line of Z-propagators, whether as for odd powers of ξ 2 it starts and ends on different Z-lines.
For two-magnon case ("double spiral") we have the structure function where G (n) 2 (x 1 , x 2 |x 3 , x 4 ) denotes the Feynman graph depicted on the right of Figs.1 and 5 and having 2n interaction vertices (black dots).
Let us relate the correlation functions G A , G B , G C and G D to G 0 , G 1 and G 2 .

Relations between correlation functions and n-magnon functions
First, we can immediately see from Figs. 1(left) and 5(left) and from the definitions (1.13) and (2.4) that Further on, the correlation function G A is given by the following linear combination of the functions G 0 where the last term takes care of double counting of the tree-level diagram in the first 4 terms. It is straightforward to verify that the linear combination on the right-hand side of (2.8) satisfies the symmetry relation (1.15). Comparing (2.8) and (2.7) we notice that the two correlation functions are related to each other as (2.9) Having determined G 0 (x 1 , x 2 |x 3 , x 4 ), we can apply (2.7) and (2.8) to find the correlation functions G A and G D . For the one-magnon correlation functions we get the following expressions through the even and odd in ξ 2 parts of the magnon function G 1 : Finally, the two-magnon correlation function (1.10) coincides with the two-magnon function (2.6): In the next section we review the general method for computing the magnon functions G 0 , G 1 and G 2 based on the Bethe-Salpeter equation and conformal symmetry. The explicit expressions for these functions are derived in section 4.

Conformal symmetry and Bethe-Salpeter equations: generalities
As we saw above and discuss in detail in the next section, the correlation functions introduced in previous sections are given by very specific type of fishnet graphs: in each graph the periodically repeating configurations of propagators are connected by pairs of coordinates of the related vertices (see Fig. 5). Each graph can be cut into two disconnected parts by splitting only two vertices. The three cases we are going to consider differ, in particular, by the values of dimensions of four external (protected) operators. For all correlation functions under considerations we have ∆ 1 = ∆ 4 and This section is based on the observation that three topologically distinct configurations G 0 , G 1 and G 2 can be written, each, in terms of a suitable "graph-building" operatorĤ. In each case we find at the level of the operatorŝ where n is a nonnegative integer, the constant χ is proportional to a fixed power of the coupling constant ξ 2 (specified below for each case) and d is the dimension of the space-time. In most of the paper we set d = 4 although, as we will see in the section 9, most of the equations discussed here have a natural generalization to general d, where the bi-scalar theory can be also formulated [47]. The operatorsĤ andĜ are represented by the corresponding integration kernels, e.g.
in such a way that  The problem of finding G(x 1 , x 2 |x 3 , x 4 ) can be split into two main steps. First, we have to solve the eigenvalue problem forĤ and, then, decomposeĜ over the complete basis of the eigenfunctions ofĤ. Fortunately, the first step is simple in our case. The eigenfunctions Φ µ1,...,µ S ν,x0 (x 1 , x 2 ) ofĤ which are defined by where x 0 , ν and S parameterise the eigenstates, are totally symmetric traceless tensors in 4-dimensional indices µ 1 , . . . , µ S . The form of Φ µ1,...,µ S ν,x0 (x 1 , x 2 ) is completely fixed by the conformal symmetry -if the operatorĤ commutes with the generators of the conformal group, its eigenstates should transform covariantly under conformal transformations acting on x i . As such, Φ µ1,...,µ S ν,x0 (x 1 , x 2 ) can be represented as a conformal "triangle" -three-point correlation function of two scalar operators at the points x 1 , x 2 and some reference operator with dimension ∆ = 2 + 2iν and spin S at the point x 0 . Explicitly where 4 t = ∆ − S and x ∆ ab ≡ (x 2 ab ) ∆/2 . In order to simplify tensor structure, we projected all Lorentz indices onto a light-like vector n µ . In the next section we verify explicitly that the functions (3.5) diagonalize the graph-building Hamiltonians and find the corresponding eigenvalues E ∆,S .
The scaling dimension ∆ in (3.5) is given by [53] ∆ = 2 + 2iν (3.6) where ν is real nonnegative. For such values of ∆ the functions Φ µ1,...,µ S ν,x0 (x 1 , x 2 ) define the complete orthonormal basis of states on the Hilbert space on which the graph building kernelĤ acts. Viewed as a function of x 0 , Φ µ1,...,µ S ν,x0 belongs to the irreducible principle series representation of the conformal group labelled by real ν and nonnegative integer Lorentz spin S. Together with (3.4) this leads to where c 1 is the normalization factor defined in (A.4) and we used thatΦ µ1...µ S ν,x0 . The values of ν in (3.7) can be restricted to be ν ≥ 0 since the states Φ µ1,...,µ S ν,x0 and Φ µ1,...,µ S −ν,x0 belong to the same representation and are related to each other through intertwining relation.
Inserting (3.7) into (3.1) we get the following expression for the correlation function G (3.8) The integral over x 0 can be evaluated explicitly in terms of four-dimensional conformal blocks [53][54][55] the expression in the brackets depends only on the cross ratios u and v defined as where z andz are auxiliary complex variables. The conformal block g ∆,S (u, v) depends on ∆ 1 and ∆ 2 and it is given explicitly in terms of the hypergeometric functions in (A.1). The coefficient c(ν, S) is a ratio of the normalisation coefficients c(ν, S) = c 1 (ν, S)/c 2 (ν, S) defined in (A.4) and (A.6). Combining together (3.8) and (3.9) we obtain 5 where G(u, v) admits the following representation 4 Following the standard conventions, we shall refer to t at zero value of the coupling constant as twist. 5 The constant c = 1 4π should not be confused with the function c(ν, S).
where ∆ = 2 + 2iν. Here we combined the two terms on the right-hand side of (3.9) to extend the integration over ν to the whole real axis and used the identity c(ν, S) = c 1 (ν, S)/c 2 (ν, S). Note that doing so we required that This property follows from the fact that the states with the same Lorentz spin S and the scaling dimensions ∆ and 4 − ∆ belong to the same representation of the conformal group. We will check (3.13) explicitly in each case.
To bring the integral (3.12) to the standard OPE form we examine the short distance limit x 12 → 0, or equivalently u → 0 and v → 1. In this limit, the conformal block scales as g 2+2iν,S (u, v) ∼ u 1+iν−S/2 (1−v) S and decays exponentially fast for Re(iν) → ∞. This allows us to close the integration contour in (3.12) to the lower half-plane and compute the integral on the right-hand side of (3.12) by residues. In terms of the OPE, the condition Re(iν) > 0 is equivalent to the restriction Re ∆ > 2 on the scaling dimension of exchanged operators.
The integrand in (3.12) has 'physical' poles coming from zeros of the denominator and two series of 'spurious' poles generated by the kinematical factor c 2 (ν, S) and the conformal block g 2+2iν,S (u, v). We show in Appendix B that the spurious poles cancel against each other provided that E ∆,S satisfies the following relation with r k defined in (B.1). Then, the correlation function (3.12) is given by the sum of residues at the physical poles (3.14). Finally, we obtain the following conformal partial wave expansion of the correlation function (3.12) where the OPE coefficients are given by the residues at the physical poles (3.17) and the sum in (3.16) runs over solutions of (3.14) with Re ∆ > 2.
In the next section, we apply the relations (3.17) and (3.16) to compute the four-point correlations introduced in the previous section. In each case we shall verify the relations (3.13) and (3.15) which we assumed in the above derivation.

Four-point correlators and the conformal OPE data
In this section, we describe Feynman graphs for the 3 types of 4-point correlation functions depicted on Fig.5, and establish the corresponding graph-building operatorsĤ. In the previous section, it was shown that the eigenvalue E ∆,S of these operators is the only input needed in order to write a representation of OPE type for the correlation function. We diagonalize the operatorsĤ and present explicit expressions for the conformal data (scaling dimensions and the OPE coefficients) for each of these 4-point functions. Figure 6. First 3 orders contributing to the G0 correlator.

Zero-magnon case and the wheel-graphs (G 0 )
The zero-magnon correlation function was studied in detail in [35]. Here we re-derive the results of [35] and show how they fit into the general scheme described in the previous section.
As we show below, the zero-magnon case corresponds to the situation when the correlation function (3.16) receives contribution from magnon-free operators of the type tr(X(n∂) S X) and tr(X (n∂) S X) 6 . We will see that only these two types of operators with arbitrary spin S contribute.
To find the zero-magnon correlation function we have to summed up diagrams shown on Figure 5(left). These diagrams contain an arbitrary number of wheels attached to the rest of the diagram at two (single-trace) vertices. It is easy to see that the integral over the position of these vertices develops a ultraviolet (UV) divergence at short distances. This seems to be in contradiction with expected UV finiteness of four-point correlation function of protected operators. We recall however that quantum corrections induce new double-trace interaction terms (1.2). In particular, the doubletrace coupling tr(X 2 ) tr(X 2 ) affects the zero-magnon correlation function. It produces a UV divergent contribution which cancels against that of the wheel graphs in such a way that the four-point correlation function remains UV finite at any order in the weak coupling expansion. Due to form of the double-trace interaction term tr(X 2 ) tr(X 2 ), it can only affect the contribution of partial waves to (3.16) with zero Lorentz spin S = 0. We therefore expect that the contribution of the wheel graphs to (3.16) is well-defined for S = 0 whereas for S = 0 the additional contribution due to double traces should be taken into account. We discuss this issue in more detail in Sec.5.
In this section we proceed without taking the double trace interaction into account and identify the contrubution of wheels graphs to (3.16). We will see that the double trace contributions will be automatically taken into account by correct treatment of the singularity at ξ → 0 in the forthcoming formulas. We start with constructing the graph building operator and identifying the parameters ∆ i , n and χ introduced in (3.1). We recall that the parameters ∆ i define the scaling dimension of the external operators. Since the wheel graphs have only one propagator attached to each external leg we have ∆ 1 = ∆ 4 = ∆ 2 = ∆ 3 = 1.
To first two orders of the weak coupling expansion the zero magnon function G 0 is given by the 6 Interestingly, in our case only single box can appear. In space-time dimension d = 4 this is not the case.
sum of diagrams shown in Fig.6. The expressions corresponding to the first two diagrams are where each scalar propagator brings in the factor of c/x 2 ij . These expressions can be represented as 1st and 2nd powers of the following graph building operatorĤ 0 Indeed, we verify that from where it is clear that for a general graph with (n + 1) wheels we getĜ 0 It is straightforward to verify that the function (4.2) transforms covariantly under the conformal transformations acting on x i . 7 As a consequence, the corresponding integral operatorĤ 0 commutes with the generators of the conformal group.
Thus the zero-magnon correlation function G 0 can be written aŝ Comparing with the general expression (3.1) we deduce that χ = (16π 2 ξ 2 ) 2 and n = 1 in the zeromagnon case.

Eigenvalue of the zero-magnon graph-building operator
In order to use the general expression for the correlation function (3.16), we have to determine the eigenspectrum of the graph building operator (4.2). In virtue of conformal symmetry, its eigenstates are given in (3.5) with ∆ 1 = 1 and ∆ 2 = 1. Substitution of (4.2) into (3.4) leads to an integral, which can be evaluated using the star-triangle identity as explained in appendix C.2. Going through the calculation we obtain the following simple result [35] It is easy to see that E 0 is invariant under ∆ → 4 − ∆, in agreement with (3.13). We can also check that (4.5) verifies the relation (3.15) that ensures the cancellation of spurious poles. In the present case, for ∆ 1 = ∆ 2 = 1, it follows from (B.1) that r 2n+1 = 0 for integer n and thus (3.15) reduces to which is indeed satisfied for (4.5).

Spectrum for zero-magnon exchange operators
We can now use (4.5) to determine the scaling dimension of the zero-magnon operators contributing to the correlation function (3.16). Substituting (4.5) into the relation (3.14) and replacing χ = (16π 2 ξ 2 ) 2 we find that the scaling dimensions ∆ = 2 + 2iν satisfy the following quartic equation The simplest way to check this is to employ inversions x µ i → x µ i /x 2 i and take into account that subject to the additional condition Im ν < 0. At finite coupling, this yields the following expressions for the scaling dimensions ∆ 2 (S) = 2 + (S + 1) 2 + 1 − 2 (S + 1) 2 + 4ξ 4 , ∆ 4 (S) = 2 + (S + 1) 2 + 1 + 2 (S + 1) 2 + 4ξ 4 , (4.8) The two remaining solutions to (4.7) are related to (4.8) by ∆ → 4 − ∆ and describe shadow operators with Re∆ < 2. At weak coupling, for ξ 2 < 1, and nonzero Lorentz spin, S > 0, the scaling dimensions (4.8) look as and the corresponding operators can be identified as twist-2 and twist-4 operators, respectively. 8 They have the following form O 2 = tr(X(n∂) S X + . . . and O 4 = tr( X(n∂) S X + . . . where dots denote similar terms with light-cone derivatives distributed between the fields. Similarly to [44], the twist 4 operators can be written, due to the equations of motion, as O 4 = tr X∂ SZ XZ + . . . . Notice that the weak coupling expansion of (4.9) goes in powers of ξ 4 which is exactly what one expects since each wheel in the graph shown in Fig. 5(left) is attached to the rest of the diagram through two single-trace vertices. Something special happens at S = 0. In this case, we find from (4.8) (4.10) and the weak coupling expansion looks as Surprisingly, for S = 0 the weak coupling expansion of the scaling dimension ∆ 2 (S) starts from the power ξ 2 , instead of the naively expected ξ 4 (the power corresponding to each insertion of the graph-building operator (4.2)).
To understand the reason for this we examine the eigenvalue of the zero-magnon kernel (4.5) for S = 0 and ∆ = 2 + 2iν We notice that it goes to infinity at small ν. Then, expanding (4.4) in powers of ξ 4 we find that the contribution of the states with S = 0 to the correlation function at order O(ξ 4n ) is proportional to dν(E 0 | S=0 ) n and it diverges for ν → 0. This is in agreement with our expectations that the contribution of the wheel graphs is well defined for all states except those with S = 0. To remove the divergence we have to include the O(ξ 4n ) contribution of double traces. We observe that, in the resummed expression for the correlation function (4.4), the contribution of the states with S = 0 to G 0 involves the integral dν/(1 − (16π 2 ξ 2 ) 2 E 0 | S=0 ) which is convergent for ν → 0 at finite ξ 2 (for the integral to be well-defined we assume that ξ 2 has a small imaginary part). It is easy to see that, at weak coupling, the integration over small ν produces a square-root singularity at the origin, G 0 ∼ ξ 4 . This explains why the weak coupling expansion of G 0 in powers of ξ 4 is divergent. At the same time, this also suggests that the double-trace contribution should be essential only in the weak coupling regime whereas at finite coupling it can be safely ignored. We discuss this issue in more detail in Sect. 5.
Arriving at (3.16) we have tacitly assumed that the physical poles (3.14) are located away from real ν−axis in (3.12). As follows from (4.7), at weak coupling, the two physical poles located at ν = ±iS/2 + O(ξ 4 ) pinch the integration contour at the origin for S → 0 and produce a divergent contribution. The role of the double-trace contribution is to subtract this divergence and, thus, make the weak coupling expansion of G 0 well defined. Turning the logic around, we can say that the double traces provide a nonvanishing contribution to the scaling dimensions because the eigenvalue (4.5) diverges as E 0 (ν, 0) ∼ 1/ν 2 for ν → 0. The relation (4.10) is in a perfect agreement with the result of explicit 7−loop calculation [35].
We recall that the correlation function (3.16) is given by the sum over the solutions to (3.14) with Re ∆ > 2. In our present case for S = 0 the solution (4.10) satisfy Re ∆ = 2 for real ξ 2 . As was mentioned above, for the correlation function (3.16) to be well-defined ξ 2 should have a nonvanishing imaginary part (see Appendix D for discussion of analytical properties of (3.16)). The expression (4.10) satisfies Re ∆ = 2 for Im ξ 2 < 0. For Im ξ 2 > 0, the scaling dimension is given by the same expression (4.10) upon replacing ξ 2 → −ξ 2 .
Let us examine the properties of the scaling dimensions (4.8). The dependence of ∆ 2 and ∆ 4 on ξ 4 is shown on Fig. 7 for S = 0 and S = 2. We observe that the two functions (4.8) represent in fact two branches of the same analytic function. It has two branch points located at For ξ 4 = ξ 4 − the two operators collide, ∆ 2 (S) = ∆ 4 (S), whereas for ξ 4 = ξ 4 + one of the operators collides with its shadow, ∆ 2 (S) = 2. 9 . The collision of operators at ξ 4 = ξ 4 + modifies analytic properties of the correlation function G 0 . According to (D.2), the correlation function has the cut for ξ 4 > 0 that starts at ξ 4 = 1/max ν∈R E 0 (ν, S). It is easy to see from (4.5) that E 0 (ν, S) is a decreasing positive-definite function of ν 2 . Therefore, the cut starts at ξ 4 = ξ 4 + (corresponding to ν = 0) and goes to infinity along real ξ 4 −axis. In the vicinity of the branch point, it follows from (4.7) that Strong coupling. At strong coupling, for ξ 2 → ∞, the relation (4.7) has four solutions iν = ξ e iπk/2 +O(1/ξ) (with k = 0, . . . , 3). Among them only two satisfy the additional condition Im ν < 0. The corresponding expressions for the scaling dimensions are where integer 0 ≤ k ≤ 3 satisfies the condition Re(ξ e iπk/2 ) > 0 and depends on ξ.

Structure constants for zero-magnon exchange operators
We apply the general relation (3.17) to find the OPE coefficient for zero-magnon operators [35] .
Stricktly speaking, C ∆,S is given by the product of (properly normalized) 3−point correlation functions OX Z OXZO ∆,S and O XZ O XZŌ∆,S . In unitary CFT they are complex conjugated to each other and, as a consequence, C ∆,S is positive definite. This is not the case for the conformal fishnet theory (1.1) and (1.2). In virtue of the symmetry (1.3) the above mentioned 3-point functions coincide and, therefore, C ∆,S is given by the square of the 3-point correlation function of two protected and one unprotected operators We will generalize this result to a more complicated structure constant involving two non-protected operators in section 8.
Weak coupling limit. Replacing the scaling dimensions in (4.15) by their expressions (4.8) we can obtain the OPE coefficients at weak coupling. First, consider twist-2 operators with the scaling dimension ∆ 2 and non-zero spin S > 0. In this case we get where ψ(x) is the Euler polygamma-function. Notice that C ∆ t=2,S becomes singular for S → 0. Similar to the situation with the scaling dimension ∆ 2 (S), this happens because the two limits S → 0 and ξ 2 → 0 do not commute. To get the correct result for C ∆ t=2,S=0 at weak coupling, we should first put S = 0 in (4.15) and, then, expand it in powers of ξ 2 . This gives In analogy with (4.11), the weak coupling expansion starts at order O(ξ 2 ) indicating that C ∆ t=2,S=0 is sensitive to the contribution of the double traces. For twist−4 operators we find from (4.15) and (4.9) In distinction with (4.18), the weak coupling expansion of C ∆t=4,S starts at order O(ξ 4 ) and runs in powers of ξ 4 . The latter property is in agreement with our expectations that twist−4 operators are not affected by double-trace interaction. The twist-4 OPE coefficient (4.19) is suppressed by the factor of ξ 4 as compared with (4.18). Due to the equation of motion, X = 16π 2 ξ 2Z XZ, the corresponding operator takes the form O 4 = tr( X(n∂) S X) + · · · = 16π 2 ξ 2 tr(ZXZ(n∂) S X)) + . . . . The reason for the above-mentioned suppression is that Strong coupling limit. At strong coupling the dimension ∆ become large according to (4.14) and we get Thus we see that the structure constant at strong coupling is exponentially small since ∆ 2ξ e iπk/2 .

Zero-magnon 4-point correlation function
Having determined the conformal data of the zero-magnon operators, we can apply (3.11) and (3.16) to compute the four-point correlation function (3.11) and (3.16) where we replace the scaling dimensions of the external protected operators by their values, ∆ 1 = ∆ 2 = ∆ 3 = ∆ 4 = 1, and the function G 0 (u, v) is given by where the sum runs over the states with the scaling dimension (4.8) and the OPE coefficients (4.15).
Here g ∆,S is the four-dimensional conformal block defined in (A.1) (with ∆ i = 1). The relation (4.22) involves an infinite sum over the conformal blocks and it is not obvious a priori that one can find a closed expression for G 0 (u, v) even at weak coupling. We show in section 5 by explicit two-loop calculation that G 0 (u, v) can be expressed in terms of special functions, the so-called harmonic polylogarithms (HPL). In section 6 we extend this result to any order of the weak coupling expansion. Also in section 7 we analyse the strong coupling limit of the expression (4.22). The analytic properties of G 0 (u, v) with respect to the coupling ξ are discussed in Appendix D.
This relation follows from (3.9), it can also be verified directly from the definition (A.1). Combining together (4.22) and (4.23) we conclude that, in the expressions for G A and G D the terms in (4.22) with odd S cancel out whereas those with even S get doubled.

One-magnon case and single spiral graphs (G 1 )
In this subsection, we consider the one-magnon correlation function G 1 described by graphs shown in Fig. 5(middle). As we will see shortly, its calculation is simpler than that of G 0 and G 2 . Since the graph contributing to G 1 have two propagators attached to x 1 and x 4 and only one to x 2 and x 3 , we identify the scaling dimensions at the external points as To identify the graph-building operatorĤ 0 , we consider the first few terms in the weak coupling expansion of G 1 (see Fig.8). The expressions corresponding to the first two diagrams are Let us show that these expressions can be represented as 2 nd and 3 rd power of the graph building operator with the integral kernel Indeed, we apply (3.3) to get It is clear from these examples that, in general,Ĝ 1 . Thus the correlation function G 1 can be written asĜ Comparing with the general form (3.1) we see that χ = 16π 2 ξ 2 and n = 2.

Eigenvalue of the one-magnon graph-building operator
As in the previous case, we can verify that the integral operatorĤ 1 with kernel given by (4.26) commutes with the generators of the conformal group acting on the external points x i and the corresponding conformal weights given by (4.24). As a consequence, its eigenstates are given by (3.5) with ∆ 1 = 2 and ∆ 2 = 1. The eigenvalue equation (3.4) reduces to an integral which can be evaluated using the star-triangle identity, as explained in Appendix C.3. The result turns out to be quite simple .

Spectrum of one-magnon exchange operators
Substituting (4.29) into (3.14) and taking into account that χ = 16π 2 ξ 2 , we determine the scaling dimension of the one-magnon operators O 1 = tr(XZ(n∂) S X) + . . . contributing to the correlation function (3.16) Interestingly, this expression coincides with the asymptotic dispersion relation for the one-magnon state previously found from the double-scaling limit of asymptotic Bethe ansatz (ABA) in Ref. [43]. Following the ABA approach, one could have expected that the asymptotic dispersion relation (4.32) should be corrected already at order O(ξ 4 ) by the wrapping corrections. The relation (4.32) implies that the wrapping corrections to the one-magnon operator tr(X 2 Z) vanish. This is indeed the case for the single wrapping contribution found in [43]. The relation (4.32) is indeed consistent with the ABA, including the known single wrapping correction! 10 Note that the factor (−1) S in the expression for E 1 plays an important role for this equation to be satisfied. Weak coupling. At weak coupling, the scaling dimension (4.31) of the one-magnon states reads Strong coupling. At strong coupling, we find from (4.31) where the branch of (−1) S+1 2 is such that Re ∆ > 2. Interestingly, at the leading order we get the same coefficient as for the zero-magnon case (4.14).

Structure constants with one-magnon exchange operators
We can also employ the general equation (3.17) to find the OPE coefficient .

(4.35)
As in the previous case, it is given by the square of the 3−point function of two protected and one unprotected operator, Weak coupling. At weak coupling, we obtain from (4.35) and (4.33) In particular for S = 0 we find Since the first few coefficients do not involve wrapping it should be possible to compare with the perturbative ABA based general expressions from [57,58].
Strong coupling. At strong coupling ∆ becomes large (4.34) and we can use it as a large expansion parameter since the expressions are more compact The structure constant is again exponentially decaying as it was in the zero-magnon case (4.20).

One-magnon four-point correlation function
Substituting (4.24) into (3.11) and (3.16), we obtain the one-magnon 4-point correlation function where and the sum runs over the states with the scaling dimensions ∆ and the OPE coefficients C ∆,S given by (4.31) and (4.35), respectively. Notice that in the above sum there is only one state for each value of the spin S and spin takes all non-negative integer values. In section 6 we study (4.40) at weak coupling and compare it with the result of perturbative calculation performed in section 5. Also in section 7 we analyse the strong coupling limit of (4.40). Finally, we can use (2.10) and (4.39) to calculate the 4-point correlation functions G C and G C .

Two-magnon case and double-spiral graphs (G 2 )
The two-magnon correlation function G 2 is given by graphs shown in Fig. 5(right). Since these graphs have two propagators attached to all four external points, we identify the scaling dimensions as As before, in order to construct the graph building operator for the two-magnon caseĤ 2 , we examine the first few orders of the weak coupling expansion of G 2 (see Fig.10). The expressions corresponding to the first two diagrams on Fig.10 are (4.42) Note that the two-loop integral entering G 2 factorizes into a product of one-loop integrals. This is not the case already at the next O(ξ 8 ) order for the right-most diagram in Fig.10. Figure 10. First 3 orders contributing to the G2 correlator.
The kernel H 2 generating the two-magnon diagrams is . (4.43) Indeed we verify that the convolution of H 2 reproduces all the diagrams depicted in Fig.10.
Thus for the sum of all two-magnon diagrams we get Comparison with the general expression (3.1) shows that χ = (16π 2 ξ 2 ) 2 and n = 1.

Eigenvalue of the graph-building operator
To compute the two-magnon correlation function (4.45) we have to diagonalize the operatorĤ 2 . We can use (4.43) to show that it commutes with the generators of the conformal group. As a consequence, its eigenstates are given by (3.5) with ∆ 1 = ∆ 2 = 2. To find the corresponding eigenvalue E 2 , we replace the eigenstates in (3.4) by their explicit expressions (3.5). This leads to a rather complicated integral on the left-hand side of (3.4). We can simply its calculation by sending x 0 → ∞ on the both sides of (3.4). In addition, we project all Lorentz indices on an auxiliary light-like vector n ν (with n 2 = 0) and obtain the following representation for E 2 where ∆ = 2 + 2iν and we put x 2 12 = (nx 12 ) = 1 for convenience. Since the integrand of (4.46) acquires the (−1) S factor under the exchange of the integration points, x 3 ↔ x 4 , E 2 vanishes for odd S. For even S the calculation of (4.46) yields (see appendix C.4 for details) where ψ (1) (x) = dψ(x)/dx. The expression for E 2 is manifestly invariant under ∆ → 4 − ∆. Let us verify the condition (3.15) for cancelling the spurious poles. In the present case, for ∆ 1 = 2 and ∆ 2 = 2, it follows from (B.1) that r 2n+1 = 0 for integer n, and the relation (3.15) reduces to for n, s = 0, 1, 2, . . . . It is easy to check that it is indeed satisfied.

Spectrum
The scaling dimensions of the two-magnon operators satisfy the relation (4.49) subject to Re ∆ > 2 and S being even nonnegative. This time the spectrum of ∆'s has a rich structure since for each value of S there are infinitely many solutions to (4.49). Indeed, as follows from the first relation in (4.47), the function (4.47) has an infinite sequence of double poles at ∆ = S + t for t = 4, 6, 8 . . . . As a consequence, for small ξ the relation (4.49) always has two solutions in the vicinity of ∆ = S + t describing operators with twist t and even spin S. 11 Indeed we see on Fig.12 that all levels are twice degenerate at weak coupling. The weak coupling expansion of ∆ 2 can be found by replacing E 2 by its expansion in the vicinity of the pole. Going through the calculation we obtain (4.51) For each t and S, the relations (4.50) and (4.51) yield two scaling dimensions that are related to each other through the transformation ξ 2 → −ξ 2 . For even t/2, the expansion coefficients in (4.50) are real. For odd t/2, the leading coefficient γ t,S is pure imaginary (see Fig.12). The relations (4.50) and (4.51) describe the scaling dimensions of an infinite set of operators, two per each twist t = 4, 6, . . . and spin S = 0, 2, . . . .
In particular, for S = 0 and t = 4, for the two-magnon operators of the form tr(X 2 Z 2 ) + . . . the scaling dimensions are given at weak coupling by The first 4 terms reproduce the prediction from ABA [59] including the first Lüscher correction.
Critical coupling. The dependence of the two-magnon scaling dimensions ∆ 2 on the coupling constant is shown on Fig.12. As can be seen from this figure, for each S the lowest level approaches the value ∆ = 2 at some finite ξ = ξ * . Expanding E 2 near ∆ = 2 we find the corresponding value of the coupling constant ξ * For ξ > ξ * , the scaling dimension ∆ develops an imaginary part. As we see in a moment, it grows linearly with ξ at strong coupling.
Notice that the real part of most of the zeros is close to integer. They determine the leading large ξ asymptotics of the scaling dimension of all but two states shown in Fig.12.
The remaining two states satisfy Re ∆ = 2 and have an imaginary part that grows linearly with ξ. They correspond to the solution to (4.49) with ∆ → ∞. Indeed, the function (4.47) decays at large ∆ as E 2 ∼ 1/∆ 4 and the relation (4.49) is automatically satisfied. By expanding E 2 at large ∆ and solving (4.49) we find (4.55)

OPE coefficients
From (3.17) we get for C ∆,S : where c 2 is given by (A.6) for ∆ 1 = ∆ 2 = 2. Replacing E 2 with (4.47), we obtain a rather cumbersome expression for C ∆,S , we do not present it here.
Here γ (0) n,S is the one-loop anomalous dimension defined in (4.51) and the coefficient c n,S is given by . Strong Coupling. To analyse the strong coupling behaviour of (4.56) it is convenient to rewrite it as C ∆,S = − π/c 4 c 2 (∆, S) ∂∆ ∂(4πξ) 4 (4.59) where we used (4.49) to get rid of E 2 . Since for most of the states ∆ approaches a constant value at strong coupling, we see from (4.59) that C ∆,S should decay as 1/ξ 8 . For the two remaining states whose imaginary part grows linearly in ξ, the OPE coefficient (4.59) decreases slower as 1/ξ for real ξ. We recall however that in order for the correlation function (3.12) to be well-defined, ξ should have an imaginary part. For ξ with large imaginary part the structure constant of the lowest level at each S decays exponentially and thus the OPE expansion at strong coupling should be dominated by the other state, for which ∆ → const at ξ → ∞. This is rather different behaviour to that of two other correlators and needs further investigation. We discuss this issue briefly in section 7.

4-point correlation function
Having determined the scaling dimensions and the OPE coefficients we can compute the two-magnon correlation function where We note that for each spin S and twist t there are two states, as one can see from the weak coupling expansion (4.50). Notice that the weak coupling expansion of (4.50) and (4.57) goes in powers of ξ 2 . However, due to the symmetry of the spectrum, the two terms in the sum (4.61) are related to each other through transformation ξ 2 → −ξ 2 so that the weak expansion of G 2 (u, v) runs in powers of ξ 4 .
In section 6, we study the equation (4.61) at weak coupling and compare the result with the predictions from the perturbation theory of section 5.

Correlation functions at weak coupling from Feynman diagrams
In the previous sections, we have derived three different types of four-point correlation functions by applying the operatorial methods. Namely, we have solved the underlying Bethe-Salpeter equations by diagonalizing the corresponding "graph-building" kernels with a help of the conformal symmetry. Doing so, we have ignored the double trace counterterms (1.2) which are nessesary for the consistent definition of the bi-scalar model (1.1) on the quantum level and for restoring the conformal symmetry of the theory.
In this section we discuss the role of the double-trace interaction terms (1.2). We show that they are necessary at weak coupling in order to make each order of the perturbation expansion of the correlation functions to be finite. At the same time, they do not affect the results for the correlation functions at finite coupling obtained in the previous section.
Let us review the role of double-trace couplings α 1 and α 2 from the action (1.2) in perturbative computations of the correlation functions G A , G B , G C and G D . Below we discuss which of the topologies of the Feynman graphs G 0 , G 1 or G 2 have to be completed with the double-trace interactions in order to have meaningful weak coupling expansion of the abovementioned four-point functions.

Double-trace contribution to G A , G B
We recall that the four-point functions G A , G B are completely defined, at least for any finite ξ, by the zero-magnon function G 0 . The latter is given by sum over the wheel graphs shown in Fig. 1. As was already mentioned, each wheel in these graphs develops a ultraviolet divergence at short distances. We expect that the double-trace contribution should remove this divergence. The double-trace interaction is described by the action (1.2). It is easy to see that, due to the cylindrical topology of the underlying planar graphs, among four different double-trace interaction terms in (1.2) only one term (4π) 2 α 2 1 tr(X 2 ) tr(X 2 ) can contribute to G 0 in the planar limit. It generates a new local quartic scalar vertex. The resulting planar graphs contributing to G 0 are shown in Fig.13. They are obtained by gluing together wheel graphs. Indeed, the insertion of the double-trace vertex 16π 2 α 2 1 tr(X) 2 tr X 2 effectively splits the planar wheel Feynman graph into two disconnected parts, with the single trace operators, tr(X) 2 and tr X 2 , attached to each part.
As we demonstrated in the previous section, the wheel graphs can be summed up by introducing the graph building operatorsĤ 0 . In the similar manner, we can take into account the graphs shown in Fig. 13 by replacing the kernel (16π 2 ξ 2 ) 2Ĥ 0 in the equation (4.22) by a linear combination (16π 2 α 2 1 ) V + (16π 2 ξ 2 ) 2Ĥ 0 of operators V andĤ 0 generating double-trace vertices and scalar loops, respectively (see Fig. 13). Since the contribution to the correlation function of individual diagram shown in Fig. 13 is divergent, we introduce dimensional regularization with d = 4 − 2 . Then, the regularized operators H 0 and V are defined as where Φ(x 1 , x 2 ) is a test function. They admit a simple diagrammatic representation, see Fig. 14(right) and (left), respectively. We would like to emphasize that the operators (5.1) are well-defined for = 0. Making use of the operators (5.1) we obtain the following representation for the zero-magnon correlation function where α 2 1 ≡ α 2 ± (ξ) is the double-trace coupling at the fixed point (1.4). Expanding (5.2) in powers of the couplings α 2 1 and ξ 2 we find that the first few terms of the weak-coupling expansion of G 0 are given by graphs depicted in Fig. 16 below. We compute them later in this section. It is easy to see that higher order terms of the weak-coupling expansion of (5.2) produce graphs shown in Fig. 13. Note that for = 0 the conformal symmetry of the correlation function (5.2) is broken. To elucidate the mechanism of restoration of the conformal symmetry of G 0 and the role played by the double traces, we present below the two-loop calculation of the correlation function (5.2).
). This means that the operator H 0 is singular on the space of functions Φ(x 3 , x 4 ) that do not vanish at short distances x 34 → 0. Examining the expression for the eigenfunctions (3.5), we find they scale at short distances as Φ ν,S,x0 (x 3 , for ∆ 1 = ∆ 2 = 1 and ∆ = 2 + 2iν. We recall that computing the correlation functions in the previous section we deformed the integration contour over ν into the lower half-plane. It is easy to see that in this case Φ ν,S,x0 (x 3 , x 4 ) vanishes for x 34 → 0 and, as a consequence, the operator H 0 does not develop UV divergences. Moreover, as follows from the definition (5.1), the double-trace operator V annihilates the eigenstates Φ ν,S,x0 (x 3 , x 4 ) with Im ν < 0 and, therefore, does not contribute. This explains why the double-trace interaction can be neglected when computing the four-point function G 0 by the Bethe-Salpeter method. The appearance of UV divergences at weak coupling is a manifestation of analytic properties of G 0 . As a function of ξ 4 , it has a square-root cut at the origin so that its pertubative expansion runs in powers of (−ξ 4 ) 1/2 . We have already observed this phenomenon on the example of the scaling dimension (4.11).

Double-trace contributions to G B and G C
The inspection of Feynman graphs defining G B and G C shows that the double-trace interactions with the coupling α 1 do not contribute in the planar limit. On the other hand, the interactions with the double-trace coupling α 2 do contribute to both correlation functions through the graphs of the type shown on Fig.15 on the example of G C . Each vertex depicted by white blobes on Fig.15 describes both single-and double-trace couplings. The contribution of each such vertex to G C is UV divergent and it is proportional to (ξ 2 − α 2 2 ). As a result, it vanishes at the fixed point (1.4), so that we are left only with the sums over UV finite single spiral graphs summed up by the UV finite structure function G 1 . This property is not surprising given the fact that the correlation function G C has to be a finite function of ξ 2 whereas the Feynman diagram in Fig.15 involves the ultraviolet divergent scalar loops.
The function G C is regular at ξ 2 → 0 and its weak-coupling expansion runs in powers of ξ 2 . The expression for G C at arbitrary coupling has been derived in section 4.2 using the Bethe-Salpeter equation in the form of conformal partial wave expansion, Eqs. (2.10), (4.39) and (4.40).
The correlation function G B also receives UV divergent contribution from planar graphs similar to those shown in Fig. 15. Their contribution vanishes at the fixed point through the same mechanism    as in the previous case. This means that G B is defined by the two-magnon function G 2 , see Eq. (2.11).
Since two-magnon graphs contributing to G 2 contain even number of single-trace vertices, the weak coupling expansion of G 2 runs in powers of ξ 4 . At arbitrary coupling, G 2 is given by the conformal partial wave expansion (4.60) and (4.61).
In the rest of this section, we compute the first few terms of the weak coupling expansion of the 4-point correlation functions and, then, compare them with the exact expressions obtained in section 4.

Type G D
In the previous sections, G D was defined in terms of the function G 0 by (2.7), and G 0 is given by the wheel diagrams shown in Fig.6. As we already discussed, this definition is not complete and has to be supplimented with the double-trace contributions. The reason for this is that each diagram in Fig.6 is UV divergent and the extra double trace contribution is needed to make each term of the weak coupling expansion to be UV finite.
To illustrate this, we examine the first few terms of the weak coupling expansion of G 0 . They are given by Feynman diagrams shown in Fig. 16 where G (n,k) denotes the contribution of diagram with n double-trace vertices and k single-trace vertices. We recall that G 0 depends on only one double-trace coupling α 2 1 whose value is given by  (1.4) at the fixed point.
The first term on the right-hand side of (5.3) has been previously defined in (4.1), 3) can be computed from Feynman graphs contributing to G D and is given by a finite four-dimensional "cross" integral (see Fig. 16) where the factor 1 2 comes from the relation (2.7) between G D and G 0 , 4 is the symmetry factor, c = 1/(4π 2 ) and D(u, v) has a simple form when expressed in terms of auxiliary variables defined in (3.10) The function inside the brackets is known as the Bloch-Wigner function. The two-loop corrections G (0,2) and G (2,0) come from last two Feynman diagrams shown in Fig. 16. The corresponding Feynman integrals are divergent and require regularization. In dimensional regularization with d = 4 − 2 we have where the notation was introduced for The integral on the right-hand side has a UV divergence coming from integration at short distances x 2 56 → 0. Applying the identity 1/(x 2 56 ) 2−2 = π 2 δ (4−2 ) (x 56 )/ + O( 0 ) we find that the residue of I(x 1 , x 2 |x 3 , x 4 ) at the pole 1/ is proportional to the same one-loop integral that enters (5.4). The same is true for the function I(x 1 , x 3 |x 2 , x 4 ). As a consequence, the divergent part of two loop correction to (5.3) takes the form Here we replaced the double-trace coupling by its value (1.4) at the fixed point α 2 1 = α 2 − . 13 For α 2 1 = α 2 + , the function G 0 (u, v) is given by the same expression with ξ 2 replaced with −ξ 2 . Notice that the O(ξ 2 ) correction to (5.9) comes entirely from the double-trace contribution. In the next section we show that (5.9) is in a perfect agreement with the exact expression for G 0 , Eq. (4.22), which was obtained by resumming the wheel graphs shown in Fig. 6. This is in agreement with our expectations that the double-trace contribution to G 0 can be ignored at finite coupling.

Type G 1
The weak coupling expansion of the one-magnon correlation function G 1 is defined by Feynman diagrams shown in Fig. 8. The contribution of the first two diagrams is given by (4.25). In distinction from the previous case, the corresponding integrals are well-defined in D = 4 dimensions and do not require regularization.
In particular, the one-loop correction G (1) 1 defined in (4.25) can be expressed in terms of the "cross" integral (5.4) Then, the one-magnon correlation function G 1 takes the expected form (4.39) with G 1 (u, v) given at weak coupling by In the next section, we reproduce this expansion from the exact expression (4.40) for G 1 and also produce explicit expressions for higher order terms.

Type G 2
The two-magnon correlation function G 2 is defined by Feynman diagrams shown in Fig. 10. Like in the previous case, the corresponding integrals are well-defined in D = 4 dimensions and do not require regularization.
The weak coupling expansion of G 2 runs in powers of ξ 4 and first two terms are given by (4.42). As was mentioned before, the O(ξ 4 ) contribution to G 2 is given by a two-loop Feynman integral (4.42) that factorizes into the product of one-loop integrals. The latter take the form (5.4) and, as a consequence, it can be expressed in terms of Bloch-Wigner function The resulting expression for the two-magnon correlation function takes the expected form (4.60) with We managed to reproduce both terms of the expansion (5.13) from its expansion over conformal blocks (4.61) by expanding it at weak coupling. For that we followed the procedure explained in the next section 6.

Prediction for the 4-point correlation functions at weak coupling from OPE data
In section 4, we determined the conformal data of the operators that appear in the conformal partial wave expansion (3.16) of the correlation functions at finite coupling. This expansion takes the form of double infinite sums over spins and dimensions. It is not obvious a priori that these sums can be evaluated in a closed form. We demonstrate in this section that, at weak coupling, these sums can be computed order by order in perturbation theory. For zero-and one-magnon functions, G 0 and G 1 , respectively, the result can be written in terms of special class of iterated integrals known as single valued harmonic polylogarithm functions (SVHPL) [60,61]. The two-magnon function G 2 has a more complicated form and it can be expressed in terms of elliptic polylogarithms.
6.1 Zero-magnon case (Type G 0 ) As we already emphasized before, an interesting fact about G 0 is that it receives corrections from the double trace interaction. As a result, its expansion goes in powers of ξ 2 rather than in ξ 4 .
In the previous section we computed the first three terms of the weak-coupling expansion of the function G 0 defined in (4.22). Taking into account (5.9) as well as the explicit expression (5.5) for the function D(u, v), it is natural to look for a general expression for G 0 in the form where z andz are defined in (3.10). The goal of this section is to compute G (n) 0 explicitly starting from the OPE expansion (4.22).
Note that the dependence on the coupling constant ξ 2 enters into (4.22) only through the scaling dimensions ∆ 2 (S) and ∆ 4 (S) given by (4.8). For general S their weak coupling expansion only involves powers of ξ 4 . We recall however, that, due to non-comutativety of the limits ξ → 0 and S → 0, for S = 0 the weak coupling expansion of ∆ 2 (S = 0) does contain powers ξ 2 powers. This means that all terms on the right-hand side of (6.1) with odd powers of ξ 2 come entirely from the contribution of the twist-2 operator with zero spin. We would like to emphasize that the scaling dimension of this operator has to satisfy the condition Re ∆ 2 (0) > 2. Together with (4.10), this implies that ξ 2 should have a nonzero imaginary part Im ξ 2 < 0. For Im ξ 2 > 0 we have to exchange the operator with its shadow whose scaling dimension is given by 4 − ∆ 2 (S) and it can be obtained from (4.10) through the transformation ξ 2 → −ξ 2 . This ambiguity exactly corresponds to the ambiguity of choosing the fixed point (1.4). In what follows we assume that Im ξ 2 < 0 and apply (4.10).
In order to find the explicit expressions for the coefficient functions G (n) 0 (z,z) we match (6.1) into the OPE expansion (4.22) in the short distance limit, u → 0 and v → 1, or equivalently z ,z → 0. In this limit, the conformal blocks scale as g ∆,S (u, v) ∼ (zz) (∆−S)/2 and their expansion in powers of the coupling generates terms of the form z nzm log k (zz) with k not exceeding the order in the coupling. The small z ,z expansion of G (n) 0 (z,z) involves the same terms and their coefficients can be computed for any finite n + m using the exact expressions for the conformal data of the operators.
A nontrivial property of G (n) 0 (z,z) is that for n ≥ 1 they can be expanded over the basis of special iterated integrals, the so-called harmonic polylogarithms (HPL) (see [62,63] where the sum runs over the two sets of indices (including empty sets) a = (a 1 , a 2 , . . . ) and b = (b 1 , b 2 , . . . ) with a i and b j taking the values {0, 1}. Most importantly, the number of terms on the right-hand side of (6.3) is finite for any n thus allowing us to find the expansion coefficients C a, b by matching the small z,z expansion of G (n) 0 (z,z) into the corresponding expansion of the basis of HPL functions. Namely, comparing the coefficients in front of z nzm log k (zz) terms on the both sides of (6.2) with sufficiently large n + m we obtain an overdetermined system of linear equations for C a, b . It is remarkable that for any finite n this system has a unique solution.
For the first few coefficient functions we find where we introduced a short-hand notation for H a1,a2,... = H a1,a2,... (z) andH a1,a2,... = H a1,a2,... (z). The same expressions can be rewritten in terms of classical (di)logarithm functions as We verify that these expressions coincide with the result of the two-loop calculation (5.9).
It is convenient to assign to each term on the right-hand side of (6.2) and (6.3) the weight equal to the total length of the sets a ∪ b. Then, G The coefficient functions G (n) (z,z) should be single-valued functions of z forz = z * . This property is not obvious from (6.2) since HPL functions have, in general, branch cuts that start at z = 0 and z = 1. We can use this information to restrict the basis of possible functions. The HPLs can enter (6.3) through special linear combinations which are free from the branch cuts. They are known as single-valued harmonic polylogarithms L a1,a2,... , their definition can be found in [60,61]. The resulting expressions for G Similar expressions up to 7 loops can be found in a Mathematica file attached to this submission.
6.2 One-magnon case (Type G 1 ) Going along the same lines as in the previous case, we were able to expand the sum (4.40) in terms of the SVHPLs. The weak-coupling result (5.11) suggests to look for G 1 in the form The explicit expressions for the first few coefficient functions are We verify that the first two terms reproduce correctly the perturbation theory result (5.11).

Two-magnon case (Type G 2 )
In this case we only managed to reproduce the tree level and 2-loop perturbation theory result (5.13). We found that it is not possilbe to express the 4-loop expression in terms of SVHPLs. We found that it is given by elliptic function. In small z limit the underlying elliptic curve degenerates and the correlator can be written in terms of HPL's ofz multiplied by powers and logarithms of z. 15

Classical (Strong Coupling) Limit of the 4-point Correlators
In this section, we investigate the strong coupling limit of the 4-point correlators. Even though the world-sheet description of this theory is still not known, it was shown in [44] that there is a classical limit of the underlying integrability construction for the spectrum where it reduces to an algebraic curve reminiscent of that of the classical strings in AdS 5 × S 5 in the full theory. Similarly, we will see that the leading strong coupling ξ → ∞ asymptotics of the correlation function is saturated by one state with large ∆, S ∼ ξ and the corresponding result scales as e −ξA(z,z) , where A(z,z) is a certain function of cross-ratios. Moreover, since ξ = ge −iγ3/2 this classical asymptotics reminds the behavior of the three point correlation functions for short operators in strongly coupled N = 4 SYM theory in planar limit (see for example [64,65]). Whereas it would be premature to conclude that this points to the dual string description for the bi-scalar model, the observed behaviour looks very much like an action evaluated on some classical solution. All that constitutes an evidence towards existence of the classical limit of the fishnet theory at strong coupling ξ → ∞. The possible relation to the recently proposed AdS sigma model description of the ground state of the fishnet theory for infinite L by [66] is also yet to be understood. In this section we assume that z takes a generic value in the range 0 < |z|< 1 16 . The convergence is not uniform and the limits z → 0 or z → 1 have to be taken with extra care. We also assume that ξ = e −iφ ξ 0 where ξ 0 is large and real and φ is a small positive phase. This assumption is necessary as the correlator has poles at real values of ξ which accumulate at infinity and the limit is not defined.

Correlation function for zero-magnon case (wheel-graphs)
First we consider the correlation function G D . It is obtained from the zero-magnon function G 0 given by (4.22) via (2.7). As we discussed above, this results in dropping terms with odd S in the sum (4.22) and doubling the terms with even S. In this subsection, we compute (4.22) in the limit when ξ → ∞.
Our main assumption, which is backed by intensive numerical analysis, is that for ξ → ∞ the sum in (4.22) is saturated by large spins S ∼ ξ. Then, replacing the conformal blocks in (4.22) by their asymptotic behavior at large S and ∆, we can evaluate the sum over S by the saddle-point method.
In what follows we only evaluate the leading exponential factor. The pre-exponent and the subleading terms can be computed by the same method, we leave it to future studies.
Before we begin we notice the following property of the structure constant (4.15) It allows us to write the 4-point function (4.22) in the following way where k(∆, z) is given by the hypergeometric function in (A.2). This expression considerably simplifies our analysis as it allows to replace the sum over S by an integral with exponential precision at large ξ and the evaluate it by the saddle point method.
The asymptotic behavior of the conformal block can be found by using a series representation of the hypergeometric functions in (A.2) Then rescaling the variables as ∆ = ξ d , S = ξ s , k = ξ w , (7.4) 16 As below we find square roots and logariphms for definiteness we also assume in the calculation that 0 < arg z < π/4. we expand the expression under the sum at large ξ and extremize in w to obtain (7.5) in agreement with [67]. Similarly we replace the structure constants (4.15) by their leading asymptotic behaviour at large S and ∆ to get from (7.2) and (2.7) (7.6) where the sum runs over even spins S, both positive and negative. We recall that for each S there are only two values of ∆ that contribute to the sum, ∆ 2 (S) and ∆ 4 (S), given by (4.8). At strong coupling, we apply (7.4) to find that they have different dependence on the spin Note that for s < 2 we get purely imaginary ∆ t=2 /ξ. In order to get a consistent strong coupling limit we should take ξ to have slightly negative phase, so that Re ∆ t=2 → +∞ at strong coupling. After that we substitute (7.4) into (7.6) and extremize the expression in the exponent over s to get where the two states (7.7) produce two different exponents wherez = z * . Notice that A t=4 is positive whereas A t=2 is purely imaginary. As a consequence, −ξA t=2 has a large negative real part and, therefore, the correlation function (7.9) is exponentially suppressed at strong coupling, as it is expected for the tunneling processes in the classical limit.
Let us also point out that the state which saturates the sum over S in the sum (7.6) has the following scaling dimension and spin with the upper sign for t = 2 and lower for t = 4. It would be interesting to find a classical model which reproduces these results. 17 As a test of our result we consider the small z limit of (7.8). In this limit, we expect that the leading contribution to the correlation function should come from the states with S = 0 leading to G D (z,z) ∼ (zz) ∆t(S=0)/2 . Expanding the relations (7.9) and (7.10) at small z we find The leading term on the right-hand side agrees with analytic result (7.9) for Re(−ξA t=2 (z,z)) and Re(−ξA t=4 (z,z)) up to 27 significant digits (i.e. within the fit precision 18 ). The relation (7.12) clearly indicates that the first subleading correction to log G D is the same for the two states. It scales as ln ξ/2 and generates √ ξ in pre-exponents on the right-hand side of (7.8). This factor could come from some zero modes in the semiclassical analysis (see very similar prefactors in the expectation values of a circular Wilson loop in N = 4 SYM [68]).

Correlation function for one-magnon case (single spiral graphs)
The strong coupling limit of the correlation function G 1 happens to be very similar to the previous case. We immediately observe using (4.31) that the states contributing to G 1 have exactly the same strong coupling asymptotics as t = 2 states in G 0 from the previous section, namely, for even S we have ∆ ξ s 2 − 4 (7.13) where s = S/ξ and the odd S have to be considered separately and their contribution can be obtained by replacing ξ by −iξ. From (7.13), it is not surprising that the leading asymptotics of the 4−point correlator appears to be also the same as that of G D , Eq. (7.8) G 1 (z,z) ∼ e −ξAt=2(z,z) + e −ξAt=4(z,z) . (7.14) where A t=2 comes from even spins and A t=4 from odd spins. Note, however, that the pre-exponential factors in (7.8) and (7.14) are different.
To verify our result we computed G 1 (z,z) numerically for fixed values of ξ = ne −iπ/6 , n = 100, 105, 110 . . . , 200 and for fixed z = 1 8 e iπ/5 . Fitting this data with n, log n, 1, 1/n, 1/n 2 , . . . , 1/n 18 we obtain the following result The leading coefficients agrees with our analytic result −Re(e −iπ/6 A t=2 (z,z)) and −Re(e −iπ/6 A t=4 (z,z)) with 27 digits (i.e. within the fit precision). * * * The analysis of the strong coupling behavior of the two-magnon correlator (double spiral graphs) is more complicated and we leave it for the future studies. It would be interesting to guess the "worldsheet" degrees of freedom leading to these classical asymptotics, similarly to what was done AdS 5 × S 5 in the full AdS 5 × S 5 duality. A possible way to further elucidate whether the bi-scalar model has 18 For presentation purposes we lowered the precision in (7.12)

Gluing triangles: general structure constants
In the previous section, we computed the OPE coefficients C ••• ∆,S for the two protected and one unprotected zero-magnon operator, see (4.16). In this section we generalize these results to the OPE coefficients of one protected and two unprotected zero-magnon operators C ••• ∆1,S1,∆2,S2 . We show below that these OPE coefficients can be easily obtained by gluing together the threepoint correlation functions defined by where n is an auxiliary light-like vector and with dots denoting terms with derivatives acting on both fields. As usual, the normalization factor N is fixed by the two-point function, The stucture constant C ••• ∆,S is given by (4.15) and (4.16) and the scaling dimension ∆ = ∆ 2 (S) is defined in (4.8).
The three-point function (8.1) resums the wheel diagrams shown in Fig.17. They contain two scalar lines connecting x 1 with the two external points x 2 and x 3 dressed by an abritrary number of wheels incircling x 1 . We can use this result as a building block for a more complicated three-point correlation function, involving two twist-2 operators with an arbitrary spin and one protected operator whereŌ n2 S2 is given by (8.2) with X replaced withX and n i (with i = 1, 2) being light-like vectors. In the planar limit, this correlation function receives contribution from the diagrams shown on Fig. 18. The main observation is that they can obtained by merging together two wheels depicted in Fig. 17.
To exemplify the idea we first consider (8.3) for S 1 = S 2 = 0. In this case, we have where the operator (− 0 ) amputates an extra propagator and the integration over x 0 glues two sets of wheels together. Here ∆ 1 = ∆ 2 are the scaling dimensions of the operators tr(X 2 ) and tr(X 2 ) given by (4.10) but it is convenient to keep ∆ 1 and ∆ 2 to be arbitrary. 20 Then, after differentiation in x 0 the integral becomes We can employ inversions to verify that the integral transforms under the conformal transformation as a three-point correlation function of scalar operators. This property fixes the form of T 0,0 (x 1 , x 2 , x 3 ) up to a structure constant The integral in (8.5) can be computed immediately using the star-triangle identity leading to Replacing C ••• ∆,S with (4.16) we obtain the following expression which is symmetric in ∆ 1 and ∆ 2 as it should be. So far, we treated ∆ 1 and ∆ 2 as arbitrary parameters. We expect that, in the limit ∆ 2 → 2, when the correlator T 0,0 looses all its wheels around x 2 , it should reduce to the wheel correlation function (8.1) multiplied by a free scalar propagator connecting x 2 and x 3 . Indeed, it is easy to see that C ••• ∆1,0;2,0 = C ••• ∆1,0 . In application to the conformal fishnet theory, we have to replace the scaling dimensions ∆ 1 = ∆ 2 in (8.8) with their exact expression (4.10). In this case, the expression (8.8) simplifies to Non-zero spins. It is straightforward to generalize (8.4) to the operators with non-zero spins where the notation was introduced for the so-called conformal triangle function
Applying (8.19) we have to replace the scaling dimensions ∆ 1 and ∆ 2 by their explicit expressions. For the twist-two operators (8.2) they are given by the function ∆ 2 (S), Eq. (4.8), evaluated for spin S = S 1 and S = S 2 , respectively. We recall that the scaling dimensions of twist-2 and twist-4 operators, ∆ 2 (S) and ∆ 4 (S) are two branches of the same function of the (complexified) coupling constant ξ 2 . This suggests that depending on the choice of the branch of the functions ∆ 1 (ξ 2 ) and ∆ 2 (ξ 2 ), the relation (8.17) should also describe the three-point correlation function of the protected operator with two unprotected operators each having twist−2 or twist−4.

Generalization to any dimension
Many of the results obtained in the previous sections for the bi-scalar theory (1.1) in d = 4 dimensions can be easily generalized to its d−dimensional version proposed in [47]. In particular, the zero-magnon four-point correlation function G D was computed there for any d, generalizing the d = 4 results of [35].
The Lagrangian of d−dimensional bi-scalar model is non-local where the differential operator in an arbitrary power is defined in a standard way, as an integral operator. For the particular "isotropic" case ω = 0, the action (9.1) should be supplemented with the same double-trace counterterms (1.2). As before, the theory has two fixed points with the corresponding values of the double-trace couplings depending on ξ 2 and computable at least at weak coupling. Since the d-dimensional theory (9.1) has the same chiral interaction vertex as (1.1), we can consider the correlation functions the same types as those shown on Fig.1. An important difference with the d = 4 case is that free scalar propagators are now given by c/(x 2 12 ) d/4 with c = 1/(2π) d/2 . To compute zero-and one-magnon four-point correlation functions, as it is done in 4 dimensions in the previous sections, we have to find, at any d, the eigenvalues of the graph-building operators E 0 and E 1 . They are computed in appendix C leading to (see also [47,75]) In addition, we also need the expression for the kinematical factor c 2 defined in (A.6). It is given in by (A.6), where we should take ∆ 1 = ∆ 2 = d/4 for the zero-magnon case and ∆ 1 = 2d/4, ∆ 2 = d/4 for the one-magnon case. Assuming that all the intermediate steps are still valid for general d, we find the the following expression for the OPE coefficients where χ 0 = (4π) d ξ 4 , n = 1 and χ 1 = (4π) d/2 ξ 2 , n = 2 for zero-magnon and one-magnon cases, respectively. The dimension ∆ of operator appearing at the corresponding pole is related to the representation label ν in the following way: ∆ = d 2 + 2iν. For d = 2, the zero-magnon spectrum looks particularly simple [47] ∆ 0 = 1 + S 2 − 4ξ 2 . (9.5) Interestingly [47], the particular case of d = 2 with ω → 1/2 in (9.1) is relevant for the BFKL approximation in high-energy QCD.
For general d the structure of the poles in ν could change. In particular the poles coming from E ∆,S = 1/χ will have a different number of solutions for different d (actually, this number is infinite for odd dimensions, see [47] for the analysis of possible exchange states in zero-magnon case), and those have to be taken into account when computing the correlation function. The computations are particularly simple in dimensions d = 4k, where k is integer. Then the spectral equations become polynomial in ∆. For example, for one-magnon case, the spectral equation for d = 8 looks as etc.
We leave the computations of a more complicated two-magnon correlation function, as well as the detailed study of the general d case, for future investigation.

Discussion and Conclusions
In this work, we computed exactly and explicitly certain 4-point correlation functions in bi-scalar CFT (1.1) which represent a specific double scaling limit (combining weak coupling and strong imaginary twist) of γ-deformed N = 4 SYM theory in 't Hooft approximation [42]. Although this theory obeys the integrability properties, related to the integrability of fishnet Feynman graphs [48] dominating in its planar perturbation theory, we concentrate here on the physical quantities which can be computed without any appeal to the integrability and utilizing only the conformal symmetry properties. The three types of such 4-point functions are named as zero-magnon, one-magnon and two-magnon cases, referring to the exchange operators which are of a type tr(Z 2 ), tr(Z 2 X) and tr(Z 2 X 2 ) + perm., as well as their non-zero spin cousins. The last operator is the analogue of the famous Konishi operator in N = 4 SYM theory. It is quite remarkable that we managed to compute explicitly the all-loop anomalous dimensions of all these operators, as well as their structure constants with the external protected operators. These formulas give a host of very non-trivial non-perturbative OPE data, with a rich analytic structure. To our knowledge, this are the first examples of such non-perturbatively and explicitly computed 4-point correlation functions in an interacting CFT in d > 2 dimensions.
The results for 4-point functions are presented in a standard form, as explicit OPE expansions over conformal blocks. We studied these functions by the weak coupling expansion and found that the results for zero-and one-magnon cases can be presented, in each order of weak-coupling expansion, in terms of special class of functions called Harmonic Poly Logariphms (HPL), thus facilitating the study of these functions in the cross channels. As for the more complicated two-magnon case, we observed that the HPL representation is only possible in the light-cone limit z → 0, or equivalently u → 0 with v fixed.
Many interesting questions related to the study of these quantities are left out of the scope of this paper. In particular, we left for the future the study of the operator content and the OPE in cross-channels of all three 4-point functions, where the analytic structure looks more complicated than in the original channel. It would be also interesting to generalize these 4-point functions to the case of external operators with spins. The correlation functions we computed could serve as a building block for more complicated n−point functions by gluining them together. Similar approach was recently developed in SYK context in [76] 22 . To demonstrate the procedure we show in section 8 how to obtain more complicated structure constants starting from the elementary blocks.
Notice that the main common feature of the 4-point functions, which allowed for the explicit computations using the combination of Bethe-Salpeter techniques and the conformal symmetry, was the presence of not more than two scalar fields of each species, Z or X, in the exchange operators. If we want to compute the multi-point correlation functions with more than two scalars of any of the species in exchange operators, we have to appeal to the integrability methods, based on the non-compact conformal Heisenberg type spin chains. All computations become then much more sophisticated, though in many cases possible. Already the computations of anomalous dimensions of tr(Z 3 ) operator and its cousins with the same R-charge, necessitated the extraction of the results from the quantum spectral curve of the full γ-deformed N = 4 SYM theory, in the corresponding double scaling limit [44]. The alternative method for computations of anomalous dimensions, is based on the integrable quantum conformal spin chain [78]. For the 3-point and 4-point correlation functions in bi-scalar theory, this last method seems to be particularly promising. Recently, in a very similar setup it was shown that the separation of variables (SOV) approach could be used to get compact expessions for the correlation function in [22]. The expressions for the structure constants obtained there are very similar to our results, which indicates that a similar SOV-based approach could work here. There are also some other options, such as the form factor approach [21], successfully applied in [52] for the computation of certain 4-point correlation function based on fishnet graphs with disc topology.
Interestingly, all our current results for 4-point correlation functions in d = 4 dimensions appear to be directly generalizable to the bi-scalar CFT at any d which was recently formulated in [47]. The d = 2 case is closely related to the BFKL model of reggeized gluons describable by the SL(2, C) Heisenberg spin chain [50,51,79,80]. It is worth noticing that the computations of 4-and 6-point correlation functions of BFKL pomeron light-ray operators presented in [81][82][83] look closely related to our current computations. It would be interesting to compare those results to our current results.
The most interesting, physical questions about the bi-scalar CFT still remain to be answered. It still remains to be understood whether this theory has a dual string description. The availability of explicit all-loop results, such as presented in this paper, should allow for some guesses in this direction. Quite intriguingly, the strong coupling limit of the 4 point correlators studied in this paper exhibits typical classical exponential scaling with the coupling, suggesting existence of the dual strong coupling classical description.
Note added: while this paper was in preparation a paper [84] appeared about SYK model, which may however have some overlap with our results.
The sum of the above contributions in general is not zero, however, there are also poles at ∆ = S + 3 + k, k = 0, 1, 2, . . . coming from 1/c 2 (∆, S) factor in (3.12). They can be also expressed in terms of r n as follows 1 c 2 (S + 3 + k + 2i , S) − 1 r k c 2 (S + 3, S + k) we have to consider We see that there will be no additional terms in (3.16) if we require r k (E 3+S+k,S − E 3+S,S+k ) = 0 , k = 0, 1, 2, . . . .

(B.8)
We verify in the main text that this requirement is indeed satisfied for each case considered.

C Eigenvalues of the graphs-building operators
In this appendix we give technical details of the derivation of the eigenvalues of the graphs-building operators H 0 , H 1 and H 2 .

C.1 Star-triangle identity
The calculation can be simplified by applying the so-called star-triangle identity where the exponents a, b and c satisfy a + b + c = −2d. A particular case of it is which can be obtained from (C.1) by sendind point x 3 to infinity. Again we can absorb the last term under the integral into the derivatives in x4 and evaluate the integral using (C.2) (C.13) Finally, computing the derivatives and comparing the result with (C.11) we arrive at (C.14) which reduces at d = 4 to (4.29).

D Analytic properties of the correlation functions
We can apply (3.11) and (3.12) to understand analytic properties of G(x 1 , x 2 |x 3 , x 4 ) as a function of the coupling constant ξ 2 . The integrand of (3.12) is a meromorphic function of χ which is just a power of ξ 2 dependent on the type of the correlation function. The integral in (3.12) is well-defined as soon as the physical poles (3.14) are away from the real axis. As soon as the physical pole (3.14) approaches the real axis, the integral (3.12) generates a branch cut in χ. The discontinuity across the cut can be found from (3.12) In virtue of (3.13), the integral localizes at two points ∆ = 2 + 2iν and 4 − ∆ = 2 − 2iν satisfying (3.14). Taking into account (3.17) we obtain According to (4.5) and (4.47), E 0 and E 2 are positive definite functions of ν and, therefore, the integral on the right-hand side of (D.1) vanishes for χ < 0 for zero-and two-magnon functions, G 0 and G 2 , respectively. Moreover, the function E 2 satisfies the relation 0 < E 2 ≤ 3ζ(3)/(128π 4 ) on the real ν−axis and, as a consequence, disc χ G 2 (u, v) = 0 for χ > 128π 4 /(3ζ (3)).

(E.4)
where A(g) and B(g) are real functions of the couplping and γ j . which is defined only through the universal quantity -the anomalous dimension γ(ξ). This type of RG flow for the double-trace coupling also first appeared in the papers [87,88]. The physical quantities of the theory have a cut along the half-axis ξ > 0 and poles. For generic complex ξ's there is no such problem. Assuming the coupling to have a little negative imaginary ξ → ξ − i the RG flow (E.12) will avoid the poles in the r.h.s. and interpolate between two fixed point α 1,+ or α 1,− (one of them is IR fixed point, another is the UV fixed point, depending on the sign of ). We illustrate the possible RG flows in Fig.19.