Separation of Variables in AdS/CFT: Functional Approach for the Fishnet CFT

The major simplification in a number of quantum integrable systems is the existence of special coordinates in which the eigenstates take a factorised form. Despite many years of studies, the basis realising the separation of variables (SoV) remains unknown in N=4 SYM and similar models, even though it is widely believed they are integrable. In this paper we initiate the SoV approach for observables with nontrivial coupling dependence in a close cousin of N=4 SYM - the fishnet 4D CFT. We develop the functional SoV formalism in this theory, which allows us to compute non-perturbatively some nontrivial observables in a form suitable for numerical evaluation. We present some applications of these methods. In particular, we discuss the possible SoV structure of the one-point correlators in presence of a defect, and write down a SoV-type expression for diagonal OPE coefficients involving an arbitrary state and the Lagrangian density operator. We believe that many of the findings of this paper can be applied in the N=4 SYM case, as we speculate in the last part of the article.


Introduction
The Separation of Variables (SoV), which traces its origin to the Hamilton-Jacobi approach and later to Sklyanin's works [1][2][3][4], is often regarded as the most powerful method to solve quantum integrable systems (for recent developments see [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]). It is based on the expected property that eigenfunctions of the integrals of motion factorise when evaluated in a special basis x| labelled by the "separated variables" x ≡ {x n } L n=1 , where L is the number of degrees of freedom. Schematically, the eigenstates decompose as where the factors can be determined by solving linear functional equations in one variable. Typically, they coincide with the so-called Q-functions (or simple combinations of them), and the functional equations are known as Baxter TQ equations or, in certain contexts, Quantum Spectral Curve (QSC) equations, which also include nontrivial analyticity conditions. When it can be worked out, the SoV gives access not only to the spectrum and eigenfunctions of the integrals of motion (IM), but also usually leads to simple expressions for correlation functions (see e.g. [21,22] for examples in two-dimensional integrable field theories and recent work [7] for spin chains). The SoV paradigm is expected to hold true for quantum integrable systems carrying any representation of the symmetry algebra, in contrast e.g. with the Bethe Ansatz (BA), which can be applied only in special cases. In particular, there is a growing body of evidence that this approach is applicable for such complicated systems as 4D N = 4 SYM. In [23][24][25][26], it was shown that some correlation functions at finite coupling, obtained by various techniques such as supersymmetric localisation or the direct resummation of diagrams, can be written as an integral of a product of Q-functions, in agreement with the type of structures expected from SoV. This is especially exciting as these examples appear in the non-perturbative regime of observables involving short operators, where the BA methods are not applicable. 1 This gives hopes that the complete non-perturbative solution of planar N = 4 SYM could be obtainable by means of SoV methods. The result for correlators should be given in terms of the Q-functions, which are, luckily, already under full control at finite coupling: they are solutions of the QSC equations [30,31] 2 , obtained in one-toone correspondence with the anomalous dimensions of primary operators. In addition, there are other examples of the applicability of SoV for studying various corners of the theory [9,12,[37][38][39][40][41], and these techniques could eventually lead to the first principles derivation of the BA-based Hexagon formalism [27,28].
Developing the SoV approach in N = 4 SYM from first principles is very challenging. The main conceptual difficulty is that we do not even have a precise definition (in terms of 1 At weak coupling or for very long operators, one can get spectacular results with BA-inspired methods [27,28], and even go beyond the planar limit [29]. 2 The QSC formulation is also known for the ABJM theory [32][33][34] and for some non-local operators [35].
A similar construction with nontrivial analyticity conditions for the Q-functions also exists for the Hubbard model [36].
a Hilbert space and Hamiltonian) of the integrable system which controls the spectrum of N = 4 SYM at finite coupling. It can be in principle obtained in AdS/CFT by quantising appropriately the worldsheet σ-model, which has not been done beyond a few quasi-classical orders. The situation improves considerably if we consider a theory that is closely related to N = 4 SYM, the so-called fishnet CFT. This model was obtained as a double-scaling limit of the γ-deformation of N =4 SYM in [42]. Remarkably, it is a non unitary, non supersymmetric, but exactly conformal theory (at least in the planar limit) defined by a very simple Lagrangian: 3 where the complex scalar fields φ i , i = 1, 2 are N c × N c matrices. The model inherits integrability from the γ-deformed N = 4 SYM: in particular, the Q-functions and QSC equations can be understood as a limit of the ones in the "parent theory" [44]. However, in the fishnet model this integrable structure can be understood much more clearly [42,[44][45][46], making it an ideal playground for developing the SoV program.
Integrability arises directly as a property of the fishnet Feynman diagrams, which were already studied by Zamolodchikov in [45]. The resummation of fishnet graphs with the topology of a cylinder, and the related Dyson-Schwinger equations, define rigorously the integrable system associated with the spectral problem. It consists in a spin chain carrying an infinite-dimensional representation of the 4D conformal group. The role of Hamiltonian is played by the graph-building operator which constructs the Feynman diagrams [44], or more conveniently its inverse, which acts as a differential operator. In the following, we will refer to the conformal spin chain with this Hamiltonian as the fishchain, following [47].
As we review in the main text of the article, the eigenstates of the integrable system can be interpreted as certain trace-trace correlators introduced in [48], and called CFT wave functions. The CFT wave functions encode all information on the single-trace operators at planar level, and moreover appear as natural building blocks in correlators, see Figure 1. One expects that, written in a basis of Separated Variables, they can be computed in terms of the Q-functions solving the QSC equations. The chain of relations between operators, wave functions, and Q-functions, is illustrated here: For the simplest family of nontrivial operators, those with length one in the presence of twists [49], the map between wave functions and Q-functions was found by the present authors with A. Sever in [50]. Generalising this result to all states is currently an important open problem, but the general methods of [4,7,20] give a clear indication of how to construct the SoV transformation at least formally. The explicit construction of such map is referred to as the operatorial approach to SoV. In this paper, we do not address this problem yet. However, we take a functional approach [15], which is based only on the Baxter TQ equations, and allows us to derive a number of mathematical results related to the structure of the SoV.
A key element of our construction is the introduction of inhomogeneities, and arbitrary scaling weights, on the sites of the spin chain. Such parameters are typically very useful in the SoV program. Here we introduce them for the first time for the single-trace operators in the fishnet model at finite coupling. Following [7,15], we show how variations of these parameters can be used to generate nontrivial operators acting on the spin chain, and to compute their expectation values in terms of Q-functions, as we discuss below.
Another important technical ingredient used in this paper are quasi-periodic, or twisted, boundary conditions along the spin chain, which break the conformal symmetry of local operators [49]. The twists are crucial for the SoV, since they lift the degeneracies of the spectrum and guarantee that solutions of the QSC are identified one-to-one with the eigenstates of the spin chain. We use the construction of the colour-twist [49], which is particularly simple for the fishnet model, and introduces the twists directly at the level of Feynman graphs.

Summary of the main results
In the first half of the paper we develop the functional SoV formalism for the fishchain. In the second part, we discuss some applications, which broadly concern three topics: scalar products, spin chain form factors, and the g-function. We also speculate on the applicability of our methods to N = 4 SYM.
Integral orthogonality relations. It was shown in [15] that one can realise the orthogonality properties of different integrable eigenstates, in terms of an integral over certain combinations of their Q-functions. This type of orthogonality relations in the context of the spin chains are very useful, e.g. made it possible [7,13] to extract the SoV measure explicitly for higher-rank models.
In the previously studied cases, the Q-functions were polynomials, but we demonstrate in this paper how the approach extends to the situations when all Q-functions are nonpolynomial. The key building blocks for the integral orthogonality relations are integrals of the product of two Q-functions, which we call Q-bilinear forms. We discuss in detail their properties. We show that some of these forms are finite precisely for "on-shell" solutions of the TQ-relations, i.e. those satisfying the quantisation condition giving a discrete spectrum of conformal dimensions. This gives a new interpretation of the quantisation conditions, which were previously quite mysterious and could be justified mainly by invoking properties of the QSC in N = 4 SYM.
Variation of parameters and spin chain form factors. We obtained an analytic formula for the variation of any integral of motion of the fishchain with respect to a parameter. The expression, similar to those found in [7,13] for HW spin chains, consists in a ratio of determinants built with integrals over products of Q-functions. Such observables can be interpreted as spin chain form factors, i.e. expectation values of operators acting on the spin chain. Since we have ∼ 4J integrals of motion, and we can use as a parameter any of the 2J inhomogeneities + weights, a naive counting suggests that we have access to ∼ (4J) × (2J) different observables, where J is the length of the chain.
The simplest example is the derivative of the scaling dimension ∆ A of an arbitrary single-trace operator O A with respect to the coupling constant. It is well known from conformal perturbation theory that such variation can be identified with the OPE coefficient [51] ∂ where L int = Tr(φ 1 φ 2 φ † 1 φ † 2 ) is the interaction vertex in the Lagrangian. Even though such correlators can also be extracted from the known spectrum, our result gives a closed answer in terms of the Q-functions at a given value of the coupling ξ -as a determinant of one-dimensional integrals of Q-functions. Based on the example of [23], we expect this to provide rich structural information for generalisations to more complicated OPE coefficients. Correlators of the type (1.4) were also recently studied in [52] with the hexagon decomposition approach. Understanding fully the relation with our result could help to advance with the resummation of the Hexagon series. We were also able to compute the expectation value of the variation of the Hamiltonian w.r.t. local weights, which is a natural observable and provides an exact result for a class of Feynman diagrams with nontrivial local insertions.
The structure of the g-function. We will also make some observations on the computation of the so-called g-functions inspired by the recent papers [53] and [6]. The g-functions capture the overlap between a spin chain eigenstate and an integrable boundary state (for a more detailed definition see section 7), and have appeared in connection to several observables in N =4 SYM [54][55][56][57]. The traditional integrable formalism for computing the g-function [58] is based on the thermodynamic Bethe Ansatz, and organises the result as a product of two factors. One, simpler to compute, is boundary-dependent, while the other is a universal factor which depends only on the state. In [53], a conjecture was presented expressing the universal factor in terms of determinants built with the Q-functions, for the case of the Sinh-Gordon model. This proposal -which still contains an undetermined proportionality constant -passes an infinite number of nontrivial consistency checks, as it reproduces a selection rule for the g-functions. Moreover it matches the SoV structure of the overlaps which was later established rigorously for the Heisenberg spin chain in [6]. The arguments of [53] are in the spirit of the functional SoV method, and are based on the analysis of the Q-function representation for the scalar product. We show that a similar analysis can be extended to our case, and present the analogous conjecture for the fishnet model.
Finally, we point out that there are many aspects of the functional SoV setup that can be generalised directly to N =4 SYM. In particular, it is quite clear how a notion of scalar product as a multiple integral of the Q-functions can be obtained in this theory, as we discuss in section 8. In turn, this gives a good starting point to construct the g-function, which is expected to correspond to the TBA expressions written down in [54,55].
The key peculiarity of the N = 4 SYM case is that, even for short operators, the underlying integrable system has an infinite number of degrees of freedom at finite coupling (which is also expected from the dual worldsheet description). So the closest analogy is [21,22], where indeed the size of the determinant was infinite in the SoV-type integrals appearing in the description of correlators. Nevertheless, at weak and at finite coupling, it should be possible to truncate the determinants to a finite size, as the number of relevant parameters in the Q-functions can be taken to be finite with very high precision when studying the spectrum.
The paper is organised as follows: In section 2, we review the general spin chain formulation, presenting a precise definition of the scalar product, and we summarise the Baxter equations and properties of the Q-functions, introducing the new ingredients of the inhomogeneities and weights. In section 3 we link the general spin chain formalism with the fishnet CFT. In section 4 we describe technical details about an important notion of Qbilinear forms, and discuss the new interpretation of the quantisation condition. In section 5, we derive the form of the variation of the integrals of motion. In section 6, we discuss the orthogonality relations and the expectations for the scalar product in the SoV basis. In section 7, we discuss the natural candidate for the universal factor of the g-function in terms of our construction. In the last section 8, we sketch how the arguments of this paper can be generalised to N = 4 SYM. We close with a short summary, and some appendices where technical details are unrolled.

Conformal spin chain
In this section we review in general the spin chains with 4D conformal symmetry at each site. We start by discussing the conformal spin chain abstractly. We generalise previous constructions, e.g. [46], by considering arbitrary conformal weights on different sites and arbitrary inhomogeneities (similarly to [59]). We construct a natural notion of conformally invariant scalar product on the spin chain, and we review the Baxter TQ equation, which is the engine behind the integrability construction.
At the end of the section we also review the quantisation condition for the wavefunctions and its counterpart in terms of the Q-functions, which singles out a discrete spectrum of anomalous dimensions as a function of one parameter, which will be identified with the 't Hooft coupling in the next section.

Representation of the conformal group
In this section we introduce some notations for the representations of the Euclidean conformal group SO(1, 5) which we use below. We found it beneficial to use the 6D formalism especially when it comes to fusion of transfer matrices [46]. Some expressions in 4D can be found in e.g. [60].
For simplicity we consider only scalar representations on each site, with a generic scaling dimension h, realised as functions f (x µ ) of points x ∈ R 4 . The space of such fields will be denoted as F h . We will treat F h as a real vector space, so that we do not have to make any assumption on the reality of h. We will often focus on the h = 1 or h = 2 cases, which are connected to the fishnet theory, and correspond (in the Euclidean signature in which we work) to a representation of the conformal group without a highest nor lowest weight vector.
In 6D space, the conformal group generators take the standard SO(1, 5) form which satisfy the standard commutation relation Here η is the 6D Minkowsky metric with signature −+++++. The action on the functions of 4 variables is defined by acting on the expression As X M X M is invariant under SO(1, 5) we can set it to zero and thus exclude X − ≡ X −1 − X 0 . As furthermore the action byq M N commutes with the re-scaling X M → αX M and its eigenvalue is measured by h, we see that the form (2.3) must be preserved under the SO(1, 5) action. In this way we induce the action of (2.1) on functions f (x µ ) of 4 variables. The generatorsq M N can be related to the dilatation D, translation P µ and special conformal transformations K µ as followŝ On can deduce from (2.1) and (2.3) the explicit action on the 4D fields It is also useful to note that for a finite conformal transformation g ∈ SO(1, 5), the representation acts as where g(x) denotes the action of the conformal transformation on a point x ∈ R 4 . The Hilbert space of the spin chain is a space of functions of J variables transforming in the tensor product representation where for generality we consider different weights at every site and we denoted {h} ≡ {h 1 , . . . , h J }. On this space we can define an action of the conformal group at each site separately, and denote it by g α The global action we denote as g is given by the product J α=1 g α as usual. We will also use the notation for the sum of all weights, (2.11)

The conformally invariant scalar product
Now we introduce a scalar product, invariant under the conformal group: The requirement that fixes completely the form of the scalar product up to an arbitrary overall constant factor, see [61] and Appendix B. The result is given by: where D = 4 and in general we define the fractional power of the d'Alambert operator by (see e.g. [62]): For the Hilbert space of the spin chain the scalar product is defined by We will see in section 2.4 that the above scalar product needs regularisation in some important cases.

Integrable fishchain with inhomogeneities
Here we introduce notations for the spin chain and build the integrals of motion via transfer matrices. Following [46,48], we define an integrable system on the spin chain introduced above. At this stage, we will keep the weights {h 1 , . . . , h J } completely generic, even though, in applications to the fishnet CFT, the relevant weights are either h α = 1 (non-magnon) or h α = 2 (magnon).
To define the integrable system, we introduce a family of nontrivially commuting transfer matrices following the construction of [44-46, 48, 60]. Their eigenvalues provide us with a complete family of integrals of motion. We will also introduce inhomogeneities in the definition of the transfer matrices, which will be labelled by {ϑ α }, 1 ≤ α ≤ J, and twisted boundary conditions.

Commuting transfer matrices
The transfer matrices will be denoted asT r with r ∈ {1, 4, 6,4,1} labelling the representation in the auxiliary space, corresponding to antisymmetric tensors with 0, 1, 2, 3, 4 indices. To build them, we start from Lax matrices, which are differential operators acting on a single site of the spin chain, and evaluated as matrices in the auxiliary space. The Lax matrices defined in [46,48] are: 4 ). An explicit realisation of such matrices, which we will assume in the following, can be found in Appendix A. The transfer matrices are defined aŝ where G r is an element of the conformal group in the representation r, introducing quasiperiodic boundary conditions, and ϑ i 's are the inhomogeneities 5 .
Since the Lax operators satisfy the Yang-Baxter equation as shown in [46,48], the transfer matrices commute for different values of the spectral parameter and for any choice of the auxiliary space representation: T r (u),T r (u ) = 0 .
(2.23) Therefore, they can be diagonalised simultaneously. In the following, we often restrict the discussion to the eigenvalues which are denoted as T r (u).

Parametrisation of the twist
We will consider twist transformations G which admit two (real) fixed points in 4D. Such transformations are diagonalisable in the 4 representation G 4 . We will consider the case where the eigenvalues λ 1 , λ 2 , λ 3 and λ 4 = 1 λ 1 λ 2 λ 3 are distinct. Accordingly, in the vector representation G 6 will have 6 eigenvectors Y M a , which correspond to 6 fixed points y µ a = Y µ a Y + a after projection to 4D. Due to the basic properties of SO(1, 5) orthogonal matrices, for generic eigenvalues we find that Y a .Y b ∝ (y a − y b ) 2 = 0 for all a, b, except for (a, b) = (1, 2) or (3,4) or (5,6). Since in Euclidean space real vectors cannot be null-separated, this means that one can only have at most two fixed points which are simultaneously real. There are indeed precisely two real fixed points for (2.24) 5 Later we also use the notation θα (rather than ϑα) for inhomogeneities with some extra shifts, see (3.16). This is done so that the limit θα → 0 corresponds to the actual fishnet CFT.
We will denote them as x 0 and x0. To get an explicit parametrisation of the map, one can send x 0 → 0 and x0 → ∞ by a similarity transformation where K is a special conformal transformation acting on coordinates in 4D as and correspondingly, on functions f ∈ F h as 27) and R is of the form with R 4×4 being a 4D rotation matrix. A 4D rotation normally has two invariant orthogonal planes. In a conventional standard position, we can take those planes to be (1,2) and (3,4).
In this case we get (2.29) We will frequently assume that in the standard frame the twist matrix G is of this "diagonal" form where D, S 1,2 and S 3,4 are defined in (2.4). In the representation 4 the twist matrix in the standard frame becomes The eigenvalues λ i 's, satisfying 4 a=1 λ a = 1, are called twist parameters, and will appear throughout the paper.

Global conformal charges
Under a conformal transformation of the wave functions connected to the identity 6 , the transfer matrices are invariant, up to a redefinition of the twist map: Therefore, they commute with the Cartan generators of the form We will consider a basis of common eigenstates of the transfer matrices as well as the Cartan charges, such that and we will denote the conformal charges as Single-valuedness of the wave functions implies that S n ∈ Z, while ∆ can be in principle a generic complex number.
In the next section, we will see how these parameters are naturally identified with the scaling dimensions and spins of twisted operators in the fishnet model. An appropriate quantisation condition and the introduction of the coupling constant will restrict the conformal dimensions ∆ to the physical values in the spectrum of the CFT.

Integrals of motion
The transfer matrices have a polynomial dependence on the spectral parameter, therefore they are parametrised in terms of a finite number of commuting operators, the integrals of motion.
There is also a trivial polynomial dependence in some of the transfer matrices, that does not depend on the state. To separate this part, it is convenient to introduce the fixed polynomials depending on the inhomogeneities: After that the eigenvalues of the transfer matrices can then be written as where P r n (u) is a polynomial of degree n and we use the notation (2.41) The remaining 4J coefficients, denoted as I (a,α) , depend nontrivially on the state. The conformal charges appear in the next-to-highest order terms: where we introduced 7 The remaining 4J − 3 coefficients are the nontrivial eigenvalues of the commuting family of integrals of motion.

The quantisation condition for the wave functions
To define a well posed diagonalisation problem for the transfer matrices, we have to specify an appropriate function space for the eigenvectors, or, in other words, a quantisation condition. The condition we are interested in, for applications to the fishnet CFT, is that wave functions are single-valued, and do not have singularities, except when all the coordinates approach the fixed points x 0 and x0 simultaneously. 8 The approach to these points is obtained by evolving the coordinates with the transformation e −iρQ 0 , for ρ → −∞ and +∞, respectively. Using the fact that the wave function is an eigenvector under such transformation, we find, in the limit, where ∝ e −|ρ| , (ρ → ±∞) is a cutoff measuring the distance from the fixed points. Notice that (2.46) distinguishes between the points x 0 and x0 by the sign of ∆. The interpretation of these equations will become clear when we relate the wave functions to a particular class of correlators in the fishnet CFT where x 0 is the location of a non-protected operator. We expect that the regularity of the wave function everywhere except for x 0 and x0, together with (2.46) at these points, are sufficient conditions to restrict the spectrum of integrals of motion to a discrete set, for any choice of the Cartan charges {∆, S 1 , S 2 }. Alternatively, we can say that the quantisation condition restricts the 4J integrals of motion to a set of curves, parametrised by a continuous variable which can be identified with ∆. In applications to the fishnet CFT, we will alternatively identify the continuous parameter on which the spectrum depends with the coupling constant.
Note that (2.46) implies that there is certain relation between a state with dimension ∆ and a state with dimension −∆, which is obtained by interchanging the fixed points x 0 ↔ x0 with a conformal map. For definiteness we fix this map to be where K is the special conformal transformation appearing in the decomposition of the twist (2.26), andĨ is the "holomorphic inversion" defined as 48) or, in 6D representation, asĨ 6 = diag{1, −1, 1, −1, −1, 1}. Note, however, that whereas this transformation keeps the diagonal form of the twist matrix, it flips the parameters α a → −α a , a = 0, 1, 2. Therefore this transformation inverts the signs of the spins in addition to ∆ → −∆ (see (F.6)), while exchanging the fixed points x 0 and x0. As discussed in section 2.7 and appendix F, this map appears in the relation between left and right eigenvectors of the transfer matrices.

Convergence of the scalar product
Due to the singular behaviour of the eigenfunctions of the transfer matrices (2.46) the convergence of the scalar product is not obvious. The general prescription, which may not be particularly practical, is to make an analytic continuation in parameter space to reach the values of ∆ where the convergence is manifest. This is essentially a ζ-function type of regularisation. This prescription allows to regularise most of the cases, except for the special case of the scalar product of two wave functions with opposite ∆'s in the same frame (i.e. with the same singular points), for example the scalar product of a state and its conjugate. In this case, near the singular point x 0 we get an integral of the type which is log-divergent and thus cannot be regularised in this way. The only option in this case is to introduce the cut-off. In this case the nontrivial object is the coefficient in front of the log, which is likely to be regularisation independent and will play an important role Without analytic continuation in parameters, for generic states Ψ, Φ is divergent near the singularities, however, we expect the expansion in the small regulator to be power-like. In this case we define Ψ, Φ as the finite part in the small expansion, which is a non-ambiguous quantity in the absence of log type of divergences.

The Baxter TQ relations
The Baxter TQ equation gives a complete reformulation for the problem of diagonalising the transfer matrices, as an equation in just one variable. The SoV basis should be possible to obtain by diagonalising B and C along the lines of [7]. Instead of doing this, in this paper we will bypass the explicit SoV construction by studying the TQ relation, which for integrable models with SL(n) global symmetry takes a known form [63,64]. For the rank-3 case relevant for the fishnet theory, it reduces to where T r are the transfer matrix eigenvalues (2.36). After the "gauge" transformation Being a fourth order difference equation, it has four linearly independent solutions that we denote as q a (u) with a = 1, . . . , 4. For what follows, it is also important to introduce the dual Q-functions q a (u) which are built as 3 × 3 determinants times an overall factor, Using the results from Appendix C we find that these functions q a satisfy a 'dual' Baxter equation, which reads We see that it is related to the Baxter equation (2.53) for q a by simply exchanging Q + ↔ Q − , P 4 J ↔ P4 J . The four functions q a are distinguished according to their asymptotic behaviour for large u: where λ a are the twist eigenvalues, andM a are charges defined aŝ (2.57) To fix the basis completely, one should furthermore impose analyticity conditions, as we discuss in section 2.6.1.
The large u asymptotics of q a reads, similarly to (2.56) with the normalisation coefficients fixed by (2.56) as follows (2.59)

Quantisation condition for the Q-functions
To complete the solution for the spectrum, one has to further constrain the analytical properties of the Q-functions, which will lead to the quantisation condition. This condition is the Q-function counterpart of the requirement that the wave functions have no singularities except for the ones in (2.46) we discussed above. The quantisation condition for Q-functions is known at least in the case of weights h α ∈ {1, 2} which are relevant for applications to the fishnet model [49,65], and in this subsection we describe how it works.
Here we also verified it for generic h α . Once it is imposed, the values of the integrals of motion are fixed to a discrete set, for any given assignment of the twists, inhomogeneities and conformal charges (∆, S 1 , S 2 ).

Analytic properties of Q-functions
Solutions of maximal analyticity. Let us define bases of solutions to the Baxter equation (2.53) with the largest possible region of analyticity.
One possibility is choosing the Q-functions to be analytic for very large positive Im u. Such solutions will be denoted as q ↓ (the direction of the arrow indicates that one starts from large positive Im(u), and iterates the Baxter equation to move down). The form of the Baxter equation implies that they unavoidably have infinitely many single poles extending into the lower half plane, located at u ∈ P ↓ , with which implies that as we decrease Im u from a large positive value, the first poles will occur at u = ϑ α − i hα−3 2 , 1 ≤ α ≤ J, and the other poles will be reached by shifts of −i. An alternative basis of solutions, denoted as q ↑ a , can be chosen by requiring analyticity for large negative Im u. The singularities of these functions are single poles for u ∈ P ↑ , These two bases are uniquely defined once their asymptotics -in the respective domains of analyticity -are specified as follows The coefficients B a,n in the above asymptotic expansions are the same for both sets and can be found systematically from the TQ-relations. The powersM a are defined in (2.57).
A useful property is that the Wronskians built with the full set of solutions are explicitly fixed functions of the spectral parameter, generalising the case of differential equations where they are constants. In the present case, defining we have as a consequence of the Baxter equation (2.53): Finding the solution with appropriate analytic properties, we find (with the same constant proportionality coefficient), so we see that the Wronskian is a state-independent function of u.
For the solutions of the dual Baxter equation we similarly define two bases These functions have poles in the sets P dual For any generic choice of the parameters entering the Baxter equation (including the values of the integrals of motion), q ↓ a and q ↑ a can be computed numerically with arbitrary precision with the method of [66] (see [46] for the fishnet case), which we review in section 4.2. To fix the values of the integrals of motion to a discrete set (depending on the continuous parameter ∆) we need to impose the quantisation condition described below.
Relating the two bases.
Since q ↓ a and q ↑ a are both bases of solutions of the same difference equation, they must be related by a matrix of i-periodic functions: This matrix of coefficients can be constructed explicitly in terms of the Q-functions, as a solution of the linear system (2.71): .

(2.72)
Notice that the denominator in (2.72) is related by an argument shift to the quantum determinant q 1234 given by (2.64). From (2.72), it is easy to see that, for generic positions of the inhomogeneities, Ω has first order poles at all points of the form −i(hα−1) 2 + ϑ α + iZ, and no other singularities. Since we have the periodicity Ω(u) = Ω(u + i), we have the pole decomposition It is also useful to analyse det Ω. For that we can repackage the equations (2.71) as and take the determinant of the two sides to get q 1234 , so that det Ω can be fixed explicitly from (2.66), (2.67) and reads For integer weights h α ∈ Z (which includes the case of the fishnet model when h α are 1 or 2), this implies that Ω has unit determinant.
Since the dual Q-functions are defined as determinants, we also have the following relations (2.77)

Quantisation condition
The quantisation condition found in [49,65] simply reads: It has been previously verified for the cases when h α take values 1 or 2, and we have also tested it numerically for non-integer values of h α in the vicinity of 1, so we expect it to be applicable for rather generic h as well. We believe that (2.78) is completely equivalent to the quantisation conditions for wave functions described in section 2.4. Indeed, this can be established for the J = 1 case where we constructed the SoV map explicitly in [49]. It should be possible to prove it in general using the operatorial SoV formalism. The conditions (2.78) correspond to setting to zero 4 × J coefficients entering the pole decompositon (2.73) 9 . It can be satisfied by tuning appropriately the values of the integrals of motion entering the Baxter TQ relation as parameters. More precisely, we find that, for fixed twists and inhomogeneities, (2.78) admits several families of solutions, each family described by a one-parameter flow in the space of the integrals of motion. As already mentioned, in the case of the fishnet theory the continuous parameter can be identified with the coupling constant.

The gluing matrix.
When (2.78) are satisfied, the matrix Ω has a further enhanced symmetry: there exists a constant matrix (called the gluing matrix ) of the form (where the constant γ is nontrivial and depends on the coupling), such that ΓΩ = Ω T Γ . (2.80) The existence of this matrix reflects the structure of the Quantum Spectral Curve of the parent theory N =4 SYM, as discussed in [49], namely the antisymmetric i-periodic matrix ω = (ΓΩ) −1 also appears in QSC. We checked numerically that this symmetry holds also for some generic values of h, beyond the fishnet model. It would be interesting to find a direct proof for the existence of the gluing matrix, using only the conformal spin chain setup. Some nontrivial consequences of (2.80), relevant for this purpose, were found in [65] and are reviewed in appendix D.
Useful identities. We notice that one can use the gluing matrix Γ ab to raise indices and Γ −1 ab to lower them. For what follows it will be convenient to introduce the notation 10 These newly introduced functions satisfy the same Baxter equations as q a and q a respectively. They only differ by the way they transform between the UHP and LHP bases i.e. they transform in the same way as q a in (2.71) but with arrows reversed. It is convenient to introduce where as a consequence of (2.80), ω is anti-symmetric when the quantisation condition is satisfied.

Orthogonality properties
In this section, we discuss orthogonality properties for the eigenstates of the integrable system and relations between left and right eigenvectors of monodromy matrices. Given the scalar product (2.16), we can construct the "transpose" of the transfer matrices (T r ) T with respect to this quadratic form, such that for any f, g ∈ F {h} . As we know, all conformal generators are antisymmetric w.r.t. this bilinear form, so that (T r ) T differ fromT r by the replacement q M N → −q M N or, as we show in appendix F, this is equivalent to first exchanging x 0 and x0 with the holomorphic inversion (2.47), and then also reversing the order of particles in the chain. As the TQ relations are not sensitive to the order of particles this shows, as expected, that the right and left eigenvalues of all the integrals of motion are the same (as would be in the finitedimensional case), but the explicit relation between the left and right eigenvectors in the case of general ϑ α 's and h α 's could be complicated (see an example of the map interchanging particles in some particular cases in [46]). In the simplest case when all inhomogeneities and the weights are equal, the left and right eigenvectors of the integrals of motion are simply related by the holomorphic inversion (2.47) and relabelling coordinates x α → x J−α+1 . With this insight in mind, we should have that the left and right eigenstates corresponding to different eigenvalues of the transfer matrices are orthogonal, just like in the finite-dimensional case. At the same time, if the eigenvalues are equal then we are in the situation discussed around (2.50) where we concluded that the scalar product is log-divergent and the meaningful object is the finite coefficient in front of the logarithm. Furthermore, in the generic situation the spectrum is expected to be non-degenerate.

Fishnet model as a spin chain
The goal of this section is to relate the abstract spin chain formulation with the field theory observables. More precisely, we relate it to the planar fishnet CFT desribed in the Introduction and defined by

Operators and wave functions
The setup we are using in this paper includes the so-called colour-twist operators, which generalise the traditional local single-trace operators. For the details of the construction we refer to [49]. Below we briefly review the construction.
Colour-twist operators. We consider the colour-twisted version of local single-trace operators, introduced in [49]. This gives a deformation of the concept of a local operator, depending on a symmetry transformation G, which in our case will be a generic conformal transformation with two fixed points. The deformation removes degeneracies in the spectrum and as a result makes the SoV construction more regular and uniform. The presence of the twists is also technically convenient and will make the form of our results much more transparent.
Perturbatively, twisted single-trace operators are constructed as 2) where the dots stand for possible mixing with similar operators and The marker T G indicates the starting point of a "twist-cut" on the worldsheet of the planar Feynman graph. Each propagator crossing this line gets deformed in accordance with the conformal group element G. For details and examples of the construction see [49].

Symmetries.
The correlators involving twisted operators transform in a more complicated way under the conformal symmetry, since one should also keep track of the transformation of the colour-twist. Under the conformal transformation C ∈ SO(1, 5), the twist map transforms as G →G ≡ CGC −1 , see Appendix E. This in particular means that all the usual degeneracies of the usual local operators are lifted -in particular the operators transform nontrivially under the translation symmetry removing the necessity of considering conformal primaries and descendants separately. Since the twist G remains invariant under 3 generators in the Cartan subalgebraQ 0 ,Q 1 ,Q 2 defined in (2.32) the twisted operators can be classified by the corresponding 3 quantum numbers Q a , which are two spins and the scaling dimension.
As we will see the scaling dimensions ∆ of such operators can be computed nonperturbatively using integrability and are nontrivial functions of the coupling constant ξ, as well as the twist angles λ a . We can recover the value of the scaling dimension for a standard local operator by taking the untwisting limit λ i → 1.
In addition, composite operators of the fishnet theory carry quantum numbers J, M ∈ Z, associated to the U (1) × U (1) charges of the two scalar fields of the fishnet model. The fishnet Lagrangian admits a Z 4 discrete symmetry generated by the elementary In accordance with the connection with the conformal spin chain (reviewed below), we will refer to J as the length of the operator (i.e. the number of φ 1 fields 12 ), and M the number of magnons (i.e. the number of φ 2 fields minus the number of φ † 2 fields). The case M = J needs some special treatment (see appendix D), and for simplicity in this paper we restrict the analysis to the choice of quantum numbers |M | < J.
It is known that there also exist infinitely many non-dynamical operators in the fishnet theory, which have zero anomalous dimension at the planar level. In this paper, we will mostly focus on the operators O with non-zero anomalous dimensions.
Finally, notice that, since the fishnet theory is non-unitary, the dilatation operator is strictly speaking not diagonalisable. Indeed, examples of Jordan blocks of the anomalous dimension matrix have been studied in [67][68][69]. In all known examples it was found that the generalised eigenvalues of the non-diagonalisable Jordan blocks are zero, which can be verified at all loops in some cases by diagrammatic arguments [67]. It is our working assumption that this is true in general (at least at the planar level), so that the only operators acquiring anomalous dimensions are proper eigenvectors of the dilatation operator, which do not mix with the non-diagonalisable Jordan blocks. We believe that all such nontrivial operators are captured by the integrability methods discussed in this paper.
The CFT wave function. All information on a single-trace operator can be encoded into a special kind of correlator called the CFT wave function [48], which will play the role of the the spin chain wave function. It is defined as the renormalised correlator where the indices I α ∈ {0, 1, −1,0}, and the fields in the second trace are In order for the correlator (3.3) to be non-zero, these indices have to satisfy M = J α=1 I α , where M is the number of magnons in the operator O. We see that there is a weight The conformal charges of the operator O can be extracted from the CFT wave-function by acting only on the x 1 , . . . , x J coordinates. Namely, the charges defined as in (2.34): are related to the scaling dimension and spins of the operator: 12 More precisely, J is the number of φ1 fields minus the number of φ † 1 fields. However, we are interested in the case of non-protected operators, which requires only φ1 to be present, as otherwise no Feynman diagrams are possible beyond tree level. This follows from the conformal invariance of the correlator, which involves transformation of all J + 1 local operators in the correlator defining the CFT wavefunction. So by transforming only the J protected legs, we indirectly measure the quantum numbers of O.
We see that there is a clear analogy between the CFT-wave function and the wave function of the conformal spin chain (or fishchain) which we introduced in the previous section. In order to complete the embedding, one needs to show that the CFT wavefunction diagonalises the transfer matrix for some values of inhomogeneities ϑ α . This was done in [46], and in the next section we will briefly review the main elements of the construction.

CFT wave functions as integrable fishchain eigenstates
In this section we review the construction of [46] which shows that the CFT wave functions are also eigenfunctions of the integrable spin chain defined in section 2.3. The basis for this consideration is the Dyson-Schwinger equation satisfied by the wave functions. For simplicity we review the case M = 0 and then present the general result.
The Dyson-Schwinger evolution equation. Based on the Feynman rules of the fishnet theory, the wave functions for physical operators are built out of Feynman diagrams which form an infinite ladder. Each rung of the ladder is an iteration of a graph-building 13 The first Jacobian under the integral comes from the rules in [49], where the twisted propagator is given by ∂G ∂y operatorB. This is an integral operator acting on a function of J variables as 13 [42,67]: (3.10) Twisted periodic boundary conditions imply that y 0 ≡ G • y J , where G is the twist map, corresponding to the convention illustrated in figure 2 for the twisted trace.
The ladder structure shown in figure 2 implies that the CFT wave function satisfies a Dyson-Schwinger evolution equation: The integral operatorB can be shown to commute with the Cartan generatorsQ a . Furthermore, we will see that the integral operatorB belongs to a family of mutually commuting operators.
A useful property found in [46,48] is that even for the case with magnons, for |M | < |J|, the inverse of the graph-building operator is a differential operator. To demonstrate this property for M = 0 case we notice that (with boundary conditions x α−1 = G(x J ) for α = 1) and the Dyson-Schwinger equation can also be written aŝ which is a partial differential equation. Similar relation holds for the case with magnons too. Note that (3.13) is satisfied everywhere except for the singular points x 0 and x0, where the behaviour of the wave function is dictated by the conformal dimension (2.46).

Integrability.
At the heart of the integrability of the fishnet CFT is the fact that the HamiltonianĤ (and its generalisation for the case with magnons) is embedded in the family of transfer matrices introduced in section 2.3. It was proved in [46] that (3.14) For (3.14) to be true one should properly adjust the values of h α and ϑ α . The correct values of h α and ϑ α depend on the field content of the given site [46]: As a consequence of (3.14) the Hamiltonian (i.e., the inverse graph building operator) commutes with all the integrals of motion for the fishchain, and we can use all the power of integrability to diagonalise it. Even though (3.14) requires a special tuning of the inhomogeneities (3.15), we will still consider the general case in order to have more integrable parameters. In the following we introduce a convenient notation for the shifted inhomogeneities: so that the undeformed Hamiltonian is recovered when all θ α → 0.
In the light of (3.14), we see that the Hamiltonian is identified with the integral of motionÎ (0,1) . For physical states solving the Schrödinger equation (3.13), we can relate its eigenvalue to the coupling constant This is the key equation for introducing the 't Hooft coupling into the integrability formulation [46,48,65]. The spectrum of this spin chain is somewhat unusual in comparison to more standard representations with highest weight. Usually the spectrum is discrete -only some particular values of the integrals of motion can appear in the spectrum. In the present case for any fixed ∆ and integer S 1 , S 2 we find an eigenfunction satisfying the correct quantisation condition (2.46) so that the eigenvalues of the integrals of motion I (n,α) (∆, S 1 , S 2 ) are in general multi-valued functions of ∆. The physical discrete spectrum ∆ n (ξ) is then found by imposing the condition (3.17).
The numerical method described in [46], based on the TQ relation and the quantisation condition from section 2.6.2, allows one to find the values of ∆ n (ξ 2 ) with high precision.

The meaning of the spin chain norm in CFT
In order to complete the embedding of the correlators of local operators into the general fishchain framework, we have to give an interpretation of the scalar product of two CFT wave functions, according to (2.16). To understand its meaning let us consider a particular case of two operators of length J = 6 with one magnon spiralling around M = 1 as in fig. 3. The scalar product (2.16) is given by the product of these two wave functions with powers of d'Alembertian in between integrated over J variables. In the case of figure 3, all weights h α are equal to 1 except for h 2 = 2, which means that for all points we get d'Alembertian to the first power except the site with the magnon, where no d'Alembertian is added. Note that this is perfect from the Feynman diagram perspective, as we get one of two black propagators annihilated by the d'Alembertian at non-magnon sites, while where the magnons are sitting we glue them into the fishnet interaction 4-vertex. Thus the scalar product gives a correlation function of two twisted operators.
The case when the operator B is related to the operator A by the holomorphic inversion (2.47) is particularly important. Naively we get a two-point correlator in this case, however, as was discussed in section 2.4.1 the scalar product is log-divergent, unlike the two point correlator which should be finite for two normalised operators. The reason for this discrepancy is that the combinatorics of the diagrams is a bit different: for example there are two diagrams with one wheel in the scalar product and only one in the two point correlator. As was shown in [48] and we review in appendix G, the relation is as follows: where in the r.h.s. we get the two-point function: which has the same kinematical dependence as in any CFT, even in the presence of twists [49] (see Appendix E). The coefficient of the two-point function itself can be understood -at least formally -as a scalar product between a nontrivial state, and the same state extrapolated to zero coupling. 14

Baxter TQ relations in the fishnet theory
In this section we speciallise the results of section 2.5 to the fishnet case, where h α is either 1 or 2.
Assuming M ≥ 0, we can always get rid of the anti-magnons by redefining the Qfunctions as was discussed in [46]. For the sake of clarity, we will pick a particular order for the magnons, and take I i = 1 for 1 ≤ i ≤ M and I i = 0 for M + 1 ≤ i ≤ J (different orderings are related by a similarity transformation [46]). In this case, the Baxter equation (2.53) can be written explicitly as containing the shifted inhomogeneities (3.16). We denote the 4 solutions of the above equation as q a (with lower index). When θ α → 0 we get the actual fishnet Baxter equation of [46]. We will keep θ α arbitrary as an important regulator. In the simplest case of no magnons, and with θ α set to zero we get The dual Q-function, constructed from q a 's as a 3 × 3 Wronskian. In the case of the fishnet, it is convenient to introduce an extra multiplier w.r.t. (2.54), J α=M +1 (u − θ α ), which simplifies the form of the dual Baxter equation for this particular choice of weights: The dual Baxter equation has a form simply related to (3.20) by the interchange of the polynomials P4 J (u) ↔ P 4 J (u). Like before in section 2.6, there are two distinguished bases for the Q-functions: q ↑ a that are lower-half-plane (LHP) analytic and q ↓ a that are upper-half-plane (UHP) analytic, and similarly for the dual q a 's. The positions of the poles (2.60) and (2.61) simplify to 23) and the same sets describe the poles for the dual q a↓ and q a↑ . This is indeed consistent with (2.69), taking into account that q a in this section are defined with an extra factor J α=M +1 (u − θ α ), which cancels some of the poles in P dual . Notice that, for real inhomogeneities, the poles of q ↓ are in the lower half plane or at most on the real axis, and the poles of q ↑ are all in the upper half plane.

SoV scalar product cookbook
Even though in this paper we have not attempted building the SoV basis, along the lines of recent papers [7,13,15] a lot of structure of the result can be predicted based on simple observations such as analyticity of Q-functions and TQ relations. In this section we describe one of the main results of the paper -the bilinear form built from two Q-functions, which should be a building block for the scalar product in the SoV basis. We will also see how the quantisation condition described in section 2.6 naturally appears from the arguments of this section.

Bilinear forms of Q-functions
Here we introduce a set of bilinear forms pairing two Q-functions. The construction is inspired by the findings in [23]. We introduce the following notation for the bilinear form of two Q-functions, in general for different two states A and B: where p a 's are defined in (2.81), µ(u) is a meromorphic i-periodic function (without poles on the integration contour),Ô is a finite difference operator of the form with polynomial coefficients o i (u), and the integration contour is defined as follows so that it contains all poles of the Q-functions. Furthermore, in section 4.3 below we show that, among all µ(u)'s so defined, there are only J distinguishable periodic functions µ α , meaning that p B aÔ • q A b µ for any µ(u) can be expressed as a linear combination (with state-independent coefficients) of p B aÔ • q A b µ β . It is convenient to take the following basis which has zeroes at all u = θ α + iZ except for α = β where it is equal to 1. For simplicity we denote In general the integrals (4.2) are not convergent, and have to be regularised by a ζfunction type of regularisation. We will give an explicit analytic prescription on how to do that in the next section, and also explain how to compute this type of integrals numerically.

Regularisation
A convenient way of computing, or more precisely defining, the integrals of the type (4.1) is by re-expressing them as a sum over residues. Assuming for simplicity that the finite difference operatorÔ does not contain the shift D −2n with n > 1, the integrand of p B↑ aÔ • q A↓ b µ β has simple poles at θ β + iZ due to the poles in the q-and p-functions. 15 In the UHP the poles are coming from p B↑ , to make them manifest one can use (2.83) and similarly q ↓A produces poles in the LHP, whose residues are obtained as res u=θ β −in q ↓ a (u) = ω β,ac p c↑ (θ β − in) , n > 0 , (4.6) so that we can write formally forÔ = u kD2m , For the cases with shift m < 0 and for the states with magnons, there could be a finite number of double poles contributing. In the observables we consider below, the relevant case is m = −1, which gives a double pole contribution at the origin -in this case one should add res to the r.h.s. of (4.7). We will discuss such terms in more detail in the next section.
Depending on the parameters of the states A and B, the UHP and LHP sums in the first and second lines of (4.7) may or may not converge in the usual sense. Indeed, at large n we can use the asymptotic expansion (2.62), so that the combination we need reads using that the formal sum ∞ n=1 λ n n α = Li α (λ) (4.10) is well defined for all α and λ, except for λ = 1, where it reduces to the ζ-function which has a pole at α = 1. 16 This means that for generic λ A b and β B c the sum is well defined in ζ-regularisation. In section 4.2 we also describe how this sum can be computed very efficiently numerically.
We see that for two generic states A and B all the Q-bilinear forms are well defined. However, in some particular situations, and especially when A coincides with B there is a danger to hit the singularity of (4.10), when one of the terms behaves as 1/n at large n's. In the next section we will be considering the diagonal form factors for which one has to face this problem.
Note that in fact the scalar product of two CFT wavefunctions is also log-divergent when ∆ A = ∆ B as we discussed in section 2.4. So the divergence we found in the Qbilinear forms is consistent with the divergence in the coordinate representation of the scalar product. However, it is still meaningful to define the finite part in front of the log-divergence described in section 2.4 in terms of the Q-functions as well.
In the next section we will show that there is still a sufficient amount of the bilinear combinations of Q-functions which are finite even when the states A and B are identical. We will then show in section 5 that those combinations can be used to express a vast class of diagonal form factors as a determinant of the finite Q-bilinear forms. Intriguingly, we will see that those combinations are finite only when the Q-functions obey the correct quantisation condition from section 2.6.2 (in addition to the Baxter TQ relations).

Conjugation under the Q-bilinear form
An important property of the Q-bilinear form is that it allows to define conjugate finitedifference operators O † , independently of the measure µ(u). For example, for the elementary operatorÔ we define its conjugateÔ † bŷ As the operators of this type, which can be also written asÔ = (u + im 2 ) kDm , constitute a basis we thus have defined (by linearity) the conjugation for all finite-difference operators with polynomial coefficients.
The key property of the Q-bilinear forms, which we will utilise in the next section, is the following: which can be seen by shifting the integration contour by −im along itself.

Finite diagonal combinations
We will now show that the quantisation condition (2.78) can be reinterpreted as the requirement that some Q-bilinear forms are finite. As we explained in the previous section the problem arises whenever the twist parameters get cancelled and also the powers combine so that we get a 1/n term at large n's. For identical states A = B, the dangerous terms are when in (4.9) we get b = c, then both twists and the non-integer parts of the powersM a get cancelled. Then the convergence will really depend on the value k − D 0 , but we will see that we will need the values of k such that k − D 0 ≥ −1, meaning that divergence will occur in the sub-leading 1/u terms. So the best way to ensure convergence is to avoid b = c. Let us show that the following combinations are finite: (with no summation over a). Indeed, on-shell the combination ω = Ω −1 Γ −1 is antisymmetric, as was discussed in section 2.6.2, and thus ω β,ab should be anti-symmetric in a ↔ b as well. From that observation, it immediately follows that in the pole expansions on the r.h.s. of (4.7) we never get both q's (or p's) with the same indices, ensuring that the sum is well defined with our regularisation.

Numerical evaluation of the Q-bilinear forms
Here we explain how the Q-bilinear forms can be evaluated numerically in practice. It would be useful to recall first how we solve the TQ-relations and impose the quantisation condition. First, we find q at large u by plugging the asymptotic expansion (2.62) into the Baxter equation (3.20). We usually keep around 20 coefficients of B a,n , B a n to get a very good (around 50 digits) precision for q at the values |u| > 100. After that we use (3.20) to descend from large u's along the imaginary axis to a finite value of u. This way we can compute both q ↓ a and q ↑ a , with the difference that for the former we use (3.20) to move down from large positive Im(u) values, whereas for q ↑ a we use (3.20) to move up from the asymptotic domain far below the real axis. Then we compute Ω using (2.72) at J + 1 different points, which then allows us to deduce Ω β , β = 0, . . . , J from (2.73). Then we have to adjust the coefficients in the polynomials P , which are the eigenvalues of the integrals of motion, to impose the quantisation condition Ω T β = −Ω β , which then gives all integrals of motion fixed as a (multi-valued) function of ξ and also allows us to determine the gluing matrix Γ (which is parametrised by one constant γ from (2.79), fixed by (2.80)). This procedure was already well established in the previous works [44,46,49] and is a simplified version of the similar method [66] developed for the spectrum of N = 4 SYM.
Next, in order to evaluate the bilinear combinations of the Q-functions we start from the expression (4.7). Consider the first term for example: (4.14) In principle we know already all ingredients such as ω β and q which we can compute at any point. However, as we discussed previously, the sum (4.14) does not necessarily converge in the usual sense and a ζ-type of regularisation is needed. What we do in practice is fix some cut-off N (in practice around 100) and split the sum into two parts. For n < N we compute the sum without any further approximation, for n > N we replace q s by their asymptotics (2.62) and re-expand it for large n, so that it takes the form ∞ n=N Λ n n −α 1 + C 1 n + C 2 n 2 + . . . .  We found that this procedure is very accurate and allows one to compute the Q-bilinear forms with 20 digit precision quite easily. We will give some examples in section 5.4.3.
In addition to the infinite sums described above, in some cases we might need to evaluate the contributions of the double poles (4.8). These contributions may potentially require evaluation of derivatives of Q-functions. In the examples we considered, the second 17 Φ is HurwitzLerchPhi in Mathematica. order pole cancels and becomes a single pole, so the term (4.8) can be evaluated on equal footing with others.

Completeness of the basis of the Q-bilinear forms
Here we discuss why the basis of periodic functions (4.3) is complete in the sense described below. We will assume that the double poles do not appear.
First, we notice that the following periodic function will result in a zero for the integral of the type (4.1) as µ 0 (u) will cancel all residues inside the integration contour. Then we notice that any polynomial in e 2πu can be represented asP (e 2πu )µ 0 (u) + A β µ β (u), for some polynomialP (u). The first term, divisible by µ 0 (u), gives zero and can be removed. Similarly, any i-periodic function F (u) with no poles inside the integration contour can be written as where R(u) is an i-periodic function, also analytic inside the contour. The last term does not give any contribution and thus this case also reduces to our basis. More complicated is the case when the periodic function F (u) has poles inside the integration contour. To discuss this case, let us assume that the twists are phases λ a = e iφa such that φ B a > φ A b meaning that the integrand, which behaves as , is decaying to the right of the contour. In this case the function which has a pole at u = u 0 + iZ, will give zero under the integral (4.1). This is because the integrand decays exponentially for large positive Re u, due to the asymptotics of q and p, and also at large negative Re u, due to the exponential decay of the function µ(u) itselfhence the contour can be pushed to infinity giving a zero result. For the opposite φ B a < φ A b case, we rewrite where now the first term will be zero under the integral and the second term is regular. In conclusion, we see that one can subtract simple poles in the measure factor F (u) in terms of such functions which do not contribute to the Q-bilinear forms, bringing us back to the case with no poles. This shows that the basis K AB ab,α of Q-bilinear forms is the most general. In the next section we will show how to use it to compute nontrivial correlation functions in the fishnet theory.
Note that the above consideration also allows us to relate the basis of forms described above, to the one defined with opposite arrow combinationsK ab,α ≡ p ↓ aÔ q ↑ a µα . For that, one can use the matrix Ω to reverse the direction of the arrows (2.76). Since Ω is i-periodic, at the level of the integrals it can be replaced by a combination of the elementary measures µ β (the coefficients are state-dependent in this case). For some observables, a basis obtained with a certain pattern of the arrows might be more natural, and lead to simplifications. One such example will be the g-function, which we consider in section 7.

Functional SoV for a class of correlation functions
In this section we show how the Q-bilinear forms, introduced in the previous section, can be combined into physical observables of the QFT.
In particular, we focus on a certain class of observables which can be obtained as a result of variation of the numerous parameters (such as inhomogeneities, twist angles and coupling). Those observables can be studied within the functional SoV approach introduced in [15] and do not require further input from the more conventional operatorial SoV method. Following [15], the variations of the integrals of motion are found as solutions of a linear system of equations.
Like in the case of spin chains [7], these results are expected to provide deep structural insights on the form of the result for more general observables, including the off-diagonal cases. In particular, they should fix explicitly the measure factor in the SoV basis, and lead to the determinant form for an even wider class of correlators.
Among these observables, we will discuss in detail the explicit expression for the diagonal 3-point functions involving the Lagrangian. We also describe the form factor of the variation of the Hamiltonian w.r.t. local weights which amounts to a nontrivial local insertion at the level of diagrams.

Conjugation property of Baxter equations
The key starting point for the functional SoV consideration is the Baxter equation satisfied by the Q-functions q a and also the dual Baxter equation, satisfied by q a or p a as defined in the previous section in (2.81).

Baxter TQ-relation as a finite difference operator
Following [15], we introduce a useful notation, and represent the Baxter equation in terms of a finite difference operator B. To define it we use the shift operators D ≡ e i 2 ∂u acting on a function of the spectral parameter as D • f (u) = f (u + i 2 ). The Baxter finite-difference operator is given by where the state-dependent factors I (b,α) are the 4J nontrivial eigenvalues of integrals of motion, D b are state-independent finite difference operators defined as and P θ , Q θ , R θ are fixed polynomials: With these definitions the Baxter equation (3.20) is written as Similarly, the dual Baxter equation satisfied by the Q-functions with upper indices q a is given by To see this, it is sufficient to notice that D † −4 = D +4 and D † −2 = D +2 . The property (5.6) is similar to what was observed in [15] and will be used intensively in the next several sections.

Variation of the Baxter equation
We now consider the variation of the integrals of motion of a physical state with respect to a tunable parameter p. A natural application is when this parameter is the coupling constant, but we can consider also varying the twists or inhomogeneities which enter the definition of the Hamiltonian. The most general variation is a combination of all these.
Any type of variation induces a change in the Q-functions q a → q a + δq a , and simultaneously in the constants appearing in the Baxter equation. In general we have so that the Baxter equation to linear order in the variation reads (B + δB) • (q a + δq a ) = 0 . (5.8) Inserting this condition inside the brackets with the dual Q-function p a (to ensure the bracket is finite we should make sure the indices coincide) we find 0 = p ↑ a (B + δB) • (q ↓ a + δq ↓ a ) µα , and expanding we find where the first term on the lhs vanishes due to (5.4) and (5.6), and we keep only the first order in the variation. Using the basis of measures µ α discussed in the previous section, these are precisely 4J linear equations as a = 1, . . . , 4 and α = 1, . . . , J, for 4J variables ∂ p I (b,β) . It is convenient to rewrite the linear system as the matrix of coefficients is defined by blocks as for 1 ≤ a ≤ 4, b ∈ {2, 0,0, −2} and 1 ≤ α, β ≤ J. The r.h.s. δ V contains the part of the variation for fixed value of I: In the next sections we consider some particular cases revealing additional features in comparison with the HW spin chain cases.

Zero mode of the variation matrix
An important observation, which we will further explore here, is that the matrix M is in fact degenerate, and in the generic situation should have one null vector. As we will see, the zero mode is related to a physically very important case of the variation in parameters -where we vary the coupling constant for fixed values of the twists and inhomogeneities. This is particularly interesting because even in the non-twisted theory the quantity ∂ ξ 2 ∆ gives a nontrivial structure constant [51].
To understand how that is related to the existence of the null vector, we recall that both ∆ and ξ are two conserved charges. Thus, in order to compute ∂ ξ 2 I a,α we are not actually changing any parameters of the system, we just follow the one dimensional space of solutions of the same spin chain. Due to this, we have the inhomogeneous part set to zero δ ξ 2 D b = δ ξ 2 V = 0, therefore (5.10) reduces to a homogeneous equation: M · ∂ ξ 2 I = 0 .
(5.14) Equation (5.14) implies that there is a one-dimensional null space. The null vector ∂ ξ 2 I has to be normalised so that its (0, 1) component is fixed according to (3.17) to be After that, the other components are fixed (according to our numerical experiments) and give indeed the variation of all other integrals of motion (which also include ∆) with respect to the coupling constant ξ.
The observation that the matrix M has rank 4J − 1 is based on a large variety of numerical tests, but the analytic proof seems to be complicated. We leave this problem to future studies. We also note that ∆(ξ) is a multi-valued function and has branch cuts at various nontrivial values of ξ. Of course, at the branch points the derivative ∆ (ξ) diverges, which should manifest in the vanishing of the (0, 1) component of the null eigenvector of M.

Explicit result for ∂ ξ 2 ∆
For SoV applications it could be useful to express the result for ∂ ξ 2 ∆ explicitly as a ratio of determinants -something which the operatorial SoV is expected to give like in [7].
To make the solution explicit, we turn the homogeneous linear system into an inhomogeneous one by moving the column "(0, 1)" to the right hand side of the equation. Using the normalisation (5.15), we find: which is a system of 4J equations in 4J − 1 unknowns. One equation out of 4J is linearly dependent so there is a unique solution. Furthermore, some of the variables (∂ ξ 2 I) (b,β) are related to each other in a simple way, since due to (2.44) we have Using these relations and selecting 4J − 3 equations out of 4J we can get a non-degenerate non-homogeneous linear system. The solution can be written as a ratio of determinants where N, D are determinants of (4J −3)×(4J −3) dimensional matrices, where each element is a Q-bilinear form.

Variations in other parameters
As discussed above in section 5.3, we can also consider more general variations, with respect to an external parameter such as the twist angles or inhomogeneities. In this case, we obtain an inhomogeneous linear system (5.10) with a nonvanishing r.h.s. and the same matrix of coefficients M as for the ξ variations. Since M has rank 4J − 1, such a linear system has a one-parameter family of solutions, which differ by the null eigenvector of M. As we saw in the previous section, the null vector represents the variation with respect to the coupling constant. Depending on which observable we are computing, this ambiguity should be fixed accordingly. For instance, in order to compute variation of the scaling dimension w.r.t. the twist, we would keep the value of the coupling fixed (as the graph building operator is well defined regardless of the twist value). This means that we should set the solution's component "(0, 1)" to zero, which fixes the solution uniquely. Below we will also consider variation of the Hamiltonian w.r.t the weights h α , in which case as we will discuss the ambiguity is fixed by choosing appropriate values for the leading nontrivial integrals of motion (2.44). By choosing a subset of 4J − 1 among the original 4J equations, we can write the solution as a ratio of two determinants of size (4J −1), built from Q-bilinear forms, similarly to (5.17).

Numerical test
We verified the above general formalism numerically. For instance, we studied the case 7 12 , and coupling ξ 2 0.100992, we considered a solution corresponding to a state with zero spins, and ∆ 1.93034. We computed numerically the coefficients of the linear system, 18

Variations as spin chain form factors
At the operator level variations of the integrals of motion w.r.t. a parameter correspond to nontrivial spin chain expectation values, as was understood in [7,13,15] for rational gl(N ) spin chains. Let us briefly recall this argument here. Consider a left and a right eigenstate of the transfer matrix Ψ L and Ψ R corresponding to the same eigenvalue. When we change the values of the parameter p → p + δp, these wavefunctions will also change, as will the transfer matrix eigenvalues, so we have Ψ L , (T r + δT r ) • (Ψ R + δΨ R ) = (T r + δT r ) Ψ L , (Ψ R + δΨ R . (5.24) Expanding it to the first order in the variation and using that Ψ L is a left eigenvector of T r , we see that several terms cancel and we are left with So we see that the diagonal form factor of the nontrivial operator δT r δp is written in terms of the variation of the eigenvalue, or, equivalently, of the integrals of motion. As a result, like discussed in section 5.4.2, this form factor can be written in terms of determinants of the type (5.17) built from Q-functions.
The operator whose form factor we are computing can be quite nontrivial. In the next section we study as an example the variation of the Hamiltonian eigenvalue with respect to the local weight.
To link the spin chain picture with the fishnet CFT, we can notice that the right spin chain eigenstates become the CFT wavefunctions (i.e. correlators of the type (3.3)) once we set the inhomogeneities to the corresponding values (3.16). Moreover, as shown in appendix F.3, the left state Ψ L also can be identified with an appropriately chosen CFT wavefunctions that involve a conjugated set of fields, with the conjugated operator sitting at the other fixed point x0 and twisted by the inverse map G −1 . The form factors act as operators on the chain in between the two wave functions. In general, the result can be interpreted as a certain off-shell observable given by diagrams with nontrivial modification or insertion of some propagators and vertices, whose precise form is dictated by the operator we consider. In other words, this way we can compute a class of rather nontrivial Feynman diagrams similar to those appearing in multipoint correlators.

Variation with respect to local weights h α
The variations in the twist angles λ a and the inhomogeneities θ α were studied already in the spin chain context in [7]. In this section we consider another type of variation corresponding to changing the weights h α that define the representation of the conformal group at each site α of our spin chain. It would be interesting also to further explore this type of variation for the usual gl(n) spin chains.
In the case of fishnet CFT, we fix h α to the values 1 or 2, however by taking a variation around these points we will be able to compute rather nontrivial quantities. Below we discuss several subtleties related to these variations and then present a particular nontrivial example -the variation of the Hamiltonian.

Variation with respect to weights and scalar product
As the weights h α explicitly enter the scalar product (2.16), it is not immediately clear that the argument leading to the form factor result (5.26) still goes through. Let us show that the result is still valid. A small variation of one of the weights h α → h α + δh α induces a change in the transfer matrix eigenstate Ψ → Ψ + δΨ, and in the integrals of motion I →Î + δÎ. At the same time, the conformally invariant scalar product (2.16) itself will change infinitesimally, which can be represented by the insertion of an operator where we assumed that the wave functions f , g have no dependence on the weights. Let us then consider a left and right eigenstates of the transfer matrices, Ψ L and Ψ R , respectively, corresponding to the same transfer matrix eigenvalue. Under a small variation of the weights, the eigenvalue equation becomes Expanding it at the first order in the variation and using the same logic as before, we see that there will be extra terms in the r.h.s. and l.h.s. involving the insertion ofÂ α , (with no summation over α). We notice that the extra terms involvingÂ α cancel against each other. Thus the expression (5.26) is still valid.

Fixing the zero mode
Let us also clarify how to fix the ambiguity in solving the linear system (5.14) (coming from the Baxter equation) for the variations of integrals of motion when we vary h α . Keeping the value of I (0,1) fixed in this case is not what we would like to do, since the relation betweenÎ (0,1) and coupling only exists at the specific values of h α corresponding to the fishnet CFT. Instead what we are trying to achieve is to find the diagonal matrix element of the derivative w.r.t. h α of the integrals of motion, which are differential operators with coefficients dependent on the weights h α . We notice that the h-dependence of the first subleading integrals of motion (2.44) arises only due to the linear in h term in the dilatation operator (2.8) (these integrals of motion are simply linear combinations of the Cartan charges) which appears in (2.44) via ∆. Thus we have These are the conditions we impose in order to fix the zero mode ambiguity in solving the linear system, which as a result allows us to compute the variations of all other integrals of motion in terms of the Q-functions.

Example: variation of the inverse graph-building operator
A particular nontrivial operator whose variation we can consider is the HamiltonianĤ, i.e. inverse of the graph-building operator, described in section 3.2 in (3.12), since it is one of the coefficients of the transfer matrix T 6 , namely the integral of motionÎ (0,1) = T 6 u=0 . Let us compute its variation w.r.t the weight on one site h α . We will restrict to the case of the state without magnons, i.e. after taking the variation we set all h α = 1 and θ α = 0. We find δĤ = δh α Tr L 6 J (0) . . .L 6 α+1 (0)V αL where the operatorV is read off from the Lax operator (2.19), in terms of the conformal generatorsq M N defined in section 2.1. The variation ∂ hq comes about due to the explicit h-dependence in the realisation of the generators in (2.4) which leads to the nonzero components of ∂ hq being (5.33) In order to write the operator (5.32) more explicitly, it is very useful to employ the 6D realisation of the conformal group we outlined in section 2.1. In this formalism, theq M N operators do not depend on h explicitly, rather the h-dependence is contained in the 6D function (2.3) on which they act. After a somewhat lengthy calculation (similar to that done in [48] to compute explicitly the Hamiltonian), we find many cancellations and the result for the variation of the Lax operator in h is surprisingly simple, namelŷ We see that the change from the original Hamiltonian (3.12) amounts to removing the Laplacian operator at the site α and replacing the two inverse propagators involving that site with a combination of x's and derivatives. Thus we can compute the form factor of this operator according to (5.26). We can further simplify the result by considering the operator ∂Ĥ ∂hαĤ −1 whose form factor is trivially related to that of ∂Ĥ ∂hα (via multiplication by ξ 2J ) since the states we consider diagonaliseĤ. Nicely, the result is a 'local' operator that acts only on the three neighbouring sites α − 1, α and α + 1, . (5.36) The expectation value of this operator corresponds to a class of nontrivially modified Feynman diagrams that are therefore computable within our functional SoV approach. It would be also interesting to try and link it with concrete OPE coefficients and other observables. In general, using variations in the local spin chain parameters such as the weights and the inhomogeneities opens the way to computing many nontrivial observables, and it would be important to explore them further.

Scalar Products in Separated Variables
In the spin chains with more regular HW representations, it was shown in [7,13] that one can obtain a lot of structural information from the functional SoV approach [15]. In particular one can deduce the scalar product in SoV basis in determinant form. In the current case, when none of the Q-functions are polynomial, there are additional subtleties which can be properly addressed with extra additional input from the operatorial SoV approach, like in [7]. The manifestations of these additional complications can be seen for example in the properties of the scalar product in the coordinate space (see section 2.4.1). A rather new feature, in comparison to the quasi-periodic systems studied with SoV methods so far, is that the scalar product has a log-divergence when both states have the same dimension ∆ (and the twists are opposite); in this case, the meaningful quantity is the coefficient in front of the logarithm. At the same time, for two generic states we expect the properly regularised scalar product to be finite and thus we may be able to follow the usual path of [7,13] and establish the link between the SoV and coordinate representations of the scalar product. In this section we will introduce the SoV-like scalar product, following the functional SoV approach. Then we discuss its key properties such as orthogonality and also show that the log-divergence presented in the coordinate representation emerges naturally in the functional SoV formalism.

General philosophy
First let us schematically repeat the general logic. As in the case of the simple spin chains we expect that there is an SoV basis x| such that the eigenvectors of the transfer matrix (as well as many other states) |Ψ factorise. For convenience in this section we will use the bra and ket notations, but the scalar product in the case of the conformal spin chain should coincide with the one introduced in section 2.2. In general we can write where Q(u) are some simple combinations of the Q-functions q a (u), corresponding to the state A, with some shifts. Similarly, the conjugate states (as defined in section 2.4) are factorised in general in the dual SoV basis |y Then, using the expected completeness of the SoV basis, we get where the SoV measure M x,y ≡ x|y −1 , can be deduced most efficiently from the functional SoV approach, like in [7]. The main point of the expression (6.3) is that it allows one to concentrate all state-dependent information in the Q-functions, and combines it together with some universal measure factor M x,y into the scalar product. Furthermore the r.h.s. of (6.3) can be usually written as a determinant, which makes it very useful in practice. The idea of [15] was that the r.h.s. of (6.3) can be deduced directly from the Baxter TQrelations and is greatly constrained by the requirement that for two different eigenfunctions of the transfer matrices the combination of Q-functions should obey orthogonality, i.e. vanish for any pair of states with distinct eigenvalues.
In this section we will follow [15] to establish the possible form of the r.h.s. of (6.3) and discuss the main properties of the resulting expression. In order to conclude with certainty about the relation to the l.h.s. of (6.3), i.e. to the conformally invariant scalar product (2.16), further studies, in particular at the operatorial SoV side, are required and will be reported elsewhere.

Orthogonality relation for generic states
Here we repeat the argument of [15] for the orthogonality relation. The starting point is the Baxter operator (5.1). We will need two copies of that for two different states A and B such that Using the conjugation property of the Q-bilinear form we have again Writing the difference of the Baxter operators more explicitly we get Combining the two we get J β=1 b∈{−2,0,0,2} where like in (5.12) we have The main difference with the section 5.3, where we were interested in the continuous variation of the conserved charges w.r.t. the internal (e.g. ξ) or external (e.g. θ α ) parameters, is that for the convergence of the Q-bilinear form constituting M AB we are no longer required to assume a and c to be related to each other. So in total we have 16 possible combinations of the indices a and c. In the case of spin chains with HW representation the natural choice for a and c is given by the polynomiality constraint of p a and q c . In that case the set of possible values of a and c is such that the analogue of the matrix M AB is a square matrix. In our case this would correspond to a selection of 4 combinations of (a, c) out of 16. Imagine this was done, then the system (6.7) becomes a homogeneous 4J × 4J linear system on 4J unknowns I A − I B . Following the logic of [15], this is only possible if The expression (6.9) represents the so-called orthogonality relation. Even though the above expression is derived for two different eigenstates of the same transfer matrix, the interpretation (6.9) goes beyond that case. In the cases of the HW spin chains, where the operatorial formalism was worked out explicitly, it was shown [7] that it gives the realisation of the expression (6.3) for the scalar product of two states whereN A , N B are the normalisation coefficients, dependent on one state only. Of course, when |Ψ A and Ψ B | are two different eigenstates of the same transfer matrix then (6.9) should hold indeed. Furthermore, in [7] it was shown that the relation stays correct in a number of less trivial situations, for example when the states correspond to the transfer matrices with different twist eigenvalues λ a , or even when the states are not eigenstates of any transfer matrix, but are factorisable in the same SoV bases x| and |y . In all these cases (6.10) gives a nontrivial r.h.s. and allows to bring together the SoV representation and the coordinate representations of the various overlaps.
Whereas the proper proof of the proposal (6.10) requires further insights from the operatorial SoV approach, like in [7], we can analyse some common features of the l.h.s. of (6.10), which is represented by a 4J-dimensional integral of two CFT wave-functions, and the r.h.s. which is given by the determinant of one dimensional integrals of bilinears of Q-functions. One key property, discussed in section 2.4.1, is the logarithmic divergence in the case ∆ A = ∆ B . This can be reproduced in the r.h.s. of (6.10) as some of the Q-bilinear forms are also divergent when ∆ A = ∆ B , due to a non-regularisable 1/n-divergence in (4.9). In order to mimic the cutoff in the coordinate space, we can introduce a slight difference in the twist eigenvalues λ a , then the 1/n-divergence will get replaced in (4.9) by Analysing the finite part of the determinant M AB under such regularisation, one can also deduce the SoV representation for the particularly important case ∆ A = ∆ B . We will leave all these very intriguing questions for future investigation. In particular, this includes the question of how to pick correctly the 4 combinations of a, c out of 16 possibilities. Let us point out, however, that a similar ambiguity/freedom can be observed in the spin chains in finite dimensional representations. In this case, all Q-functions are (twisted) polynomials and multiple combinations of indices a and c would produce the correct expression for the norm in SoV. Different possibilities correspond to using different reference states in the algebraic Bethe ansatz.

On the structure of the g-function in separated variables
The g-function is an important object appearing in studies of integrable systems with boundaries [70], and recently connected to interesting observables in N =4 SYM [54,55,57,71]. In our context it can be interpreted as an overlap between a CFT wave-function and a fixed boundary state. Whereas the exact form of the overlap would depend on the details of the boundary state, it also contains a universal part, which is, however, very hard to calculate. This universal part satisfies special selection rules, and in this section we propose a construction based on the Q-bilinear forms which nontrivially obeys these properties and is a suitable candidate for the universal part of the g-function.
The construction here is inspired by [6,53,72], where it was observed that the universal part has a very suggestive structure, which is related in a simple way to the expression for the norm -both in the case of spin chains and for some field theories.

Review: the g-functions and their parity selection rules
We start with some introductory remarks where we will review the definition of these observables, and describe the parity selection rules that they have to satisfy.
The g-functions measure the overlap between a generic state |Ψ and an integrable boundary state B|. Usually the boundary state B| can be obtained as a "Wick rotation" of an integrable boundary condition, which satisfies a suitable boundary-Yang-Baxter equation [73]. For example, a version of the Wilson line in the fishnet theory was shown to be an example of an integrable boundary condition [74] and the corresponding boundary state can be written explicitly as (see figure 4) 1) where for simplicity we assume that the contour is in the (1, 2) plane. 20 Other examples from N =4 SYM could arise in various situations: for instance, when we consider the expectation value of a single-trace operator in presence of a domain-wall defect [55,56], or a Wilson loop [57], or its contraction with two determinant operators [54], [75]. In N =4 SYM, all these quantities have been related to a g-function with an appropriate boundary state.
A simple characterisation of integrable boundary states was proposed in [76] for integrable 1+1 dimensional QFT and generalised in [77] for lattice systems. This condition states that quasi-cyclic invariant integrable boundary states are annihilated by all the conserved charges of the spin chain that are odd under a certain "chain-reflection" operation Π. Assuming that Π also maps the integrals of motion into linear combinations of the integrals of motion, they can be organised in two families, denoted as H + and H − , which are respectively even or odd with respect to Π: Then the nontrivial integrability condition of boundary states proposed in [77] reads: namely the boundary states has to be annihilated by all odd charges. Usually (7.3) is a nontrivial condition, which we verified for |B W L explicitly at J = 2 to be true for trivial twists λ a = 1. In our case, a possible choice of the space-reflection is the transformation of a CFT wave function where in addition we perform a reflection in the space-timex α = (−x α,1 , x α,2 , −x α,3 , x α,4 ). Note that the boundary state B W L | is invariant under this transformation, as reordering changes the orientation of the WL, which is then restored with the spatial reflection. The g-function is defined as the normalised overlap between an integrable boundary state and an eigenstate of the integrals of motion |Ψ : Note that (7.3) implies that there is a selection rule for g-function namely, the overlap is non-vanishing only for symmetric states. The g-function was studied quite intensively recently. One can encode the information on the boundary state |B through the boundary reflection operator [78], which in many interesting cases reduces to a reflection phase Θ(u). The first proposal for a method to compute the g-function in integrable systems was made in [79], as a convolution involving the reflection phase and the Y-functions. However, it was first pointed out in [80] that the g-function must also contain an additional universal factor, that does not depend on the boundary state, but only on the state |Ψ . A correct proposal for this universal factor in rank-1 system was made for the first time in [58]. Later, path integral arguments to prove this result were presented in [81] (reproducing part of the result), and [82] (reproducing the whole result). Recently a rigorous derivation of the general form of the g-function, extended to generic rank, was given in [83], see also [54] for a different argument.
Schematically, the structure of the g-function is the following: where we are omitting sums over different types of Y-functions for simplicity. The universal factor is a ratio of Fredholm determinants of integral operatorsĜ ± , which are defined in terms of the Y-functions (see for example [54] for details).
Here, we want to discuss the structure of the universal factor in terms of the Qfunctions, following similar arguments to [53] in the case of the Sinh-Gordon model. In that work, it was suggested that the Fredholm determinants can be related to factors of definite parity of the determinant defining the norm of the state in separated variables. 21 An important guiding principle in formulating such a proposal is that any g-function has to satisfy the selection rule (7.6), which should be a property of the universal factor. Below we demonstrate that there is such combination of Q-functions in the fishnet theory, which we propose as a natural candidate for the universal factor.

Chain-reflection symmetry in the fishchain
We start by discussing the chain-reflection symmetry Π in the fishnet spin chain. For simplicity we consider the case with no magnons M = 0. In fact the presence of magnons explicitly violates parity, as in the fishnet theory the magnons spiral in one direction around the "worldsheet" of the Feynman diagram. The chain-reflection, acting on the CFT wavefunction, is defined in (7.4).
Now we need to see how this symmetry acts on the integrals of motion. In order to map the integrals of motion to each other we will additionally require that the twist eigenvalues satisfy λ 1 = 1/λ 2 , λ 3 = 1/λ 4 and that the inhomogeneities are such that θ α = −θ J+1−α . A reversal of the order of spin chain sites is equivalent to transposition in the auxiliary space in the definition of the transfer matrices. From the definition of the Lax matrices in appendix F, we see that (L 6 (u)) M N = (L 6 (−u)) N M , (L 4 (u)) a b = −(L4(−u)) a b , (7.8) namely, transposition in auxiliary space sends u ↔ −u and 4 ↔4. The reflection x α →x α is a conformal transformation which in the 4 representation interchanges 1 ↔ 2 and 3 ↔ 4, and thus transforms the diagonal twist matrix Λ 4 = diag (λ 1 , λ 2 , λ 3 , λ 4 ) → diag (λ 2 , λ 1 , λ 4 , λ 3 ) = Λ4. At the same time the Λ 6 is not diagonal but the transposition is exactly compensated by the twist interchange. Thus we conclude that which then results in the following transformation of the polynomials P r (u): As the parity changes the eigenvalues of the integrals of motions, it also acts nontrivially on the eigenstates and consequently on the Q-functions. If the initial q(u) was solving the TQrelation B • q(u) = 0, then the Baxter equation with transformed coefficientsB •q(u) = 0 is solved simply byq(u) = q(−u). More precisely, taking into account the conventions for the asymptotics (2.56) and analyticity we getq where σ is a permutation acting as {1, 2, 3, 4} → {2, 1, 4, 3}.
Parity-even and parity-odd integrals of motion. Finally, from (7.10) we introduce even and odd combinations of the integrals of motion: Parity-symmetric states are the ones for which H (−,α) = H (−,α) = 0, which implies in particular that they must have zero spins S i = 0, i = 1, 2, whereas ∆ is parity even and does not transform.

Expression in terms of the Q-bilinear forms
In section 6 we introduced a linear system of equations associated to any pair of states A, B (see (6.7)), the solution of which is the difference of their integrals of motion. Now let us apply the same argument for a pair of states chosen as follows: a generic state A, and the stateÃ which is the image of A under the parity transformation Π. Note that in this case the values of ∆ for both states coincide and we have to take a = c in (6.7) in order to keep the matrix elements M AÃ finite. These two states must have the same value for all the even charges H + , H + , while the odd charges differ by a sign In this case the null eigenvector I A − I B in (6.7) has only 2J nontrivial components, since when written in the basis of H ± the components corresponding to H A + − HÃ + should be zero by construction. After this change of basis of the integrals of motion, we get where the blocks m ± are constructed explicitly in terms of Q-bilinear forms in appendix H (in particular, see (H.6)). Now let us consider a non-symmetric state, whereÃ = A. In this case the system of equations (7.16) has a nonzero solution, which implies that any 2J × 2J minor extracted from the left half of the linear system has to vanish. In particular, for a non-symmetric state we have: The reason we pick this particular minor is because, as we show in appendix H, for the partity symmetric states A =Ã the matrixM AÃ becomes block diagonal: 19) and the only minor which is non-zero in this case is |M − |. For example, we computed these determinants in the case J = 2, M = 0, for twists which has |M + | 0, while its null vector is proportional to ∂ ξ 2 ( H + , H + ).
The structure of the g-function.
We see that the determinant |M − | vanishes if and only if the state A is parity-symmetric. This is the same selection rule (7.17) exhibited by the overlaps with integrable boundary states. That makes it a natural candidate for the universal part of the overlap appearing in the g-function, The meaning of this proposal is that we expect the r.h.s. to appear naturally in the SoV construction of the overlap, in the same spirit as the discussion of the scalar product in section 6. Crucially, the relations (7.17) imply that the r.h.s of (7.22) satisfies the selection rules, which provide an infinite number of consistency conditions. However, one would still have to figure out explicitly how to build the non-universal part of this overlap, in terms of the boundary reflection operator corresponding to a given boundary state.
Finally, let us point out that for the parity symmetric states, due to the block diagonal structure of (7.19), we also find that the determinant corresponding to the finite part of the norm has a factorised form: the same way as was found in [6,53,72]. As we discussed in sections 5 and 6, for the standard choice of Q-bilinear forms the determinant (7.23) vanishes, which in this case is due to |M + | = 0. To find a connection with the (finite part of) the norm of the state, a regularisation prescription is required, which would be clarified in the operatorial SoV construction. Here, we point out that a natural possibility is to consider |M + | * , defined as the product of the 2J − 1 non-vanishing eigenvalues of |M + |. This leads to a proposal for the structure of the universal factor of the g-functions, similar to [53]: which is also very similar to the structure found rigorously in [6] in the case of the Heisenberg spin chain.

N = 4 super Yang-Mills: generalisations and speculations
In the previous sections we demonstrated that the functional SoV approach to spin chains can be also applied in such integrable field theories as the fishnet CFT. Here we speculate on how this method can be further extended to more complicated cases such as N =4 SYM. In this section we assume the most general deformed version of the theory, which also includes the twist parameters λ a in AdS 5 . In this case, the theory is described by a more complicated P SU (2, 2|4) Q-system. Four of the Q-functions, usually denoted as Q i , are directly connected with the Q-functions of the fishnet model in the double scaling limit of [42]. They represent AdS 5 degrees of freedom in the dual string picture. In addition, in this case we also have the P a functions, which represent the S 5 degrees of freedom. In the fishnet limit they simplify drastically and become rational functions with a pole at the origin, and the coefficients in the numerator of these rational functions encode the integrals of motions I (b,α) .
The P and Q functions have a characteristic analytic structure illustrated in figure 5. The size of the cuts is given by the 't Hooft coupling g = √ λ 4π and tends to zero in the fishnet limit, with cuts of Q a becoming the already familiar poles of q a , whereas P a become rational functions with a pole at the origin.
The Baxter equation satisfied by q a generalises to a fourth order finite-difference equation [90]: where the coefficients are P a dependent combinations defined as 3) Figure 5. Analytic structure of P and Q functions. In the fishnet limit g → 0 and the cuts collapse into poles.
and D k are some combinations of P a 's defined in appendix I. In the fishnet model, the transfer matrix eigenvalues appearing in the TQ equations were polynomials, leading to a natural basis for the integrals of motion. The coefficients A k (u) have a more complicated analytic structure, see appendix I. Crucially, they are asymptotically constant at infinity, and are analytic outside a certain radius R * in the complex plane. This means that they admit a convergent Taylor expansion around infinity: The coefficients I (k,n) should contain a basis of integrals of motion, generalising I (b,α) of the fishnet theory. At any fixed order in the weak coupling expansion, only a finite number of integrals of motion are linearly independent: as it is known from the scaling of the coefficients of P a , at order O(g 2l ) such functions truncate to finite Laurent polynomials in u with ∼ l coefficients, see e.g. [84]. However, at finite coupling, we have an infinite number of integrals of motion, as expected since the theory corresponds to a sigma model with infinite number of degrees of freedom (unlike the fishnet theory, which is dual to a system of J particles with 4J degrees of freedom [47] 22 ). Finding the most convenient basis of the coefficients I (k,n) , which behaves nicely at weak coupling and takes into account the cut structure of the coefficients A k (u), is an important problem for future studies. One possibility is to use a spectral representation of A k (u), which has a finite number of cuts, with the discontinuities expanded into Zhukovsky variables x n -this basis will behave well at weak coupling, but could still contain some linear dependencies.
In order to advance further in the analogy with the fishnet theory, we have to introduce Q-bilinear forms and make sure that the Baxter equation (8.1), written as a finite difference has a simple conjugation property under this bilinear form. Fortunately, this step can be realised in N = 4 SYM too. We introduce the dual Q-function Q i by which in the QSC Q-system notation correspond to the Q-functions Q ∅|ijk . This dual Q-function does indeed satisfy the dual TQ-relation as we demonstrate in appendix C, where B dual is the conjugate operator to B in the same sense as in the fishnet CFT case, i.e. as defined in section 4.1.2. This means, for example, that we can repeat the derivation of the orthogonality relations from section 6, giving the SoV form of the scalar product starting from the equation where as before we define the bracket f µ as a formal integral over the where we assume that c > 2g. Such an integral will encircle the infinite ladder of cuts of the Q-functions. In practice, say for numerical calculations, one would write this integration as a sum over the infinite set of contours encircling the branch-cuts of Q i and Q j . Then for each cut one can use a strategy very similar to section 4.2 -one can flip the location of poles by using the i-periodic, anti-symmetric matrix ω ij (u), which then should allow us to have full analytic control over the tail of the sum over cuts and enable efficient ζ-regularisation just like in the fishnet case. A convenient set of periodic functions µ(u) is which generalises naturally the corresponding basis in the fishnet CFT (4.3). Furthermore, like in the fishnet case, there are exactly 4 combinations of the indices (i, j) in the Q-bilinear forms (8.8) which give a finite result in the important case when ∆ A = ∆ B and also the twists λ A a = λ B a . Exactly the same combinations of the Q-function indices (i, j) ∈ {(1, 2), (2, 1), (3,4), (4, 3)} lead to a finite result for the integral (8.8) when the Q-functions are on-shell i.e. all QSC analyticity constraints are satisfied including the gluing condition [87,88]. For the more general deformed setup of two states with different twists, we expect again that the integrals are convergent for generic pairings the same way as was discussed in section 6.
Of course the plan outlined in this section deserves very intensive further studies, but we believe the main ingredients presented here will become instrumental in the future. The main difference with the fishnet case is that the number of integrals of motion is now infinite, leading to infinite size determinant for the SoV scalar product. In order to overcome this additional complication, one should develop a good truncation at perturbative level, where the number of independent integrals of motions becomes effectively finite. Then this perturbative method can be extended to finite coupling by an appropriate extrapolation procedure, similarly to how that is done for the spectrum [66].
One of the applications of this approach would be the calculation of the overlaps between local operators and integrable boundary states, given by the g-function [53,54]. For that one should follow the steps done in section 7 for the fishnet model. The main message of this section is that for the applications of the SoV methods to models like N = 4 SYM one can essentially follow the same steps as for the fishnet model. This makes the development of the SoV for spin chains in general paramount for the progress with the non-perturbative description of correlators in AdS/CFT.

Summary
In this paper we initiated the program of studying the multiple wrapped diagrams in fishnet theory by means of Separation of Variables methods. This approach promises to extend the power of the Quantum Spectral Curve and alike methods, initially developed for the spectrum, to more general classes of observables. In contrast with other integrability based methods for correlation functions the presented approach resums all wrapping diagrams and the calculations are valid non-perturbatively for finite length operators 23 .
In this paper we develop the so-called functional SoV approach introduced in [15]. The next obvious step would be to join our results with the insights coming from the operatorial SoV approach [4,13,20], which would allow us to construct more general overlaps and form factors in fishnet theory like in [13]. Moreover, an explicit construction of the CFT wave functions (which should be obtainable with the operatorial approach) would be a very significant step towards the full solution of the planar theory. This is because the wave functions contain all the nontrivial coupling dependence of generic correlators.
Within the functional SoV approach, we made a proposal for the expression for the scalar product in the SoV basis. Furthermore, we derived a closed expression for the derivative of the dimension w.r.t. to the coupling constant ξ, which is written in terms of the Q-functions evaluated at one fixed value of ξ. This result, which is related to a 4J × 4J determinant of one dimensional integrals of Q-functions, gives a prototype for the type of structures which should appear in more general correlators. We also developed general numerical methods for efficient evaluation of such integrals of Q-functions. This method is based on the core properties of the Q-functions and the QQ-relations.
Another application of our method is the computation of the g-function, which is closely related to the scalar product from the SoV point of view [6,53]. We discuss this in detail in section 7.
One unexpected observation we made is that the quantisation condition, which comes as an extra condition on the Q-functions, is tightly related with the finiteness of the norm in SoV representation.
In section 8 we speculated on how our methods can be extended and applied for the more complicated theories like N = 4 SYM. There, we argued that all main ingredients of our construction find their counterparts in this more general case, even though more detailed future studies are required.
Let us also mention that in a particular (Gaudin) limit the spin chain we considered here should be linked closely to the computation of multi-point conformal blocks in any dimension, in the light of the results from [91]. It would be important to explore the implications of SoV techniques developed here in this context, where further simplification is expected.
Another important direction is to apply the methods of this paper to the boundary problems like [74]. This could lead to generalisation of the initial observation of [35] for the case of the simple cusp 3-point correlator.

A Sigma matrices
The Sigma matrices Σ M N in our conventions read

B Derivation of the form of the scalar product
The scalar product is a bilinear functional. Without loss of generality, we can represent it in the form whereK is a linear operator, which can be written in terms of an integration kernel (which might be a distribution): Let us now impose the covariance of the scalar product under a generic conformal transformation C −1 , acting on the site n as The covariance condition reads: where D is the spacetime dimension, i.e. D = 4 in our case. Since this must be satisfied for all sites n, we can look for a solution in terms of a factorised kernel: Imposing the conditions (B.4) for all sites n is equivalent to: This is the same transformation rule as for a two-point function with scaling weight (D−h). The solution is therefore fixed to: Multiplying by an appropriate normalisation factor, this integral kernel can be related to the fractional operator

C Dual Baxter equation
Here we describe in general the form of the dual Baxter equation satisfied by 3 × 3 determinants of Q-functions. Consider a finite difference equation of the form A −2 (u)Q a (u−2i)+A −1 (u)Q a (u−i)+A 0 (u)Q a (u)+A +1 (u)Q a (u+i)+A +2 (u)Q a (u+2i) = 0 (C.1) where a = 1, . . . , 4 labels 4 independent solutions of this equation. Define dual Q-functions by P a (u) = F (u) aa 1 a 2 a 3 Q a 1 (u + i)Q a 2 (u)Q a 3 (u − i) .

(C.2)
Then the 4 functions P a (u) satisfy the dual finite difference equation To show this we plug the definition (C.2) into (C.3) and then use (C.1) to exclude Q a (u + 3i), Q a (u + 2i) and Q a (u + i) from the resulting equation.
To make the result simpler, we can make use of the arbitrary overall multiplier B −2 (u) which we set to be B −2 (u) = A +2 (u − 2i), and also adjust the factor F (u) in the definition of P a (u) so that it satisfies . (C.9) Then we see that the dual Baxter equation (C.3) takes a very simple form,

D Additional properties of the Q-functions on-shell
In this appendix we review a reformulation of the quantisation introduced in [65], and already used in [49]. We will consider vanishing rapidities θ i → 0 in this section. The most important aspect of this alternative method, is that it provides a way to compute the eigenvalue of the graph-building operator, which also applies in the more complicated case of quantum numbers |M | = J. Consider where we assume that the Q-functions are on-shell, namely satisfy the quantisation condition, and Γ ab is the gluing matrix. We want to show that, as a consequence of the quantisation conditions, this function has an analyticity strip [−in/2, in/2]. To see this, we rewrite and using the antisymmetry of ω ab = (Γ · Ω) ab , we find where the r.h.s. is explicitly analytic at u = 0, . . . , −in. Furthermore one has Among these quantities, a particularly important one is Q + (u) ≡ T (u), which has a pole of order J at u ∼ 3i 2 . From the coefficient of this pole, one can extract the eigenvalue of the Hamiltonian (identified with the coupling constant) [65]: .

(D.5)
This equation applies to all states of the theory, including those with |M | = J for which the simple relation I (0,1) = ξ 2J is not valid.

E Covariance of twisted correlators
In this appendix we review the covariance transformation rules for colour-twisted correlation functions in a large-N CFT, and the insertion points x i are fixed points of these maps, see [49]. Consider a correlator of colour-twisted operators, where G 1 , . . . , G n are conformal transformations used in the definition of the twisted operators. Under any conformal transformation C, the correlator is covariant in the sense that wherex i ≡ C(x i ), and the twist maps are also conjugated asG i ≡ CG i C −1 . We considered scalar operators for simplicity.

Two-point function at opposite twists.
A particularly important case, already studied in [49], is the 2-point function of two operators with opposite twists, G 2 = G −1 1 . In this case, the 2-point function has the same kinematical dependence on the points as in the standard CFT without twists. Let us recall how this works. We consider twist maps with fixed points x 0 , x0. Such maps always admit the following decomposition (we use the notations explained in the main text): where Λ is the diagonal part, r is a 4D rotation, and K is the special conformal transformation defined in (2.26), satisfying K(0) = x 0 , K(∞) = x0.
It is simple to show that the 2-point function O G (x 0 )O G −1 (x0) does not depend on the form of the rotation r. In fact, consider the covariance rule (E.1), for the map C = K • r • K −1 , where r is any different SO(4) rotation. Since such map C has unit Jacobian at the fixed points, we conclude that the 2-point function is the same with r → r • r.
Therefore, the only relevant parameters are the fixed points and the diagonal part of G: we can write O G (x 0 )O G −1 (x0) ≡ F(x 0 , x0; λ). Now, (E.1) just becomes the familiar so we can proceed as in the standard CFT case. Considering the effect of dilatations, rotations around an arbitrary point, and translations, we conclude that Finally, consider the covariant transformation property for the map C ≡ K • e ρD • K −1 , which leaves G invariant. Then, (E.1) implies that the 2-point function has to vanish unless Transformation of the CFT wave function.
For the CFT wave function, which is defined in terms of the twisted operator O G and a twisted trace (see (3.3)), the covariance property (B.4) gives where ∆ is the scaling dimension of the operator O. Notice that we took into account that, as a correlator, the wave function depends on J + 1 insertion points, one of them being the fixed point of the twist x 0 where the operator O sits.

F Properties of the transfer matrices under transposition
In this appendix we derive some properties of the adjoints (or transposes) of the transfer matrices with respect to the conformally invariant scalar product, which are denoted as We will derive a convenient explicit equation for the transpose, as a combination of reversal of the order of sites in the chain, and conjugation with a special conformal map. This representation gives a simple map between right and left eigenvectors in the case of homogeneous chain, and allows us to prove directly that the transfer matrices and their transposes have the same spectrum. At the end of this appendix, we also use this relation to show that, in the fishnet theory, the left eigenvectors can be interpreted as conjugate CFT wave functions.
In this section we will obtain the formula: whereT r rev denotes the transfer matrix with reversed order of sites inside the trace: and F = K •Ĩ • K −1 is the conformal map defined in (2.47), withĨ the holomorphic inversion (2.48). 24 For clarity, in this appendix we occasionally denote the Lax matrices and conformal generators aŝ L r x k ,h k ,q M N x,h , to mark the space-time variable on which they act, and the weight of the representation.
Properties of F . Before presenting the proof, let us list some important properties of the conformal map F , which can be easily verified.
-It is equal to its inverse, F = F −1 ; -In the frame reached by this transformation, the twist is inverted: where we usedĨ • Λ •Ĩ = Λ −1 for the diagonal part of the twist; -The map changes the sign of the Cartan generators: -The map swaps the two fixed points of the twist map x 0 = F (x0), x0 = F (x 0 ).

F.1.1 Derivation
Transposition of Lax matrices. From the very definition of the conformally invariant scalar product, the adjoints of the conformal generators satisfy: which is valid locally at every site. In virtue of this relation, the Lax operators satisfy where we kept explicit the indices in auxiliary space, and we remind that O T denotes the transpose of an operator w.r.t. the scalar product, i.e. the adjoint.
Intertwining with F . There are two key identities that we will use in the following: where we kept explicit the representation indices 1 ≤ a, b ≤ 4, and the holomorphic inversion acts in 6D asĨ 6 = diag{1, −1, 1, −1, −1, 1}. Combined with (F.8),(F.9) , these identities give Notice that the terms in round brackets on the right hand side of (F.12),(F.13) involve a transposition of the auxiliary space indices.
Transfer matrices. We will work in the frame where the twist map is the diagonal transformation G ≡ Λ. Consider the case of the vector representation, defined aŝ T 6 = Tr 6 L 6 X J (u − ϑ J ) . . . L 6 X 1 (u − ϑ 1 ) Λ 6 , (F.14) where the Lax matrices without explicit mention of the indices are assumed to have the index structure (L 6 ) M N ≡ (L 6 ) M M η M N . Using (F.8), for its transpose we find (T 6 ) T = Tr 6 L 6 which is a similar to the definition of the transfer matrix, but with the order of sites reversed, and a different twist. Above, by A Taux we denote transposition in the auxiliary space indices. The twist matrix appearing in (F.15) can be rewritten as η · Λ 6 · η −1 Taux = (Ĩ 6 ) Taux · Λ 6 ·Ĩ 6 . (F. 16) Using the fact that (Ĩ 6 ·Ĩ 6 ) N M = δ N M andĨ 6 = (Ĩ 6 ) Taux , we can use (F.16) to write where we used (F.10) in the last line. This matches the announced result (F.3). Now we considerT 4 . Identity (F.13) tells us that, upon integration by parts, L 4 (u) gets transposed in auxiliary space, and conjugated with the holomorphic inversion in quantum space. The twist matrix Λ 4 is not affected by such transposition, since it is diagonal in this representation. Therefore, repeating the steps described for the previous case, we find immediately The derivation for the case of4 is a simple generalisation. We have presented the above results (F.17),(F.18) for the cases where the twist matrix is G ≡ Λ. The general case is related to this by a conformal transformation denoted as K in the main text. Using the covariance of twisted transfer matrices (2.31) under this transformation, we see that this has simply the effect of replacing the twist Λ → G, as well as changing the holomorphic inversion to the map F . This leads to the announced result (F.3).

F.2 Spectral properties
Now we establish that the spectrum ofT r on the subspace of wave functions with fixed conformal charges Q R a = {i∆, S 1 , S 2 }, is the same as the spectrum of (T r ) T , on the subspace of wave functions with opposite conformal charges Q L a = − {i∆, S 1 , S 2 }. We have seen that the adjoint of the transfer matrix is given by the chain of Lax matrices applied on the sites in reverse order, and conjugated by the map F defined in (2.47): (T r ) T = F •T r rev • F, (F. 19) with T r rev ≡ Tr r L r x 1 ,h 1 (u − ϑ 1 ) · · · · ·L r x J ,h J (u − ϑ J ) · G . (F.20) First, we notice that the Baxter equation is invariant under synchronised permutation of the inhomogeneities and the order of magnons. Therefore, the spectrum is invariant under this operation. This means that, for every right eigenstate satisfyingT r • Ψ R A = T rR A Ψ R A , there is an associated wave function Ψ A,rev satisfyingT r rev • Ψ A,rev = T rR A Ψ A,rev , with the same eigenvalue. For fixed conformal charges, the map between states Ψ R A and Ψ A,rev is one-to-one since one expects the spectrum is non-degenerate. This wave function has the same Cartan charges Q R a as Ψ R A , since they can be read from the integrals of motion and are unaffected by permutations of the sites. Now, we can construct Ψ L A ≡ F • Ψ A,rev , which, by the identity (F.19), will be an eigenstate of (T r ) T with the same eigenvalue. Moreover, if Ψ R A (and therefore Ψ A,rev ) has Cartan charges Q R a , due to property (F.6) we find that the left eigenvector Ψ L A has exactly opposite Cartan charges.

F.3 Left eigenvectors in the fishnet case
In this section we consider the fishnet theory, specifying to zero inhomogeneities.
As we discussed in section 3.1, the CFT wave functions, which are right eigenvectors of the integrals of motion, are defined as certain correlators: Notice that going from right to left eigenvectors has the effect of inverting the twist of the operator, and changes the insertion point from x 0 to the second fixed point x0 = F (x 0 ). Using the discrete symmetry of the fishnet model under conjugation of the elementary fields, φ i ↔ φ † i , the left eigenfunction can also be written as This relation is useful to derive the interpretation of the norm in CFT, presented in section 3.3. In particular, in Figure 3, the sum of Feynman diagrams on the left of the green cut can be identified with the correlator on the r.h.s. of (F.26); while the diagrams on the right side of the cut give the CFT wave function (F.21).
G Details on the ∂ ξ ∆ calculation and relation between scalar product and 2-pt function In this appendix we consider the eigenfunctions for zero inhomogeneities, corresponding to CFT wave functions of the fishnet theory. We show how they are normalised according to the scalar product. This argument was already implicit in [46] and we thank A. Sever for discussing it with us. We report it here for completeness. For simplicity, we will consider the state for lowest scaling dimension at any J, with M = 0.
We will start from the perturbative theory, and then resum the diagrams. We take the propagator to be where W (n) represents the fishnet diagram with n wheels 25 26 : (G.3) In the coincident points limit, (x i − x 0 ) 2 ∼ (y i − y0) 2 ∼ 2 , this quantity will be denoted as W , (x 0 |x0). It exhibits a power-like divergence, which we can remove by multiplicative 25 We extract a prefactor with the coupling as B already contains it 26 The J vertices in each wheel contribute a (16π 2 ) J factor and the propagators linking them give 1/(4π 2 ) J , leaving 4 J . It is partially cancelled by the corresponding ladder (spokes) propagators that give 1/(4π 2 ) J which leaves only the 1/(π 2 ) J factor that we have included in B. Looking at the explicit expression for the scalar product, we notice that the factors first remove one layer of propagators, and after that the ladders are glued together. This produces some combinatorial factors, so that the result is given by the sum of diagrams: a = ∞ n=0 (n + 1) ξ 2Jn W (n) , (x 0 |x0). (G.9) Notice that this has the form of a derivative, in fact we can obtain the same sum as a = 1 J ξ 2 ∂ ξ 2 + 1 (W , (x 0 |x0)) ∼ 2 J (ξ 2 ∂ ξ 2 ∆) log x 00 2∆−2J N x 2∆ 00 + . . . , (G. 10) where the second term contains subleading singularities (without the log for example). Comparing with (G.8), we conclude that which is the result stated in (3.18).

H Parity structure of integral orthogonality relations
In this appendix we follow the approach of section 6 to derive a linear system of equations associated to a generic state A, and its Π-transformed state denoted byÃ, where Π is the chain-reflection symmetry discussed in section 7.
With similar arguments we find in general, in analogy with the case studied in [53].
Comment on the selection rule. As proved in the main text, |M − | = 0 whenever the state is non-symmetric, A =Ã. Let us show that, instead, when A =Ã, we have |M − | = 0 while |M + | = 0. This is a consequence of the variation equation (5.14). Since for a symmetric state the odd integrals of motion have the vanish, the latter reduces to a square system of 2J equations M + · ∂ ξ 2 H + H + = 0, (H.14) which implies (since the system has a nontrivial solution) that |M + | = 0. (H. 15) As discussed in section 5, the kernel ofM is expected to be one-dimensional in generic points of parameter space. This, together with (H.15), implies that, for symmetric states, |M − | = 0. This determinant therefore gives a precise characterisation of the symmetry of the states: it vanishes precisely for states which are not symmetric.

I Determinants in the N=4 SYM Baxter equation
Here we present the determinants appearing in the coefficients of the N = 4 SYM Baxter equation we presented in (8.2): (I.5) Since the P functions have a single short cut on the first Riemann sheet (see figure 5), these determinants all have a finite number of cuts, and are analytic outside a certain radius in the complex plane. We also notice that, using the QQ relations, these coefficients can be written as the same expressions, but where all the P functions are replaced by Q functions with the same index. The expression in terms of Q functions makes it easy to evaluate the determinants for large u when the Q functions have generic twists λ i . In this case, the determinants are non-vanishing constants at infinity. These analytic properties justify the expansion (8.4). The fact that we can alternatively write the D coefficients in terms of the Q or P functions, implies that they are functions with a finite number of either short or long cuts, on two connected Riemann sections. This suggests that they might have a simple analytic structure. Understanding it more fully might be relevant to find a good basis of independent integrals of motion in N =4 SYM.