Completeness of SoV Representation for $\mathrm{SL}(2,\mathbb R)$ Spin Chains

This work develops a new method, based on the use of Gustafson's integrals and on the evaluation of singular integrals, allowing one to establish the unitarity of the separation of variables transform for infinite-dimensional representations of rank one quantum integrable models. We examine in detail the case of the $\mathrm{SL}(2,\mathbb R)$ spin chains.


Introduction
The field of quantum integrable models takes its roots in the seminal work of Hans Bethe [2] on the XXX Heisenberg chain who developed the so-called coordinate Bethe ansatz method allowing one to construct the eigenvectors and eigenvalues of the mentioned Hamilton operator by means of combinatorial expressions involving auxiliary parameters. In order to obtain an eigenvector, one needs to impose certain constraints on these parameters, the so-called Bethe ansatz equations. Over the years, the method was refined and applied to numerous other models, such as the XXZ Heisenberg chain [31], the δ-function Bose gas [26], or the Hubbard model [27], so as to name a few. In the late 70s, the method was raised to a higher level of effectiveness by Faddeev, Sklyanin, Takhtadjan [42], thus becoming known as the so-called algebraic Bethe ansatz. This new approach provided an algebraic setting allowing one to connect quantum integrability to the representation theory of quantum groups, which had several advantages. To start with, the construction of the eigenvectors of a given integrable model was significantly simplified, hence allowing to address more involved problems such as the calculation of norms [21] and scalar products [43] of Bethe vectors and, subsequently, the one of correlation functions [17,20]. Moreover, the method allowed one to significantly enlarge the family of known integrable models, see, e.g., the review [24], and in particular efficiently and systematically address the question This paper is a contribution to the Special Issue on Mathematics of Integrable Systems: Classical and Quantum in honor of Leon Takhtajan.
The full collection is available at https://www.emis.de/journals/SIGMA/Takhtajan.html of constructing the eigenvectors of the higher rank integrable models [23]. However, it soon turned out that the method had also its limitations in that not all quantum integrable models were within its grasp, the quantum Toda chain being a prominent example thereof. In 1985, Sklyanin pioneered a new technique allowing one to address the calculation of the spectrum of this model: the quantum separation of variables [40]. He developed several aspects of the method in [39,40], and this progress was subsequently continued by Kharchev and Lebedev [18,19] in the case of the Toda chain. Derkachov, Korchemsky and Manashov [5,6], Bytsko and Teschner [3], Silantyev [36] and, more recently, Maillet and Niccoli [28] pushed the development of the method in the case of other, more involved, models (see also [13,33,34]). Recently, many important results have appeared in this area [4,11,12,29,30].
In fact, there are nowadays many indications that the quantum separation of variables is a much more general technique for solving quantum integrable models that encompasses the algebraic Bethe ansatz [5] and provides one with the quantum analogue of the classical separation of variables technique.
In precise terms, the quantum separation of variables consists in exhibiting a map U between an auxiliary Hilbert space h sov and the original Hilbert space h org on which a given model is formulated. This map should be unitary so as to ensure the equivalence of Hilbert space structure and, above all, such that it strongly simplifies the form taken by the spectral problem associated with a given quantum integrable Hamiltonian. More precisely, integrability of a given quantum Hamiltonian means that there exists a commutative subalgebra in the space of operators {H k } containing the Hamiltonian. Thus, the spectral problem associated with the original Hamiltonian is, in fact, a multi-parameter spectral problem, in that each eigenvector is associated with the tower of eigenvalues of the H k s. Now the role of the map U is to realise the unitary equivalence between h sov and h org in such a way that the original multi-dimensional (because the Hamiltonians have a non-trivial structure) and multi-parametric spectral problem on h org is reduced into a multi-parametric (because one has to keep track of all the eigenvalues) one-dimensional spectral problem on h sov . This thus explains the separation of variables terminology. In fact, this one-dimensional spectral problem corresponds to the resolution of the so-called Baxter T − Q equation associated with the model, proving in this way a remarkable bridge between the spectrum, the T -part, and the eigenvectors, the Q part.
Several ingredients are needed so as to implement the separation of variables program as described above. First, one should construct a map U satisfying to the desired requirements and then show that it indeed corresponds to a unitary map between the Hilbert spaces. In fact, the construction of U can be dealt with by exploiting the Yang-Baxter algebra underlying the integrability of the model. A first method was suggested by Sklyanin in [40]. Later, an alternative construction was proposed in [5], where, in particular, it was pointed out that U can be constructed by using the Baxter Q operator associated with the model. For the first time, the idea of a connection between the Baxter Q operator and separation of variables was apparently formulated in the work [25].
This last idea was later generalised in [28] to other conserved quantities, in the case of models having finite-dimensional local Hilbert spaces. To be more precise about the construction of U, we recall that the original Hilbert space h org , where the model is formulated and h sov , where the separation of variables takes place can be identified with appropriate L 2 spaces h org = L 2 (X , dν org ) and h sov = L 2 (Y, dµ sov ). This is a very general setting which allows X , Y to be finite, discrete or continuous. Upon such an identification of the Hilbert spaces, the map U is defined as an integral transform acting on smooth, compactly supported functions on Y: The functions Ψ y (x) describing the integral kernel of the transform can be thought of as the ana-logues of the function e iyx giving the integral kernel of the Fourier transform. That case, in fact, corresponds to X = Y = R and both dν org and dµ sov coinciding with the Lebesgue measure. In the above case, just as {x → e iyx } corresponds to the system of generalised eigenfunctions of translation operators, {x → Ψ y (x)} corresponds to the system of generalised eigenfunctions of a commutative operator subalgebra of the representation of the Yang-Baxter algebra which gives rise to the original model of interest. The construction of U hence boils down to the construction of this system of eigenfunctions, which become possible since it is reduced to solving hypergeometric like problems [40], viz. first order finite difference equations in several variables. In fact, the very structure of the Yang-Baxter algebra which allows one to construct Ψ y (x) in the first place, does also ensure that, by construction, U fulfills the desired requirement of simplifying the original spectral problem. However, unitarity is a completely different issue. It boils down to proving the orthogonality and completeness of the system Ψ y (x) which can be framed as the following relations understood in the sense of distributions Above δ sov (y , y), resp. δ org (x , x), corresponds to the generalised function which represents the integral kernel of the identity operator on Y, resp. X . Moreover, dµ sov /dy, resp. dν org /dx, is the Radon-Nikodym derivative of µ sov , resp. ν org , in respect to the canonical measure on Y, resp. X . The technique for proving unitarity of U, viz. (1.1)-(1.2), strongly depends on the dimension of the original Hilbert space h org . If h org is finite-dimensional, checking unitarity amounts to a simple comparison of dimensions between h org and h sov . However, many of the physically interesting quantum integrable models are defined on an infinite-dimensional Hilbert space h org . There, unitarity is a much more delicate issue. In fact, unitarity was first established for the Toda chain case by using harmonic analysis of Lie groups techniques [35,45]. However the methods which were used to establish this were quite sophisticated and hardly generalisable to the more complex quantum integrable models. The first step towards proving unitarity, in a simpler and systematic way, was achieved in [5], where a quantum inverse scattering based technique for proving the isometry of U was invented. Then, [22] developed a technique, solely based on the use of natural objects for the quantum inverse scattering, allowing one to prove rigorously the isometry of U † in the case of the Toda chain. All together with the results of [5], this construction ensures the unitarity of U. An even more efficient method allowing one to establish the isometry of U † was proposed recently by the authors in [8]. Around the same time, the work [9] has connected certain scalar products of functions being the building blocks of U to Gustafson integrals [14].
In the present work, we push further this link and use the relation to the Gustafson integrals, along with the closed formula for the latter, so as to propose a novel and remarkably simple method for proving the unitarity of the map realising the separation of variables for the higher spin, non-compact, XXX chains. While focusing on this example, we are deeply convinced of the method's generality and hence applicability to many other quantum integrable models possessing infinite-dimensional local Hilbert spaces and which are solvable by the quantum separation of variables. In order to illustrate the main features of the method without obscuring them by technicalities of the model, as a warm up to our main result, we illustrate how it works in the case of the Toda chain.
The paper is organised as follows. Section 2 outlines, on formal grounds, the key ideas of our method in the case of the Toda chain. Then, Section 3 introduces the XXX non-compact spinchain model along with the main notations. In particular, it defines the operators U of interest to the analysis and establishes their isometry. Finally, Section 4 establishes the isometry of their adjoint, and hence completeness of the underlying system of functions giving rise to their integral kernels.

Preliminaries
In this section we illustrate some details of our approach on the example of the open Toda chain [15,37]. which is a one-dimensional system of N particles on the line associated with the Hamiltonian The model is integrable and can be solved by the quantum inverse scattering method (QISM) [41,42]. For further discussion it is important that eigenfunctions can be constructed iteratively [18,19], The measure µ N (λ) -the Sklyanin measure -is given by a product of Γ-functions Finally, the one-particle eigenfunctions are given by plane-waves Ψ λ 1 (y) = e iyλ . Note, that equation (2.1) is nothing other as the expansion of the N -particle function Ψ λ N (x 1 , . . . , x N ) over the product Ψ γ N −1 (x 1 , . . . , x N −1 ) · Ψ Λ−Γ 1 (x N ). The expansion coefficients are given by products of Γ functions. This property -the possibility to find the expansion coefficients of N particle functions over N −1 particle functions -is very important since it allows one to prove orthogonality and completeness relations for the eigenfunctions using induction on N . 1 Indeed, for N = 1 the eigenfunctions obviously form an orthogonal and complete system. Let us assume that the following identities hold for a certain N ≥ 1 and try to prove that these identities hold for N + 1 as well.
Let us start with equation (2.2a). Substituting Ψ λ N +1 in the form (2.1) into this equation and integrating over x one gets that the orthogonality condition is equivalent to the following equation Next, replacing N → N +1 in (2.2b) and projecting both sides on the functions Ψ γ N (x) and Ψ γ N (x) one gets that the completeness relation is reduced to the following identity Of course, both of these relations should be understood in the sense of distributions. Thus the problem of establishing the orthogonality and completeness relations for the eigenfunctions of the open Toda chain is equivalent to proving these two integral identities. The problem is greatly simplified by the following remarkable result due to R.A. Gustafson [14,Theorem 5.1] where Re(α k ) > 0, Re(β k ) > 0 for all k. Note, that the integral in the l.h.s. (2.6) is exactly the integral appearing in (2.4) which, therefore, can be brought to the form The proof of this identity is already rather straightforward. Some details can be found in Appendix A. It takes a little more work to prove the identity (2.5). Having put α N +1 = L and β N +1 = tL, t > 0, in (2.6) and sending L → +∞ one arrives at the reduced version of the integral (2.6) which takes the form [9] where iγ N +1 = −iγ N +1 = L, one can evaluate the integral using equation (2.8). Then, after some algebra, one can show that (2.5) is equivalent to the following identity We recall here that equations (2.7) and (2.9) have to be understood in the sense of distributions and relegate further details to Appendix A.

Spin chains: operators and eigenfunctions
In this section we construct the representation of separated variables for generic spin chain models and recall some elements of the quantum inverse scattering method (QISM) relevant for our purposes. The spin chain is a quantum system of interacting spins S ± k , S 0 k , k = 1, . . . , N , where the index k enumerates the nodes of the chain. The spin operators are the symmetry generators of the SL(2, R) group [10] which are determined by a real number (spin) s k > 1/2, Operators with the index k act in a Hilbert space associated with the k-th site, H k , which is the Hilbert space of functions holomorphic in the upper complex half-plane, H + . The scalar product in the Hilbert space The measure takes the form where θ(x) is the Heaviside step function. The operators S α k are anti-hermitian with respect to the scalar product (3.1).
Function in H k can be represented by Fourier integrals where the integration runs only over positive momenta Then, in Fourier space, viz. in the momentum representation, the scalar product takes the form One of the main objects of QISM is the monodromy matrix. For the closed/open spin chain of our interest, the monodromy matrix is given by a product of L-operators [42] which are two by two matrices, where u ∈ C is the spectral parameter. The monodromy matrix for the closed chain of length N has the form [42] T while, for the open spin chain, it is given by the following expression [38] T where σ 2 is the Pauli matrix. The entries of the monodromy matrices form commuting polynomial operator families [38,42] which act on the Hilbert space of the model, H N = N k=1 H k . Operators in each of the commuting families share the same eigenfunctions. These systems of functions have proven to be very useful for analysing the properties of spin chains. They determine the so-called Sklyanin representation of separated variables [40]. For the homogeneous chains, viz. ξ a = 0, the corresponding systems for B N , B N and A N operators were constructed in [1,6,7], respectively. Below, we recall these constructions and, on the occasion, extend them to the general case of inhomogeneous spin chains where the ξ a 's are generic. Since the technical details are essentially the same in all three cases, we consider in some detail the B N -system and only quote the results for the other two.
All three families of eigenfunctions can be represented as a convolution of functions of a special type. Namely, let us define a function of two complex variables, z, w ∈ H + , and the variable α ∈ C which is called index, This is a single valued function 2 of z, w which is fixed by the condition arg i/(x + iy) → 0 for x → 0. Some properties of the function D α can be found in [9].

Layer operators
The eigenfunctions can be most conveniently written down in term of the so-called layer operators. Let Λ n+1 (γ, x) be an operator which maps functions of n complex variables to functions of n + 1 variables. It depends on the spectral parameter x ∈ C and the complex vector γ = (α 1 , . . . , α n , β 1 , . . . , β n ) ∈ C 2n . Its action takes the form The weight function µ s (w j ) has been defined in equation In the momentum representation obtained by taking the Fourier transform . . , p n ) e i n k=1 p k z k dp 1 · · · dp n the action of the layer operator can be expressed as where k are the "loop" momenta, p k = q k + k − k−1 and 0 ≡ 0, n ≡ q n+1 and the factor λ n reads .
It can be checked that the operators Λ n satisfy the permutation identity The derivation is based on integral identities for the functions D α which can be found in [7]. Next, for the spin chain of length N , we introduce the following combinations of spins and impurities and define the vector γ N ∈ C 2N −2 : It can be shown, see [7], that the operator These two properties of the layer operators are crucial for constructing the eigenfunctions of the operator B N (u).

Eigenfunctions
Let us define a function of N complex variables Following the lines of [22] one can show that the integrals arising from the action of the Λ k 's in (3.8) converge absolutely and that integrations can be performed in an arbitrary order. Due to the properties of the layer operators, equations (3.5) and (3.7), the function Ψ N p,x is a symmetric function of x 1 , . . . , x N −1 which satisfies the equation In the momentum representation, the function F Ψ N p,x (q 1 , . . . , q N ) is given by the convolution of the layer operators in the momentum representation, equation (3.4), acting on the function δ(p − q).
Our aim is to show that the functions Ψ N p,x , p > 0, x ∈ R N −1 form a complete orthogonal set in the Hilbert space H N . It is straightforward to check this statement for N = 1, 2. The proof for general N is more involved and presents the main task of this paper.
For real x and p the functions Ψ N p,x (z) do not belong to the Hilbert space H N . However, they allow one to define a linear transform from the Hilbert space H B N defined below into H N : where we agree upon . (3.9) To start with, given a smooth, compactly supported function ϕ on R + × R N −1 , one introduces the transform Theorem 3.1. For any smooth, compactly supported function ϕ on R + × R N −1 , T B N ϕ ∈ H N and the following relation holds As such, T B N extends to a linear isometry T B N : One may already draw several consequences from this theorem. First of all, (3.11) ensures that the system of functions Ψ N p,x , p ∈ R + , x ∈ R N −1 forms an orthogonal system in H N , viz. that where the multi-dimensional Dirac delta-function has been introduced in (2.3) while µ B N (x) has been defined in (3.9). This corresponds to the orthogonality relation for the system Ψ N p,x (z) . The r.h.s. of this equation, I(z, z ), is the kernel of the unit operator in H N -the so-called reproducing kernel, see, e.g., [16] -which takes the form I(z, z ) = N k=1 I k (z k , z k ), where It means that for any Ψ(z) ∈ H N the following identity holds The unitarity of T B N , viz. that the map has dense range, will be established in Section 4, hence leading to the main result of the paper.
Proof . In order to prove the theorem, we first establish that (3.11) holds for smooth, compactly supported functions ϕ on R + × R N −1 . For that purpose, let us introduce the regularised function Ψ N, p,x (z) which is obtained from Ψ N p,x (z) by giving small positive imaginary parts to all variables x k , x k → x k + i k , k > 0 and by changings N →s N =s N + N k=1 k in the definition of the vector γ N , equation (3.6). Further, let One may readily check that T B, Let us demonstrate that one can invoke Fubini's theorem to get The scalar product of the regularised functions Ψ N, p,x may be computed in closed form as [7] Ψ N, p ,x , Ψ N, , (3.12) and kj = k + j and κ N = κ N (s ).
In order to obtain this result it is convenient to perform calculation in the momentum space representation. Using the momentum representation for the layer operators (3.4) one can obtain an expression for C (3.13) Here the function f is a product of linear combinations of momenta ij and p raised to some powers. It is important to note that all these combinations are positive and that the integration runs over a compact region X. Performing integrations in a special order using the integral identities for the product of the propagators, see, e.g., [6,7,9], gives the expression (3.12).
Of course one has to justify that the order of integrations does not influence the answer. To this end we note that the integral of |f | can be written in the form where R(x, x , ) is some nonsingular factor given by a product of Γ functions. The function f ( , p, 0, 0, { ij }) ξ 1 =···=ξ N =0 is positive and the integral is a particular case of equation (3.13). Thus this integral can be evaluated, as was discussed before, in a closed form, see equation (3.12 where we recall that ϕ is a smooth function with a compact support. For , → 0 + the integral in the r.h.s. can be easily estimated, see Appendix A for the details, resulting in and µ B N −1 is as introduced in (3.9). Since at → 0 + , one has that T B, N ϕ → T B N ϕ (z) almost everywhere, it follows from Fatou's theorem that Thus, the function T B N ϕ belongs to the Hilbert space H N . Finally, taking into account that
The layer operators satisfy the following permutation relation [1], and define the function Φ (σ) By virtue of equation (3.16) Φ N σ,x is a symmetric function of x 1 , . . . , x N . It can be shown, see, e.g., [7], that the operator A N (x 1 ) + σB N (x 1 ) annihilates the layer operator Λ (σ) and, hence, the function Φ N σ,x satisfies the equation Thus the function Φ N x ≡ Φ N σ=0,x diagonalizes the operator A N (u). For the separated variables x with small positive imaginary parts the function Φ N σ,x has a finite norm. Indeed one can find for the scalar product of two Φ N functions [1] .

(3.18)
Here X = N k=1 x k , Y = N k=1 y k . For real x the functions Φ N x are orthogonal to each other (see Appendix A for more details) where .
Upon repeating the argument given in the previous subsection, one can prove the following statement: The transformation T A N defined for smooth, compactly supported functions χ on R N : extends into a linear isometry from H A N into H N . In particular it has unit operator norm T A N = 1 and satisfies It will be show in Section 4 that T A N is, in fact, an unitary map between the corresponding Hilbert spaces.

The layer operators (3.21) satisfy the permutation relation [6]
and is nullified by the operator B N (x) (3.22) where the vector η N is given by equation (3.17). Given z = (z 1 , . . . , z N ), define the function where S = N k=1 s k , x = (x 1 , . . . , x N −1 ) and the normalisation factor is The function Υ N p,x is a symmetric even function of x 1 , . . . , x N −1 which is well defined for | Im(x k )| < min k s k . By virtue of equation (3.22), the function Υ N p,x diagonalizes the operator B N The functions Υ N p,x are orthogonal to each other for real separated variables, x ∈ (R + ) N −1 . Namely, We are now in position to formulate the theorem: The transformation T B N defined for smooth, compactly supported functions φ on R + × (R + ) N −1 that are symmetric in respect to the last N − 1 variables as extends to a linear isometry from H B N into H N . In particular, it has unit operator norm T B N = 1 and satisfies We will show in Section 4 that T B N is an unitary operator.

Completeness
In the previous section we constructed three systems of functions, Ψ The backward arrow is dispensable here, but we consider it first because its proof is most transparent and all other proofs follow the same scheme. It was shown in the previous section that R Since the map T A N is an isometry, by assumption, it maps ker(S BA ) → ker T B N . Our immediate aim is to show that ker(S BA ) = 0.
We prove the following statement: Proof . First, we calculate the action of the operator S BA on the space of smooth functions with a compact support, χ(x). The action of T A N on a function χ(x) is given by equation (3.20). In full similarity with the construction in Section 3.1.2, we define the regularized function T A,σ, Thus we write Moreover, one has The further analysis depends on the remarkable fact that the scalar product of the functions Ψ N p,y and Φ N σ,x+i can be obtained in a closed form [1]: where X = N k=1 x k and Ξ = N k=1 ξ k , and E = N k=1 k . That is By assumption the function χ is nonzero only in a compact region. Therefore the function ϕ σ, (p, y) grows no faster that some power of y for large y while at large p it decays exponentially fast ∼ exp{− Im(σ)p}). Taking into account that the measure µ B N −1 (y) decays exponentially fast for large y one concludes that the normalisation integral for ϕ σ, converges Moreover, substituting the expression for ϕ σ, , equation (4.6), one can change the order of integration and integrate first over p and y. The momentum integral is trivial and produces the factor Γ(i(X − X) + 2E)(2 Im σ) i(X−X )−2E , while the integral over y can be calculated in closed form [14,Theorem 5.1], see also equation (2.6). Namely, Thus we get for the norm of ϕ σ, Note that the expression in the bracket is nothing else as the scalar product, Φ N σ,x +i , Φ N σ,x+i H N , see equation (3.18). Finally, taking into account (3.19), see also Appendix A, we obtain that for any smooth function χ with a compact support Since the space of smooth, compactly supported functions is dense in H A N this equation holds on the whole Hilbert space.
The identity (4.2) implies that ker S BA = 0 and hence R T B N = H N , which guarantees the unitarity of the map T B N . The proof of the unitarity of the maps T B N and T B N +1 follows the same lines and is based on the following result: Provided the map T A N : H A N → H N is unitary the following identities, In the proof of these assertions, the main difference from the proof of Lemma 4.1 lies in the type of the Γ-integrals arising in the process. We briefly discuss these differences below.
For the S BA operator the problem is reduced to calculating the norm of the function which is an analogue of the function ϕ σ, , see equation (4.4). The relevant scalar product takes the form 3 .
Above, the symbol ± stands for Calculating the norm one substitutes the function φ σ, in the form One can change the order of integrations and take the integral over p and y first. The integral over p is exactly the same while the other integral takes the form of second Gustafson integral [14, Theorem 9.3], where kj = j + k . We also extended the integral over y k from the positive half-axis to the real line using the symmetry of the integrand with respect to the reflection y k → −y k . Finally, collecting all factors, one finds that the norm φ σ, The last factor in the above equation is the unitary map from is the Hilbert space of holomorphic functions in the upper complex half-plane discussed around equation (3.1). Namely, similar to equation (4.3) we define Here z, x, and z, x, are N and (N +1)-dimensional vectors, respectively, e.g., x = (x 1 , . . . , x N ), x = (x 1 , . . . , x N +1 ), etc. Note also that the parameter σ is the same for the functions Φ N and Φ 1 .
Completely similar to the previous consideration one can show that Again we define the function ϕ σ, (p, y) = Ψ N p,y , Ψ σ, where y = (y 1 , . . . , y N ). The scalar product of the function Ψ N p,y and Φ N σ,x+i ⊗ Φ 1 takes the form (see, e.g., [9]) .
Calculating the norm of ϕ σ, we change the order of integration and first take the integral over y.
It takes the form of N -fold Gustafson's integral [14,Theorem 5.1] that we encountered earlier, see equation (4.7). After some algebra we obtain The analysis of the above expression in the limit → 0 + and σ → 0 is exactly the same as before, see Appendix A. It results in the following expression for the norm that completes the proof of the lemma.
It follows from Lemma 4.2 ker S BA = 0 and ker S B = 0 and, hence, that the operators T B N : H B N → H N and T B N +1 : H B N +1 → H N +1 are unitary provided that T A N is. The final step required to complete the induction on N is to show that the unitarity of the map T A N follows from that of the map T B N . As it was argued earlier in this section in order to prove this statement it is enough to show that the kernel of the operator is trivial. Proof . Let ϕ(p, x) be a smooth function with compact support having the factorised form The linear span of these functions is dense in H B N . The action of the operator S AB on a function ϕ can be represented as follows where y +i = (y 1 +i , . . . , y N +i ) and the scalar product of two eigenfunctions is given by equation (4.5). We also denote the function given by the integral in the above equation by χ (y), i.e., It can be shown, see [22,Lemma 3.1], that the function χ(y) i<k y ik is a smooth function. Our final aim is to show that χ ∈ H B N and that the χ 2 In contrast to the previous cases, the function χ(y) does not decrease fast enough for large y k to justify changing the order of integration over x, x and y in the norm integral. To overcome this difficulty we proceed as follows. Let us define a regularized function χ L (y) as The factor g L (y) has the following properties: (i) g L (y) < 1 for all y, (ii) g L (y) → 1 monotonically as L → ∞ for fixed y, It follows from (ii) that for any bounded region D ∈ R N D |χ(y) − χ L (y)| 2 dµ A N (y) → 0 as L → ∞.
Due to (iii) one concludes that, for finite L, the integral of |χ L (y)| 2 over R N converges Then one derives the following inequality which holds for any L greater than some L M . Since M is arbitrary we get the following inequality which holds for an arbitrary bounded region D. Therefore R N |χ(y)| 2 dµ A N (y) ≤ I. Since due to (i) I ≤ R N |χ(y)| 2 dµ A N (y) we conclude that Thus one has to find the limit of I L at L → ∞. First we note that I L can be written in the form where χ L (y) = g L (y)χ (y). At the last step, the limit → 0 + is taken after the integration. It is possible to do so since the function χ L (y) is bounded, |χ L (y)| < C L for all y, and the measure µ A N (y) decays exponentially fast at large y, so that the measure of the whole space is finite Γ(s k +s j ).
The calculation of the integral of |χ L (y)| 2 follows the familiar pattern: one substitutes the function χ L (y) using (4.11) and then perform first the integral over y. This integral is the reduced version of Gustafson's integral, equation (2.8). Making use of this result one can write I L in the form (4.12) We recall that X(X ) = N −1 k=1 x k (x k ). Due to our assumptions on the function ϕ the integral in (4.12) is restricted to a finite region hence we can expand the function M (L, p, x, p , x ) in series in L −1 where we put → 0 in non-singular terms. At large L, the dominant contribution to the integral over p, p comes from the region p = p and can be easily estimated as see equation (4.10). The factor in the first line of equation (4.13) has the form we encountered earlier, see, e.g., (4.8), and can be handled in the same way as before, see Appendix A for more details. Collecting all factors we obtain Thus one concludes that S AB is a norm preserving map, S AB ϕ 2 hence, ker S AB = 0. Therefore one concludes that T A N is a unitary operator between the Hilbert spaces H A N and H N .

Summary
This work devised a very effective inductive scheme allowing one to prove the completeness of Sklyanin's separated variables which arise in the analysis of the closed and open SL(2, R) spin chain magnets. The method we proposed heavily relies on the use of multidimensional Mellin-Barnes integrals which were calculated in closed form by R.A. Gustafson [14]. The attractive feature of our approach is that it does not depends on the details of the spin chain -spins, s k , and inhomogeneity parameters, ξ k . Moreover, the core identities which are to be used for the closed spin chain or for Toda chain are exactly the same, what stressed a certain generality of our method. Since the Gustafson integrals can be viewed as a special case of the elliptic hypergeometric integrals, see, e.g., [32,44], we believe that our method can be adapted to such models as spin chains with the trigonometric and elliptic R-matrices or non-compact magnets with the SL(2, C) symmetry group.
A Some representations for multi-dimensional Dirac δ-functions .
In the following we show that, in the sense of distributions, it holds .
In other words, given In order to establish the result, one starts by reorganising the integral in (A.2) as p, x, x ϕ p, x ϕ * p, x , .

Thus, U
( ; ) N is antisymmetric in x, x taken singly, and it is smooth and compactly supported in x, x ∈ R N −1 and smooth in a small neighbourhood of zero in respect to , .
Expanding the determinant as a sum over the permutation group and using the antisymmetry in x, x and the smoothness in , of U .
It is thus enough to study the ε → 0 + , → 0 + limit of the model integral in which χ is antisymmetric in x, x taken singly. Observe that, for fixed x, by the Stone-Weierstrass theorem, there exists a sequence of smooth, compactly supported functions ϕ k,a on R such that Next, one observes that Thus, inserting the expansion in the integral, summing up, setting a = 0 in the regular part of the integrand and using the antisymmetry in x, x of χ, one gets that where ∆ (s) x,x is a composite of operators acting on the variables x 1 , . . . , x s , x 1 , . . . , x s ∆ (s) x,x = .
Apart from the term arising in the last line, the integrand is a smooth function. Thus, by changing the variables x → y = This entails the claim. II. Let us define .
We will show that in the sense of distributions the following identity holds where W N is defined in equation (A.1). Namely, given Repeating the same arguments as above one can rewrite (A.4) in the form Finally, taking into account that and using the Stone-Weierstrass theorem one gets the necessary result.