Quantum mechanics of null polygonal Wilson loops

Scattering amplitudes in maximally supersymmetric gauge theory are dual to super-Wilson loops on null polygonal contours. The operator product expansion for the latter revealed that their dynamics is governed by the evolution of multiparticle GKP excitations. They were shown to emerge from the spectral problem of an underlying open spin chain. In this work we solve this model with the help of the Baxter Q-operator and Sklyanin's Separation of Variables methods. We provide an explicit construction for eigenfunctions and eigenvalues of GKP excitations. We demonstrate how the former define the so-called multiparticle hexagon transitions in super-Wison loops and prove their factorized form suggested earlier.


Motivation
Recently, it was suggested that null polygonal super-Wilson loop may serve as a generating function for all scattering amplitudes in (regularized) planar maximally supersymmetric SU (N ) Yang-Mills theory. The equivalence was first established for the so-called Maximally Helicity Violating (MHV) gluon amplitudes, that involve two particles with helicities opposite to the rest, and the holonomy of a gauge connection circulating around a closed contour composed of straight light-like segments [1,2,3]. These segments are identified with momenta of incoming gluons involved in the scattering process and one regards the resulting loop as living in a dual coordinates space [4,2]. The conventional soft and collinear divergences of massless gluon amplitudes emerge in this picture from the ultraviolet singularities stemming from virtual corrections in the vicinity of the cusps where the two adjacent segments meet in a point. However, not only the divergent parts were shown to agree, but also finite contributions to both coincide as well. By uplifting the contour to superspace, the duality establishes the equivalence 1 between the vacuum expectation value of the super-Wilson loop W n in N = 4 SYM and the reduced N -particle matrix element A n of the S-matrix of the theory [5,6] that is expected to hold nonperturbatively in 't Hooft coupling a = g 2 YM N c /(4π 2 ). The former is determined by bosonic A αα and fermionic F αA connections that live on a piece-wise contour in superspace C n = [X 1 , X 2 ] ∪ [X 2 , X 3 ] ∪ · · · ∪ [X n , X 1 ] with cusps located at superpoints X i = (x i , θ i ). Both superfields A αα and F αA receive a terminating series in the Grassmann variable θ, with their expansion coefficients being determined by the field content of the N = 4 SYM (see, e.g., Ref. [7] for details). The dynamical content of the duality is summarized by the following equation W n (X i ; a) = A n (Z i ; a) , (1.2) where the right-hand side corresponds to the reduced n-particle S-matrix element upon extraction of the conservation laws for the supermomentum (p αα , q A α ) S n = i(2π) 4 δ (4) (p αα )δ (8) (q A α ) 12 23 . . . n − 1n A n (Z i ; a) . (1. 3) The kinematical variables entering the super-Wilson loop are related to the momentum supertwistors Z A i = (λ α i , µ iα , χ A i ) defining the superscattering amplitude via the relation 2 In addition to the conservation-law delta functions, we traditionally factored out the denominator of the Parke-Taylor tree amplitude [10], where the angle brackets between bosonic twistors are conventionally defined as ij = λ α i λ jα . A brute force calculation of the super-Wilson loop represents a challenge on its own even at lowest orders in coupling a, thus exposing inefficiency of usual techniques based of Feynman diagrams. It was noticed however that the left-hand side of the above relation (1.2) is advantageous in computing the S-matrix since the Wilson loop formalism allows one to break free from limitations of perturbation theory, being more amenable to nonperturbative techniques [1]. Such a techniques was put forward a little while ago and is based on the operator product expansion (OPE) [11]. Recently, a complementary view on this framework was offered by exhibiting its connection to the expansion in terms of light-cone operators with boundary fields set by light-like Wilson lines, or the Π-shaped Wilson loop with field insertions [12,13]. The renormalization group evolution of the latter was shown to be governed by a noncompact open spin chain [12], whose eigenenergies provide corrections to leading discontinuities of the Wilson loop, while the corresponding eigenfunctions define the coupling of the Gubser-Klebanov-Polyakov (GKP) excitations [14,15] to the Wilson loop contour as well as their multiparticle transition amplitudes. In this paper we explore the emerging open spin chain. We will focus on the component formalism leaving a fully supersymmetric superspace formulation for future publication. Notice that quite recently a framework that combines the OPE and integrability was proposed in Ref. [16]. It is based on a few conjectured equations obeyed by building blocks of the Wilson loops that are matched to explicit calculations of amplitudes and bootstrapped to all orders in 't Hooft coupling. The next section follows a similar construction.

Soft-collinear expansion of polygons
To start with, let us recall the emergence of non-local light-cone operators in the OPE of the super-Wilson loop. We will discuss it in a simplified set-up of restricted kinematics [17], i.e., when the bosonic part of the Wilson loop contour is embedded in a two-dimensional plane 3 . The first nontrivial Wilson loop in this kinematics is the octagon. Let us focus on NMHV amplitudes and choose a specific Grassmann component that admits a clear OPE interpretation, say χ 2 χ 3 χ 6 χ 7 . It was calculated diagrammatically in Ref. [18] and also making use of the so-calledQ-equation in Ref. [19]. The subtracted superloop, with ultraviolet divergences eliminated in a superconformally invariant fashion, reads to two loop accuracy W R 8 = χ 2 χ 3 χ 5 χ 6 2367 a + a 2 ln u + (1 + u + ) ln u − (1 + u − ) − 2 ln(1 + u + ) ln(1 + u − ) + O(a 3 ) , (1.5) where the conformal cross ratios are (1. 6) The leading contribution in the OPE of this channel stems from the scalar exchange between the cusps at points x 3 and x 4 , see Fig. 1 (a). The soft-collinear expansion corresponds to the limit x + 32 → 0 and x + 67 → 0, which results in flattening of the octagon to a square, see Appendix A.2. The subleading perturbative corrections come from a number of Feynman graphs encoding the interaction of the scalar field with the contour, Fig. 1 (b), vertex corrections and independent propagation of a gauge field along with the original scalar, Fig. 1 (c). A natural set of the cusp coordinates, or equivalently twistors, to discuss the OPE is introduced in Appendix A. 3. In this frame, the remainder function admits the form W R 8 = χ 2 χ 3 χ 5 χ 6 4 cosh τ cosh σ a + a 2 τ ln(1 + e 2σ )(1 + e −2σ ) + ln(1 + e −2τ )(1 + e 2σ ) + O(a 3 ) .
The afore-introduced soft-collinear asymptotics translates into the τ → ∞ limit, with the remainder function corresponding to a single scalar GKP excitation (also known as a hole) propagating in the OPE channel 4 Here the energy and momentum of the hole excitation read to one-loop accuracy [14,15]  .
As above, the leading asymptotic is governed by a single hole excitation where the hexagon transition kernel reads to this order in coupling H h (u|v; a) = H h (u|v) a −1 + 2ψ ( 1 2 + iu) + 2ψ ( 1 2 − iv) − π 2 + O(a) , (1.13) 4 Notice that our definition of the rapidity variable differs by a sign from the conventions of Ref. [16]. Figure 1: Single (a,b) and two-particle (c) contributions to OPE of the octagon. The one-loop graph in (a) given by the scalar propagator exchanged between the cusps produces a components of the tree NHMV amplitude. The graph in panel (b) displays one of the perturbative corrections due to the Hamiltonian acting on the light-cone operator. In (c) we show the two-particle contribution that produces subleading effects in the OPE.
with the leading term being . (1.14) As we recall in the next section, these leading contributions in the asymptotic expansion, i.e., O(e −τ 1,2 ), come from the renormalization of a nonlocal operator with a Π-shaped contour, formed by two Wilson lines building up the contour in a given OPE channel and an elementary scalar field insertion, while the subleading effects correspond to more than one GKP excitations sandwiched between the boundary Wilson lines.

Light-cone operators
Let us cast the above heuristic discussion into an operator language. We will do it for the octagon at two-loop order and then generalize it to an arbitrary number of GKP excitations exchanged in a given OPE channel. As one takes the limit τ → ∞, the top and bottom portions of the contour of the octagon Wilson loop get flatten out and one ends up with a scalar field Z inserted into the top and bottom straight-line segments. Each of these correspond to a Π-shaped light-cone operator of the form where the boundary fields W (x + ) stand for the Wilson lines stretched along the light-like direction x − starting at 0 and going to the null infinity 1/ε → ∞, localized at the position x + on the tangent light-cone, By a suitable gauge choice, i.e., the light-cone gauge A − = 0, one can even eliminate the gauge links between the fields in the composite operator, i.e., [. . . ] − → 1. The contribution of (1.15) to the the leading e −τ term in the OPE reads For a = 0, this can be easily identified as a propagation of a free scalar from the bottom to the top of the loop, while the O(a) effect induces the τ -enhanced contribution to the two-loop remainder function and arises from the interaction of the scalar field with the Wilson loop contour by means of the light-cone Hamiltonian  20) respectively. Here we presented them in a generic form suitable for studies of elementary fields of conformal spin s. While for the case at hand, s = 1/2. For a one-particle GKP excitation, the diagonalization of the above Hamiltonian can be easily performed with a plane-wave eigenfunction that carries the momentum p(u 1 ) = 2u 1 and the energy cf. Eq. (1.8). Analogously, we can perform the soft-collinear expansion for the decagon in terms of the light-cone operators (1.15). The formalism can be extended to account for multiparticle effects as well. As one collapses the top and bottom of the loop in the soft limit it gives rise to curvature corrections, actually an infinite series of them in increasing powers of the gluon field strength F +− . The first correction to the asymptotic behavior stems from the diagram shown in Fig. 1 (c) for the octagon. At even higher orders, the OPE involves multiparticle light-cone operators where we stripped the + superscripts off the x-coordinates for brevity since it is the only component that will appear in the analysis that follows. Under renormalization group evolution, these operators mix and one has to solve the resulting eigensystem problem. In this work we will consider a simplified setup when all GKP excitations are of the same type X n = X and possess a generic conformal spin s. The Hamiltonian built-up from nearest-neighbor interactions only, in planar limit, is the one-loop dilatation operator acting on the operator (1.23) Here the pair-wise Hamiltonians (including the boundary ones) can be encoded in a single formula Their complementary representation in terms of differential operators reads where ψ(x) = d ln Γ(x)/dx is the Euler digamma function. It is important to notice that this is not the form naturally coming from the computation of Feynman graphs contributing to the renormalization of (1.23). Namely, the one-loop light-cone Hamiltonians for interaction of GKP excitations among themselves are not sensitive to the boundary Wilson lines and thus have to enjoy full SL(2, R) invariance [20,21]. However, it is not the case in the above form due to the presence of the logarithms. The latter emerge from 1/α-factor in the second integral in Eq. (1.26) that was introduced in order to match it to the boundary Hamiltonians (1.19), with all intermediate logarithms canceling telescopically with one another. Therefore, in order to restore the conformal symmetry property, one can shift the logarithms into the boundary terms, such that the reshuffled Hamiltonians will become Here we used identity (B.6) from Appendix B to recast the boundary Hamiltonians in terms of logarithms.
At this point we would like to point out that similar Hamiltonians emerged in the discussion of the Regge limit of scattering amplitudes in Ref. [22].
Thus, the goal of this paper is to solve the spectral problem for the multiparticle Hamiltonian (1.24) (1.29) The problem turns out to be integrable since the Hamiltonian H N possesses enough integrals of motion with eigenvalues u = (u 1 , . . . , u N ) to be exactly solvable. Thus we have the powerful machinery of integrable spin chain at our disposal for constructing the explicit eigenfunctions and finding the corresponding eigenvalues. For the case at hand we are dealing with an open spin chain that does not possess SL(2, R) symmetry due to boundary interactions. In addition, as we can conclude from the one-particle eigenfunction (1.21), the system is not endowed with a vacuum states. As a consequence of this last fact, the traditional methods based on the Algebraic Bethe Ansatz [23,24] are not applicable and instead one has to rely on the method of the Baxter Q-operator [25] and the Separation of Variables (SoV) [26].
Our subsequent consideration is organized as follows. After introducing a natural Hilbert space for the Hamiltonian and a scalar product that makes it self-adjoint in the next section, we find a complete set of integrals motion commuting with H N . As a next step on the way to solve the spectral problem, we recall the factorization of R-matrices [23,24] into intertwiners that play a crucial role in our construction. In Section 3, we use them to construct the Hamiltonians and demonstrate their equivalence with the ones arising in the problem of renormalization of the light-cone operators (1.23). Then in Section 4, we construct the Baxter Q-operator as well as its conjugate making use of available techniques developed for SL(2) invariant magnets. Both operators are then used to generate commuting Hamiltonians. In Section 5, we devise an iterative procedure to construct the eigenfunctions of the Hamiltonian. We then cast the latter in a form of a contour integral representation as well as an integral representation in the upper half-plane. The latter is instrumental in demonstration of their orthogonality under the scalar product introduced earlier. To adopt our construction for multiparticle contributions to the null polygonal Wilson loops, we define a new scalar product on the real half-line in Section 6 and establish its relation to the one for the analytically continued functions in the upper half-plane. Then we explicitly compute the square and hexagon transitions in Section 7, where we also prove the factorized ansatz for the latter suggested in Ref. [16]. Finally, we conclude by pointing applications of the current construction to the operator product expansion of polygonal Wilson loops and outline the generalization to include GKP excitations of different spins. Several appendices contain conventions and computational details that we found inappropriate to include in the main body of the paper as well as a complementary construction the wave function in the representation of Separated Variables and its relation to the eigenvalues of the Baxter Q-operator.

Open spin chain
The Hamiltonian (1.24) defines a non-periodic one-dimensional lattice model of interacting spins S n = (S 0 n , S + n , S − n ) acting at each site defined by the coordinate of the GKP excitation. They form an infinite-dimensional representation of the sl(2, R) algebra, Our choice of the representation S ±,0 n for the spin generators S ±,0 on the fields X is driven by the form of the pairwise Hamiltonians (1.28) acting on the elementary fields X(x), i.e. [S ±,0 , X(x n )] = iS ±,0 n X(x n ) [20,21]. They are taken in the form of differential operators acting on lattice sites where all labels of the representation s are the same for any n. The latter condition defines a homogeneous open spin chain. The spin variable s is assumed to be real and bounded from below by 1/2.

Scalar product
It is convenient to define a scalar product on the space of functions Ψ u (x 1 , . . . , x N ). Its choice is a matter of convenience but it has to be adapted to the physical problem under consideration. As was explained in Section 1, we are interested in the τ -dependence of correlation functions of the light-cone operators O N which read schematically Here the functions Ψ u are determined by the operator matrix elements between the eigenstate |E u of Hamiltonian and the vacuum 5 , Ψ u (x 1 , . . . , x N ) = E u |O N (x 1 , . . . , x N )|0 , and the sum runs over all eigenstates. Although the variables x n are the light-cone coordinates of the fields X entering the Lagrangian of the theory and are, therefore, real, nevertheless distribution amplitudes Ψ u (z 1 , . . . , z N ) can be regarded as functions of z n in complex plane. It is known that they are analytic functions in upper half-plane and vanish at infinity. As a consequence of this change, the sl(2, R) generators in Eq. (2.2) depend and act on the z n variables accordingly.
In what follows we will consider the spectral problem (1.29) on the space of functions of N complex variables z n analytic in the upper half-plane and endowed with a standard SL(2, R) invariant scalar product [27]. Such a formulation appeared to be very useful for constructing the SoV representation for the closed and open SL(2, R) spin chains [28,29]. The scalar product on the space of functions holomorphic in the upper half-plane 6 is defined as follows [27] where z n = x n + iy n . The integration measure reads Dz n = 2s − 1 π dx n dy n (2y n ) 2s−2 θ(y n ) (2.5) and the integration runs over the upper half-plane due to the presence of the step-function θ(y n ). Under this scalar product the generators (2.2) are anti-self-adjoint with hermitian conjugation defined conventionally as The Hilbert space of the model is given by the tensor product of Hilbert spaces at each site ⊗ N n=1 V n , such that the generalization of the scalar product (2.4) to multivariable functions is straightforward, Dz k (Φ(z 1 , . . . , z N )) * Ψ(z 1 , . . . , z N ) .
(2.8) Here J n,n+1 is an operator related to the two-particle Casimir via the formula J n,n+1 (J n,n+1 − 1) = (S n + S n+1 ) 2 . (2.11) These are explicitly self-adjoint. In the following section, we show that there exists a complete set of integrals of motion (commuting with the Hamiltonian) that are hermitian with respect to this scalar product and thus their eigenstates form an orthogonal set. Closing this section we remark that to any operator A acting on the Hilbert space (2.8) we can associate a function A of N holomorphic and N anti-holomorphic variables, i.e., the so-called integral kernel, in a unique way via the relation (2.12) It turns out that integral kernels of operators which we construct in the following sections take a form of two-dimensional Feynman diagrams. All operator identities can be cast into the identities between Feynman diagrams which can be verified with the help of a standard diagram technique that we remind in Appendix B.

Integrals of motion
We start the construction of the integrals of motion following the conventional R-matrix approach [23,24]. The Lax operator acting on the direct product C 2 ⊗ V n of an auxiliary two-dimensional space C 2 and the quantum space at the n-th site V n is defined as with a complex spectral parameter u. We also displayed its dependence on the spin parameter s. The product of N copies of this operator in the auxiliary space determines the monodromy matrix T (s) (u), with its elements acting on the quantum space of the chain ⊗ N n=1 V n . As can be easily established from the Yang-Baxter equation with a rational R-matrix [23,24], each element of the monodromy matrix commutes with itself for arbitrary spectral parameters, i.e., [D N (u), D N (v)] = 0 etc. As we will demonstrate below, the operator D N (u) commutes also with the Hamiltonian (1.24) [D N (u), H N ] = 0 .
(2.15) It can be easily seen from the form of the Lax operator that the entries of the monodromy matrix are polynomials of degree N in the spectral parameter u with operator valued coefficients. In particular, Here the operator u n are the so-called Sklyanin's operator zeroes [26], defined by the relation As follows from Eq. (2.13), the integrals of motion d n are polynomials in the spin generators S n (n = 1, . . . , N ), for instance, with the rest d n being given by n-th order differential operators. Thus, Eq. (2.17) defines a system of N differential equations that once solved yields the eigenfunctions of the open spin chain. By virtue of the property (2.6), the operator D N is self-conjugate (for real u) and thus d n (and as a consequence u n ) are hermitian as well. This implies that all eigenvalues u n are real.

Factorized R-matrices and Hamiltonians
Let us now demonstrate that local Hamiltonians, arising from the R-matrices acting on the direct product of two copies of the quantum space V ⊗ V, indeed coincide with the one coming from diagrammatic calculations alluded to in Section 1.2.

Intertwiners
To start with let us recall relevant facts about the aforementioned quantum R-matrices. Let us represent the R−matrix in the form R 12 = Π 12Ř12 , where Π 12 is a permutation operator on the product of two spaces, Π 12 Ψ(z 1 , z 2 ) = Ψ(z 2 , z 1 ). For the operatorŘ 12 the conventional RLL−relation [23,24] takes the following form Here L 1(2) are the differential operators in z 1 (z 2 ), according to their expression (2.13) continued in the upper half-plane. In the above equation, we have clearly displayed all parameters which these operators depend on (notice that we will not assume in this section that s 1 = s 2 ). Equation (3.1) shows that the operatorŘ 12 interchanges the parameters of the Lax operators acting on the first and the second spaces, (u, s 1 ) ↔ (v, s 2 ). Making use of the explicit form of the generators it is possible to demonstrate that the Lax operator can be factorized into a product of triangular matrices As it is obvious from this representation, it is convenient to introduce the following combinations of the spectral parameters u, v and the spins s 1 , s 2 A fruitful approach to solving Eq. (3.1) was pioneered in Ref. [32]. It was suggested there to look for the solution of Eq. (3.1) in the form of a product of two operatorš that exchange corresponding spectral parameters (3.3) independently, i.e., These operators change the spin of the sl(2, R) representations and map the Hilbert spaces of the spin chain sites as follows where we explicitly displayed the parameters labelling the representation. The solution for intertwiners was found in Ref. [32] and they read where z 12 = z 1 −z 2 . These can also be easily cast in the integral form making use of the definition of the Euler Beta function. Note also that the operators R ± 12 depend only on the difference of spectral parameters, i.e.

Hamiltonians
The local SL(2, R) invariant Hamiltonians arise as coefficients of the leading power in the expansion of the R-matrix in Taylor series in the vicinity of u − v = ε → 0 [31]. The representation in terms of the Euler Gamma functions (3.7a) and (3.7b) then yields, Comparing their sum h 12 = h + 12 +h − 12 +2(ψ(1)−ψ(2s)) with the intermediate Hamiltonian (1.28), we immediately find that they coincide up to an additive constant. The commutation relation of these Hamiltonians with the monodromy matrix (2.14) can easily be established from the Yang-Baxter equations (3.5a) and (3.5b) for the intertwining operators by expanding both sides For N -site Hamiltonians, the commutation relations with the monodromy matrix can be established from the above equations (3.12) to be (3.14) There are leftover boundary terms in the above equations which have to be cancelled against some boundary Hamiltonians h 01 and h N ∞ to enforce the commutativity. We are going to find their form next. To project out the D N (u) entry from the monodromy matrix, one sandwiches the above commutation relations between the two-dimensional vector |↓ = (0, 1) T and its transposed ↓| = (0, 1). These equations in turn can be reduced to a set of local equations for the boundary Hamiltonians, namely, They receive the following solutions Summing these expressions up, we recognize the boundary Hamiltonians (1.28) that emerged from the gauge theory analysis (again up to an additive constant). Thus, finding the eigenfunctions of the Hamiltonian is equivalent to the diagonalization of the operator D N (u), which we address in Section 5 below.

Baxter operators and Baxter equations
Having in mind a goal of the present study to construct the eigenfunctions for the spin chain in SoV [26], we find in this section the Baxter operators Q ± . The latter form a conjugate pair with respect to the scalar product (2.8). We start our consideration by establishing their explicit expression as operators commuting with the element D N of the monodromy matrix. Then we devise finite difference equations that they obey and, finally, we find the local Hamiltonians stemming from them.

Baxter operator Q +
To start with, let us derive a defining relation that will play a distinguished role in the analyses which follow. Recalling that the intertwiner R + exchanges u + and v + spectral parameters, we can immediately find that for a product of Lax operators. Namely, the string of the R + -operators 7 moves the spectral parameter v + from left to right in the product of the Lax operators This can easily be proved by using the defining relation (3.5a) N times. Notice that we cannot immediately identify this product of the Lax operators with the monodromy matrix T N since the first and last Lax operator on the left-and right-hand side, respectively, depend on an alien spectral parameter v + rather than u + . However, the correct dependence can be recovered owing to the (right) factorization property of L Then we find Next we notice that the intertwiner restores the correct spectral-parameter dependence in the 21-and 22-elements of the Lax operator on the left-hand side, where the second-row entries of the matrix operator N + 1 are the same as in the Lax operator Substituting this relation into the left-hand side of Eq. (4.4) and freely dragging (R + 1 ) −1 to the right of T N −1 , since the latter does not depend on the first lattice site, we project on the 22 entry of the resulting matrix equation with the two-vector ↓|, to find after multiplication by R + 1 from the right Having established the commutativity of S + N R + 1 with D N (u), we demonstrate in Section 4.3 that the former taken at the point v + = 0 obeys a finite-difference Baxter equation (4.28) thus allowing us to identify it with the Baxter operator The integral kernel (see Eq. (2.12)) of the Baxter operator Q + (u) takes a rather simple form where z 0 = 0. The diagrammatic representation for the kernel (4.10) is shown in Fig. 2. Using the Feynman trick to combine the propagators attached to the same integration vertices, w n , and taking into Eq. (B.4) for the reproducing kernel one gets the following representation for the Baxter operator which is particularly useful in deriving the light-cone Hamiltonian (1.24). 0

Baxter operator Q −
Making use of the scalar product (2.8), we can find the conjugate Baxter operator Q − (u), A kernel of this operator can be easily obtained from (4.10) where w 0 = 0. In the diagrammatic form it is shown in Fig. 3. At the same time, the operator Q − can be constructed explicitly by means of the defining relation (3.5b) for the intertwiner R − . The main steps in its construction mimic the derivation of the operator S + N R + 1 commuting with the matrix element D N of the monodromy operator T N that we deduced in the preceding section. By analogy with (4.1), we define the operator that moves the spectral parameter v − from the right-to the left-most argument in a string of Lax operators In the right-hand side of this equation, we can recover the N -site monodromy matrix by using the left factorization of the Lax operator such that To restore the correct spectral parameter in the right-most Lax operator in the left-hand side of Eq. (4.15), we notice that this can be achieved for the (relevant) second column of L N via the equation where Comparing it with the conjugate to Q + , we can verify that they coincide for the choice of the spectral parameter v − = 0 with the diagrammatic representation for its integral kernel shown in Fig. 3. In the integral form, it reads where we used the results of Appendix B.

Baxter equations
In this section we will show that the operator Q ± satisfy the first order difference equations in the spectral parameter u. We perform the analysis for the Baxter operator Q − and merely state the result for its hermitian conjugate. In our consideration 8 we will follow the procedure developed in Ref. [32]. The stating point is the relation (3.5b) and the factorized form of the Lax operator [33,34] analogous to (3.2) where Multiplying (3.5a) from the left by Z −1 1 and from the right sequentially first by Z 2 z iv + −iv − 2 and then by W −1 2 followed by z iv − −iv + 2 , we find that the defining equation takes the form We can also derive the Baxter equation for the operator Q − making use of the invariance of the monodromy matrix under gauge rotations of the Lax operators, For detail, see Ref. [34,28]. . (4.23) As can easily be seen from the explicit form of the matrix in the right-hand side of Eq. (4.23), it can be brought into a triangular form by setting the spectral parameter v − to zero. This defines the main building block of the recursion, Taking the product of N − 1 of these local relations with increasing value of indices enumerating the quantum spaces in the chain from left to right, we can identify the result with the Baxter operator up to boundary multiplicative factors. To complete the right boundary, we use the following result such that sandwiching the emerging matrix equation between the states ↓| and |↓ , we immediately find Introducing the Baxter operator via Eq. (4.19), we recover a finite-difference equation that it obeys, i.e., where we used the commutativity of the two operators in the right-hand side that was established in the previous section. Taking the hermitian conjugate of this equation, it yields the Baxter equation for operator Q + (u), Being first order finite-difference equations, both of them can be easily solved in the operator form with the result .
Here we chose a normalization that matches the integral representations (4.11) and (4.20) for Q ± .
Herefrom it is obvious that Q ± are conjugate to one another (4.12) for real spectral parameter u by virtue of the fact that the Sklyanin's zeroes u n are hermitian, u † n = u n .

Hamiltonians from Baxter operators
Finally, having constructed the Baxter operators, we can find Hamiltonians from their expansion in the vicinity of the points u = ±is. A simple calculation yields the result of interest where An immediate inspection demonstrates that their sum is indeed equivalent to (1.24), An alternative representation that makes the diagonalization of the Hamiltonians particularly straightforward can be read off from Eq. (4.29), with its eigenvalue coinciding with the well-known one-loop energy of N -particle spin-s GKP excitation [14,15] when computed on the eigenfunction Ψ u (2.17) of the operator D N that solves the spectral problem (1.29) in question. We are turning to the construction of the wave functions Ψ u next.

Eigenfunctions
There is a brute force method to find the eigenfunction of the Hamiltonian H N by solving differential equations stemming from D N , see Eq. (2.17). However, this endeavor is hopeless for generic values of N and thus we take a route of devising an algebraic construction of the former.

Recurrence relation
The formalism is based on a recurrence relation originating from the defining relation (4.4). Namely, projecting both sides of the latter on the 22-component of the matrix equation, we immediately find As it stands, this relation does not define a recursion for D N since it gets a contamination from the element B N −1 . We notice however, that we can eliminate the unwanted first term in the left-hand side of this relation by an educated choice of a function Ψ N (z 1 , . . . , z N ) that both sides of this equation act on. Let us take it in the form Ψ N (z 1 , . . . , Choosing the spectral parameter as v + = u − u 1 and introducing the notation The solution to this one-term recursion relation is simple and we find for the eigenfunction and the operators S + n are obtained from S + N by removing n of the rightmost R + -factors, i.e., We introduced an overall normalization constant in Eq. (5.5) for later convenience. The above and is labelled by the zeroes u n .

Integral representation of eigenfunctions
Having found a compact representations for the eigenfunctions (5.5), we are now in a position to cast them in a more explicit integral form that will be extremely useful in proving their orthogonality. This discussion essentially follows an analogous one devised for SL(2) invariant spin chains in Ref. [28].
To start with, let us demonstrate the key manipulations using the two-site eigenfunction as an example, Recalling the integral representation of the Euler Beta function, we immediately find This integral is proportional to the hypergeometric function z −iu 2 −s 2 2 F 1 (s − iu 1 , s + iu 2 , 2s|1 − z 1 /z 2 ). Applying the well-known connection formulas for hypergeometric functions, we can rewrite Ψ (u 1 ,u 2 ) as follows, agreeing with finding of Ref. [16], n−1 ending on sites from z 1 to z n .
The wavefunction admits the form of two plane waves z −iun n with momenta p(u n ) = 2u n scattering on each other by exchanging their momenta with the two-particle S-matrix of GKP excitations 9 [14,15,35] Another integral representation for Ψ (u 1 ,u 2 ) (z 1 , z 2 ), regarding it as a holomorphic function in the upper half-plane, can be deduced by means of the formula (B.7) derived in Appendix B such that This construction can easily be extended to N -site eigenstates. The latter admit a convenient diagrammatic form shown in Fig. 5 for the line-integral representation generalizing Eq.  9 We would like to point out that due to our choice of momenta in the definition of the plane waves, the S-matrix enters with exchanged rapidities since S(u 2 |u 1 ) = S(−u 1 | − u 2 ) compared to Ref. [16]. Figure 5: Diagrammatic representation for N -particle eigenstate (5.15). Here α ± n ≡ s ± iu n stands for the power of the propagator (see Appendix B for conventions). Each vertex corresponds to a coordinate in the pyramid that is integrated with the measure (2.5). The top leftmost vertex is located at w = 0.
where the integrand contains a product of vertices along the left edge of the pyramid (see Fig. 4) which branch down in trees ending on z n 's sites. A few low-order trees read The rest can be easily constructed from the above Feynman graph using the fact that any internal vertex of the tree is parametrized by τ -parameters as V Making use of the representation (5.13) it is easy to verify that as a function of real variables the eigenfunction Ψ u satisfies the following relation where all x k > 0.
Finally, the representation for the eigenfunction, determined by the pyramid in Fig. 5, in terms of multiple integrals over the upper half-plane yields where we introduced the function Since these are eigenfunction of a hermitian operator D N , they are orthogonal to each other by default. However, let us demonstrate this explicitly.

A simple calculation making use of integral identities from Appendix B yields
It can be regarded as a chain rule (B.8) for end-points at w = 0 and w = 0, i.e., where β ± = s ± iv. This will be an instrumental building block of the construction that follows. The orthogonality can be shown inductively, so it is sufficient to demonstrate it for reduction of the scalar product from N to N − 1 sites. First, we have to define the complex conjugate pyramid. It is represented by the same graph as in Fig. 5 where, however, one reverses the direction of arrows on all propagators. As a consequence of this, the integral is accompanied by an additional phase where the first factor in it comes from the complex conjugation of each of N (N −1)/2 Y -functions Y * u (z n , z n+1 |w n ) ≡ e 2πis Y u (w n |z n ,z n+1 ) = e 2πis (w n −z n ) −iu−s (w n −z n+1 ) iu−s , while the second one originates from the conjugation of propagators connecting vertices with the point w = 0, namely, Figure 6: The scalar product of the eigenfunctions (5.15) for N = 3. The points w and w are set to zero. This chain of diagrams exhibits the first level of reduction that yields the delta function δ(u 1 − v 1 ). The parts of the graph shown in red are the once that are involved in application of chain rule (for vertices) or permutation identity (for bonds).
Without loss of generality let us perform the reduction from N = 3 to N = 2 since it can be easily visualized diagrammatically. We will split the induction in two step, i.e., the first one that yields the delta-function and then second one that induces an overall normalization factor. They are shown in Fig. 6 and 7, respectively. We start with the right-most vertex z 3 in Fig. 6 and integrate it out making use of the chain rule (B.8) from Appendix B such that after this step the middle graph gets multiplied by the factor e −iπs a(α + 1 , β − 1 ). Next, we use the permutation identity (B.10) twice and move the (red) bond produced in the previous step all the way to the left such that it connects the points w and w , both set at zero. If we represent this bond using the chain rule backwards as in Eq. (5.18), we recognize it as the scalar product of one-particle eigenfunction that factorizes off from the graph. Therefore, the starting graph in Fig. 6 is equal to the left-most graph in Fig. 7 multiplied by e −iπ(s−iu 1 ) Ψ v 1 |Ψ u 1 . (5.20) Notice that the delta function in the inner product Ψ v 1 |Ψ u 1 forces us to set v 1 = u 1 from now on in the starting graph of the second level of reduction examined in Fig. 7.
In the following step, we start by integrating out the two right-most vertices of the left graph in Fig. 7 located at w (2) 1 and w (2) 1 and produce the graph to its right times the factor Next we move the vertical rightmost propagators (shown in red in Fig. 7) from right to left. This step exchanges the momenta u 1 ↔ u 2 and v 2 ↔ u 1 , respectively. Finally, integrating the rightmost top and bottom vertices at w (1) 2 and w (1) 2 , we acquire extra factors emerging from the chain rules e −2iπs a(α + 3 , α − 1 )a(α + 1 , β − 3 ) . Putting everything together we find that the leftmost graph in Fig. 6 is equal to the rightmost one in Fig. 7 up to an overall factor e −iπ(s−iu 1 ) e −4iπs a(α Figure 7: The second level of reduction of N = 3 scalar product yielding the scalar product product of two-site eigenfunctions. Here w = w = 0. Considering the fact that for a generic N , these two steps yield a phase e −iπ(s−iu 1 ) e −2iπs(N −1) , when combined with the overall one stemming from the complex conjugation, it reduces to the overall phase to the one for N − 1-site inner product, namely, The above construction demonstrates the induction that can be used to prove the orthogonality. For instance, to finish-up with the three-site case, repeating the same steps for Ψ (v 2 ,v 3 ) |Ψ (u 2 ,u 3 ) inner product, one finds Without further ado, we just quote the orthogonality relation for arbitrary N where we symmetrized with respect to the set of quantum numbers v to account for situations with u n = v k for n = k.
In spite of the fact we found an explicit form of the eigenfunctions Ψ u (z 1 , . . . , z N ) in Eq. (5.5), we find it instructive to provide yet another representation which is based on the SoV formalism [26] that is discussed in Appendix D. where we explicitly displayed them as superscripts. A proof of this intertwining relation is presented in Appendix C. This allows us to propose the following definition of the inner product on the real positive half-line

Scalar product on real line
where the measure is Notice that the integrand is positive-definite in the integration domain.
In order to prove that the eigenfunctions Ψ u of the previous section are orthogonal with respect to the inner product (6.3) we have to demonstrate that the operator D N is self-adjoint, for any s ≥ 1/2. This can be verified by assuming that the functions entering (6.3) vanish sufficiently fast as x n goes to infinity and confirming that all boundary terms are absent for x n = x n+1 .

Relation between scalar products
Let us establish a relation between the scalar products (2.8) and ( where the reproducing kernel reads (see Appendix B for details) K(x n ,z n ) = e iπs (x n −z n ) −2s . (6.7) Analogous relation holds for Φ. Using this representation for the functions Ψ and Φ and substituting then into the scalar product (6.3), we get (Φ|Ψ) = Φ|X|Ψ , (6.8) where X is an operator with the integral representation Dz n X (w 1 , . . . , w N ;z 1 , . . .z N )Ψ(z 1 , . . . , z N ) (6.9) and the kernel X (w 1 , . . . , w N ;z 1 , . . .z N ) = (K w 1 ,...,w N |Kz 1 ,...,z N ) , (6.10) written in the form of the scalar product (6.3) of the products of reproducing kernels Taking into account the action of sl(2, R) generators on the reproducing kernel is symmetric with respect to the interchange of its arguments S ±,0 w K(w, x) = S ±,0 x K(w, x), we immediately establish the commutativity of the operator D N (u) with X D N (u) X = X D N (u) .
(6.12) Therefore, it shares the same eigenfunctions Ψ u with the operator D N that was diagonalized in previous sections. One can easily calculate its eigenvalues X(u) from the equation 14) comes from the integration in the vicinity x k ∼ w k . It implies that we can replace the eigenfunction by the leading term in the asymptotic expansion. Obviously, since x n−1 x n , we can replace accordingly the intertwining factor by its leading order form x 2s−1 n such that the integrals in Eq. (6.14) decouple and can be easily calculated individually yielding the result that we sought for Thus for the scalar product of an arbitrary function Φ with the eigenfunction Ψ u , we get This equation provides a relation of the scalar products on the line and half plane and it will be instrumental in calculating various integrals which are hopeless to be computed analytically otherwise. Let us note an important difference between the two scalar products. The scalar product on the upper half-plane is an SL(2, R) invariant scalar product. Namely, the group transformations T s (g) = (⊗T s (g)) N of the wave function is defined as with elements of the group g −1 = a b c d being real unimodular matrices, ad − bc = 1. So that T s (g) are unitary operators with respect to this scalar product, On the other hand, contrary to this, the scalar product on the real positive half-line (6.3) is invariant only under scale transformations g λ , with the group elements g −1 λ = λ 0 0 λ −1 .

Multiparticle transitions
While in the above discussion we analyzed the eigenfunction of the matrix elements E u |O N |0 of the light-cone operators (1.23), to apply our findings to the analysis of N -particle contributions in the Wilson loops, we define an N -particle state in the following fashion with ψ u (x 1 , . . . , x N ) being the wave function of the GKP excitations in the position space [16]. These are related to the eigenfunctions (5.13) found in Section 5 via the relation where N u is a normalization coefficient. Indeed, the state |E u takes the form of the scalar product (6.3) of the operator O N with the eigenfunction Ψ u Obviously it is the eigenstate of the operator D N (u) and, hence, the Hamiltonian. In Ref. [16] the function ψ u was normalized in a manner consistent with the Coordinate Bethe Ansatz, such that for large separation between the GKP excitations, z k /z k+1 1, the initial state consists of noninteracting plane waves with unit amplitude The ellipses stand for the terms z 1 −iu k 1 . . . z N −iu k N accompanied by the GKP S-matrices (5.11) for each permutation P as well as power-suppressed O(z k /z k+1 ) contributions. The coefficient N u that corresponds to the normalization (7.4) takes the form It is evident that for physical quantities all dependence on the normalization coefficient drops out. From a technical point of view it is more convenient and natural to accept normalization which preserves the symmetry function under permutation of the separated variables u 1 , . . . , u N . However, to make comparison with the results of Ref. [16] easier, we present our calculation of the square and hexagon transitions in the normalization (7.5).

Square transitions
To start with, let us address the square transitions. Within the context of the operator product expansion for the octagon Wilson loop discussed in Section 1.1, they define its leading asymptotic behavior. For the χ 2 χ 3 χ 5 χ 6 component, it reads )) . (7.6) When generalized to multiparticle spin-s GKP excitations, the square transition is determined by the overlap integral of the D N eigenfunction and its complex conjugate both evaluated in the same conformal frame connected point-by-point by the propagators for the spin-s excitation which read It is a generalization of the one obtained from the OPE of the Wilson loop for the hole excitation, which is a Fourier transform of the measure µ h (u) in Eq. (1.10). The integration over the bottom wave function with attached propagators can be performed explicitly. Indeed, taking into account Eqs. (6.14), (6.13) and the relation (5.14), we find where the right-hand side is accompanied by the product of the one-particle measures for spin-s GKP excitations It reduces to Eq. (1.10) for s = 1/2 corresponding to hole excitation. The above equation allows us to relate the square transitions to the inner product on the real positive half-line This in turn can be rewritten via the scalar product in the upper half-plane, see Eq. (5.27), and immediately yields the result where the sum goes over all permutations accompanied by scattering matrices of GKP excitations. For N = 1, 2 our expressions for the square transitions (7.12) reproduces those in Ref. [16].

Hexagon transitions
The N -particle hexagon 10 transitions are defined as follows One can simplify the expression for H s (u|v) by performing the integrals with respect to unprimed variables x n as was discussed in the preceding section, see Eq. (7.9), Making a change of variables x n = x n /(1 − x n ) and taking into account the definitions (7.14) and (7.2), we can cast H s (u|v) into the form where g + = +1 0 −1 +1 and the transformation T s (g) is defined in Eq. (6.17). In explicit form, the transformed wave function reads Using the relation (6.16) adopted to the present case we can derive the hexagon transition in terms of the scalar product in the upper half-plane Thus the hexagon transition is determined by the matrix element of the group transformation T s (g + ).
In what follows we consider a more general matrix element where g + γ = +1 0 −γ +1 . In order to calculate it, it is convenient to perform the inversion first, where g I = 0 +1 −1 0 . We immediately find that g I g + γ g −1 I = g − γ = 1 γ 0 1 and, therefore, T s (g − γ ) generates a translation [T s (g − γ )Ψ](z 1 , . . . , z N ) = Ψ(z 1 − γ, . . . , z N − γ) . (7.22) Under inversion, the eigenfunctions Ψ u transforms as follows The diagrammatic representation 11 for Ψ u (z 1 , . . . , z N ) is shown in Fig. 8. Its explicit analytic expression can easily be read off according to the Feynman rules of Appendix B. Taking into account (7.22) we find that up to an overall prefactor, the diagram for the matrix element T γ (u, v) (see the left-most graph in Fig. 9) is given by the mirror reflected diagram in Fig. 6 for the scalar product of two wave functions, where however w = −γ is moved away from zero. The analysis goes along the same lines as the calculation for the norm of the eigenfunctions though the reduction starts from left to right, in the direction opposite to steps taken in Section 5.3 since the initial graph defining the current matrix element is a mirror reflection of the one alluded to above, as shown in Fig. 9. A crucial difference from the norm computation arises after the first level of reduction (cf. the right-most graph in Fig. 6 and the middle graph in Fig. 9) and stems from the fact that the (red) line connecting the points w and w produces (−γ + i0) i(un−vn) rather than the delta function δ(u n − v n ). Taking this into account and collecting all factors due to chain integrations, we finally obtain Sending γ → 1, we find a factorized expression for the N -particle hexagon transition where the leading order one-particle hexagon is a generalization of the corresponding expression for the hole excitation, see Eq. (1.14), . (7.26) This proves the conjecture 12 for multiparticle transitions put forward in Ref. [16].

Conclusions
To conclude, in this paper we systematically studied the open spin chain emerging in the analysis of the OPE for the null polygonal Wilson loop. Using the factorized R-matrices we constructed the commuting Hamiltonians and demonstrated that they are identical to the one emerging from the calculation of the renormalization of the Π-shaped Wilson contour with elementary field insertions. We also demonstrated how they emerge from Baxter operators. The main goal of this analysis was to find the eigenfunction of the Hamiltonian and we constructed them via an iterative procedure for any number of sites. We cast them in multi-dimensional integral form that is particularly useful for the proof of their orthogonality as well as in explicit applications to multiparticle transitions. Making use of the Baxter operators, we cast the eigenstates in the SoV form as well. Finally, we demonstrated that our finding yields the factorized form of the multiparticle hexagon transitions suggested earlier.
There are multiple ways in which we can generalize our results. The consideration can be extended to inhomogeneous case when the spin labels of representation at different sites are not identical. Further, it should be relatively straightforward to extend our analysis to supersymmetric case making use of the factorization property of supersymmetric R-matrices in terms of intertwiners following the formalism of Ref. [36]. that we will do in the next subsection. The subsequent two subsection are dedicated to particular parametrizations of Wilson loop polygons that are suitable for their OPE analysis.

A.1 Coordinates ↔ twistors
Here we will present an embedding of the Minkowski space-time in a six-dimensional projective space and introduce momentum twistors by passing to its spinor representation. This provides with us with a dictionary and an easy conversion of components between different parametrizations.
Let us start with a six-dimensional Euclidean space where the sigma matrices are defined in terms 't Hooft symbols as [37] Σ I AB = (η I AB , iη I AB ) , (A.1) with η andη being the self-dual and anti-self-dual tensors 1 2 ε ABCD η I AB = η I CD , 1 2 ε ABCDηI AB = −η I CD , respectively. In the matrix form, these read where the two-by-two blocks are determined by the conventional Pauli matrices. We pass to the R 2,4 signature of the SO(4, 2) conformal group in the following fashion. Define such that the invariant interval (squared) is given by the determinant .

(A.4)
Pulling X 5 + X 6 from Γ M X M we can identify x µ = X µ /(X 5 + X 6 ) with the Poincare coordinates and X 6 −X 5 = x 2 for the light-like six-dimensional vector X M X M = 0 that provides an embedding of Minkowski space. The above matrix representation is particularly useful to relate it to momentum twistors. The latter are conventionally defined as [38] such that a line in twistor space corresponds to a point in (dual) space-time (A.4) Figure 10: Reference square.
or explicitly where we use a standard notation for the two-dimensional Levi-Civita tensor normalized as ε 12 = −ε˙1˙2 = 1 and the four-dimensional coordinate is expressed in terms of components of the momentum twistors as Notice that in the above construction the infinity bitwistor is defined as such that

A.2 Reference square
The use of the restricted kinematics for scattering amplitudes implies that the bosonic part of the Wilson loop contour resides on a two-dimensional plane with two-dimensional vectors parametrized by light-cone coordinates x = (x − , x + ). Conformal parametrization of higher point polygons relies on the definition of a reference square that we choose to have the following coordinates of its vertices where ε → 0. We would like to take full advantage of linear realization of conformal transformations in the twistor formulation, so that we use the embedding coordinates where the top-right corner is a two-by-two matrix of the Poincare coordinates in question and read Since the twistor space is projective, we can rescale them individually and observe that in the limit ε → 0, they are invariant under the conformal S-transformation (A.15). If we choose not to stick to the twistor ↔ coordinate equivalence and thus loose the direct link to Poincaré coordinates, one can use projective invariance to pull out ε −1 from the Z 3 and Z 4 and set ε = 0 from the get-go. This gives Though we are not able to get the absolute x ± n coordinates from these twistors because the denominator in Eq. (A.12) taken literally is singular, we nevertheless can extract the distances from the square of the interval where the last equality is valid for two-dimensional kinematics only. In the parametrization used below, the x − coordinates will be a function of τ , while x + of σ parameters of the conformal transformation, so the identification will be straightforward (up to an overall constant factor that cancels in the product of x + and x − ). Notice that with the twistor defined as in (A. 16 and are trivially equivalent to the original ones (A. 16).
x 7 x 8 x 9 x 10 Figure 11: Octagon and decagon with reference squares.

A.3 Octagon and decagon
Let us now address the tessellation of the octagon and decagon used in Section 1.1. Starting with the octagon, we choose the twistors defining its contour (see Fig. 11) as with the reference square formed by Z 1 , Z 4 , Z 5 and Z 8 . It is shown by the red dashed contour in Fig. 11 and coincides with the one discussed in the previous section. It is invariant under the transformation (A.14) as was shown above. By acting with the S-transformation in the top-right corner cusp, we can get all inequivalent octagons. To get the same parametrization, we can act instead with an inverse of this S-transformation on the bottom-left cusp of the polygon, with the transformed momentum twistors being The coordinates of the cusps corresponding to the above twistors can be read of from Eq. (A.12), These are given in Euclidean space. To convert them to the Lorentzian signature (which is shown in the figure), one has to drop all minus signs. For the decagon, we choose the twistors as In this case we have two reference squares, the red and the blue one. The red one is invariant under the transformation (A.14), while the blue one is invariant under the following matrix multiplication (A.23) Therefore, we can parametrize all inequivalent decagons by acting with S on the cusps in the top-right corner of the decagon, or which is equivalent to it, the inverse S −1 of the bottom-left cusp x 3 , also we act with S on the cusp x 9 , i.e., Such that the corner coordinates are

B Feynman graphs
The basic ingredient of the Feynman rules is the propagator, illustrated graphically in Fig. 12. All calculation involving it are easily performed making use of the Schwinger representation and Fourier representation (and its inverse) for a holomorphic function Ψ(z) in the upper complex half-plane, make up for all results required for performing explicit computations efficiently.
By virtue of these results, we immediately find the so-called reproducing kernel K that obeys the following property which valid for α + β = 2s. Two properties are particularly instrumental in manipulations performed in the main body of the paper. They are: • the chain rule shown in Fig. 12, and Figure 13: Diagrammatic representation of the permutation identity (B.10) (on the left) and its form in the limit when z 1 → ∞ (on the right).
While the former can be easily obtained by first applying the Schwinger parametrization for the propagators (B.1) and subsequent use of the integral representation for the delta function (B.2), the latter requires a few extra steps. Namely, first using Eq. (B.7) with the Ψ-function given by the product of the remaining two propagators and joining the latter via the Feynman parametrization with a parameter σ, we can conveniently change these two integration variables to ξ and ρ, such that τ = 1/(1 + ξ) and σ = 1/(1 + ρ) with both of them varying in the semi-infinite interval [0, ∞). This allows us to rescale them as follows ξ → ξ(w 2 − z 2 )/(w 2 − z 1 ), ρ → ρ(w 2 − z 2 )/(w 1 − z 2 ) and evaluate the resulting ξ integral in terms of the Euler Beta function B(α − , α + ). The final ρ integral depends on a single cross ratio R = (w 2 −w 1 )(z 2 − z 1 )/[(w 2 − z 2 )(w 1 − z 1 )] and yields the hypergeometric function 2 F 1 , This representation immediately yields (B.10). Another identity used in the main text emerges when one sends the point z 1 to infinity. It is shown graphically in the right panel in Fig. 13.

C Intertwining relations
In this Appendix we provide proofs of intertwining relation (6.2) and eigenvalue relation (6.13).

C.1 Intertwining relation (6.2)
Let us perform a similarity transformation V −λ N T (s) (u) V λ N on the monodromy matrix (2.14) with V N = z 21 z 32 · · · z N,N −1 . In order to proceed, we will use the following factorized expression for the Lax operator where the two-by-two matrices are A key ingredient of the calculation can be clearly illustrated on a simplest two-site example. Namely, the product of two Lax operators admits the following factorized representation Since the left-and right-most lower-triangular matrices M Namely, we find by moving z 21 from right to left. Therefore, by setting λ = 1 − 2s we realize that the similarity transformation yields the change of the spin variable s → 1 − s in the middle matrix block namely, while the left-and right-most lower-triangular matrices remain unchanged, for every internal matrix we have the replace s → 1 − s. Thus, we have where left-and rightmost lower-triangular matrices in this equation are M   C.2 Formula (6.13) To start with, let us point out that the eigenfunction Ψ u (x 1 , . . . , x N ) defined in Eq. (5.13) can be written in the nested form corresponding to the diagram shown in Fig. 4, Here the left-most Λ N is an operator with the integral kernel acting on N − 1 coordinates of a test function Ψ(x 2 , . . . , x N ) as follows where the W N factor is defined in Eq. (6.1). The other operators Λ n<N in Eq. (C.10) act only on the last n − 1 variables, i.e., from x N −n+2 to x N . We can prove the following identity where the propagator G s is defined in Eq. (7.8), the one-particle GKP measure µ s is given in Eq. (7.10) and the integration measure D N x is introduced in Eq. (6.4) with D N x = dx 1 D N −1 x. This equation is the main building block in constructing a recursion procedure to prove Eq. (6.13). The action of the Λ N operator on the product of propagators can be calculated explicitly yielding so that the right-hand side of Eq. (C.12) can be rewritten as In the left-hand side of (C.12) we use the definition of the Λ N -operator given in the second line in Eq. (C.11) and obtain x k x k−1 dy k (y k − x k−1 ) s+iu−1 (x k − y k ) s−iu−1 Ψ(y 2 , . . . , y N ) .

(C.15)
Here the integration variables are ordered in the following fashion This constraint can written in terms of step function in two equivalent ways Making use of the second representation, it is possible to interchange the order x and y integrals in (C.15) and get where all J n integrals entering this formula (with 1 ≤ n ≤ N ) can be calculated explicitly yielding with 1 < n < N . Collecting everything together we find that the left-hand side the expression This coincides with right-hand sides of (C.12) and thus proves the validity of this identity. Thus, starting from the relation and using recurrence (C.12) we obtain the relation (6.13) in the following form Changing the variables as x n → −x n and using the relation (5.14) along with K w 1 ,...,w N (x 1 , . . . , x N ) = e iπsN N n=1 G s (−w n , x n ), it proves the relation (6.13) in question.

D Separation of Variables
In this Appendix we provide a complementary construction of the wave functions that has its roots in the SoV formalism by Sklyanin [26]. It is based on the construction of eigenfunctions of the element B N of the monodromy matrix and its subsequent use in transformation to a factorizable representation. with explicitly factorized overall momentum operator p = N n=1 iS 0 n , whose commutation relations follow from the Yang-Baxter equation with a rational R-matrix, Then we find 5) with the N − 1-site monotromy matrix T N −1 (u + , u − ) = L 1 (u + , u − ) . . . L N −1 (u + , u − ). Projecting out the 12-component of this matrix equation, we find an operator identity relating N -and N − 1-site operators Acting on a function with appropriate properties that allows us to get rid of the first term in the left-hand side, it will define a recursion relation for B N . Namely, acting on a function independent of the z N coordinate does the job of eliminating A N −1 out from the equation. Then, setting v − = u − x 1 , we get the first level of the recurrence relation Repeating this N − 2 times, we finally come to the equation for the last coordinate z 1 , where the differential operator acting on it admits a very simple form, i.e., B 1 (u) = iS − 1 . The eigenfunction of the B 1 operator with the eigenvalue p is Ψ(z 1 ) = exp(ipz 1 ). Therefore, the eigenfunction of B N is U p,x (z 1 , . . . , z N ) = e −iπsN (N −1)/2 (D.9) where we used the invariance of intertwiners under the shift of of their arguments. This expression can be cast in the integral form and, again, admits the structure of a pyramid . . in terms of the Y -function (5.16). This eigenfunction of B N coincides with the one found in Ref. [28] dedicated to SoV of the sl(2, R) invariant spin chain: in both cases here and there, one diagonalizes the same off-diagonal element of the monodromy matrix thus yielding identical results. Therefore, we can borrow the result of that work to our present needs. As can be verified in a straightforward fashion, this eigenfunction obeys the following multidimensional Baxter equation D N (x n )U p,x (z 1 , . . . , z N ) = (x n − is) N U p,x−ien (z 1 , . . . , z N ) , (D. 10) with the unit vector e n possessing only one nonvanishing component (e n ) k = δ nk . As we will see below, this will result in similar equations obeyed by the wave function in the SoV representation.

D.2 Transformation to SoV
The completeness condition for the wave functions U p,x reads  and completes the construction of the normalized eigenfunction in the SoV representation.