Baxter operators in Ruijsenaars hyperbolic system II.
Bispectral wave functions

N. Belousov$^{\dagger\times}$, S. Derkachov$^{\dagger\times}$, S. Kharchev$^{\bullet\ast}$, S. Khoroshkin$^{\circ\ast}$

$^\dagger$Steklov Mathematical Institute, Fontanka 27, St. Petersburg, 191023, Russia;

$^\times$National Research University Higher School of Economics, Soyuza Pechatnikov 16,
St. Petersburg, 190121, Russia;

$^\bullet$National Research Center “Kurchatov Institute”, 123182, Moscow, Russia;

$^\circ$National Research University Higher School of Economics, Myasnitskaya 20,
Moscow, 101000, Russia;

$^\ast$Institute for Information Transmission Problems RAS (Kharkevich Institute),
Bolshoy Karetny per. 19, Moscow, 127994, Russia

In the previous paper we introduced a commuting family of Baxter $Q$-operators for the quantum Ruijsenaars hyperbolic system. In the present work we show that the wave functions of the quantum system found by M. Hallnäs  and S. Ruijsenaars also diagonalize Baxter operators. Using this property we prove the conjectured duality relation for the wave function. As a corollary, we show that the wave function solves bispectral problems for pairs of dual Macdonald and Baxter operators. Besides, we prove the conjectured symmetry of the wave function with respect to spectral variables and obtain new integral representation for it.

§ INTRODUCTION

§.§ Ruijsenaars system and Baxter $Q$-operators

In [2] we defined Baxter $Q$-operators for the Ruijsenaars hyperbolic system and proved their commutativity. Now we apply the results of [2] to the study of the wave functions of this system.

Denote by $T^{a}_{x_i}$ the shift operator

\begin{equation}
T^{a}_{x_i}:=e^{a\partial_{x_i}}, \qquad \left( T^{a}_{x_i} \, f \right)(x_1,\ldots,x_i,\ldots,x_n) = f(x_1,\ldots,x_i+a,\ldots,x_n)
\end{equation}

and define its products for any subset $I\subset[n] = \{1, \dots, n\}$

\begin{equation}
T^{a}_{I,x}=\prod_{i \in I} T_{x_i}^a.
\end{equation}

The Ruijsenaars system [15] is governed by commuting symmetric difference operators

\begin{equation}
\label{I2}
H_r(\bx_n;g|\bo) = \sum_{\substack{I\subset[n] \\ |I|=r}}
\prod_{\substack{i\in I \\ j\notin I}}
\frac{\sh^{\frac{1}{2}}\frac{\pi}{\o_2}\left(x_i-x_j-\imath g\right)}
\cdot T^{-\imath\o_1}_{I,x}\cdot \prod_{\substack{i\in I \\ j\notin I}}
\frac{\sh^{\frac{1}{2}}\frac{\pi}{\o_2}\left(x_i-x_j+\imath g\right)}
\end{equation}

Here and in what follows we denote tuples of $n$ variables as

\begin{equation}
\bm{x}_n = (x_1, \dots,x_n).
\end{equation}

One can also consider gauge equivalent Macdonald operators

\begin{equation}
\label{I2a}
M_r(\bx_n;g|\bo) = \sum_{\substack{I\subset[n] \\ |I|=r}}
\prod_{\substack{i\in I \\ j\notin I}}
\frac{\sh\frac{\pi}{\o_2}\left(x_i-x_j-\imath g\right)}
\cdot T^{-\imath\o_1}_{I,x}.
\end{equation}

Both families of operators are parametrized by three constants: periods $\bm{\omega}=(\omega_1, \omega_2)$ and coupling constant $g$,
which originally are supposed to be real positive. The equivalence is established by means of the measure function

i≠j^nμ(x_i-x_j),  where  μ(x):=μ_g(x|)=S_2(x|) S_2^-1(x+g|). Here $S_2(z|\bo)$ is the double sine function, see Appendix <ref>.

M_r(_n;g|)   1/√(μ(_n))=
Ruijsenaars operators are symmetric with respect to the sesquilinear scalar product defined by the Lebesque measure $d\bx_n$, while the Macdonald operators are symmetric with respect to the scalar product related to the measure
μ(_n)   d_n,    x_n∈^n.

In this paper, as well as in [2] and unlike the original Ruijsenaars setting, we consider periods $\bo$ and coupling constant $g$ to be complex valued,  assuming that
ø_1 > 0,     ø_2 > 0,    0< g<ø_1+ø_2  and
ν_g= g/ø_1ø_2>0.Denote by $\K(x)$ the following function of a complex variable
(x):=_g(x|)=S_2^-1(x +g^∗/2|)S_2^-1(-x+g^∗/2|)written in terms of the double sine fuction $S_2\bigl(z\big|\bo\bigr)$ and
g^∗=ø_1+ø_2-g.We also use the products of this function

\begin{equation}\label{I6a}
\K(\bx_n,\by_m)=\prod_{i=1}^n \prod_{j = 1}^m \K(x_i-y_j).
\end{equation}

In  [2] we introduced a family of operators $Q_{n}(\lambda)$ parameterized by $\lambda \in \mathbb{C}$ and called Baxter $Q$-operators. These are integral operators

\begin{equation}\label{I14}
\left( Q_{n}(\lambda) f\right) (\bm{x}_n) =  d_n(g|\bo) \cbk\, \int_{\mathbb{R}^n} d\bm{y}_n \, Q(\bm{x}_n, \bm{y}_n; \lambda) f(\bm{y}_n)
\end{equation}

with the kernel
Q(_n,_n;ł)= e^ł(_n-_n)
(_n,_n)μ(_n),     x_j,y_j ∈and normalizing constant

\begin{equation}\label{dconst}
d_n(g|\bo) = \frac{1}{n!} \left[ \sqrt{\omega_1 \omega_2} S_2(g|\bo) \right]^{-n}.
\end{equation}

Here and below for a tuple $\bx_n=(x_1,\ldots,x_n)$ we use the notation $\bbx_n$ for the sum of components
We remark that the constant $d_n(g|\bo)$ was absent in our previous paper. We added it here to simplify formula for $Q$-operator's eigenvalue in the next subsection. The main result of [2] is given by the following theorem.
[2]  Under assumptions I0a, I0b $Q$-operators commute

\begin{equation}\label{QQcomm}
Q_{n}(\lambda) \, Q_{n}(\rho) = Q_{n}(\rho) \, Q_{n}(\lambda).
\end{equation}

The kernels of the operators in both sides of QQcomm are analytic functions of $\l,\rho$ in the strip
|(λ- ρ)| < ν_g.
The product of two $Q$-operators is a well-defined integral operator on the domain of fast decreasing functions, see Proposition 5 and remark after it in [2]. This is also true for the various products of $Q$- and $\Lambda$-operators below.

We also proved that these operators commute with Macdonald operators

\begin{equation}\label{QMcomm}
Q_{n}(\lambda) \, M_r(\bm{x}_n; g|\bo) = M_r(\bm{x}_n; g|\bo) \, Q_{n}(\lambda), \qquad r=1, \ldots, n
\end{equation}

under condition

\begin{equation}\label{I13a}
\Re g < \Re \o_2.
\end{equation}

This suggests that they should have common eigenfunctions.

§.§ Wave functions and their symmetries

Denote by $\LL{n}(\l)$ the integral operator similar to the Baxter $Q$-operator
(n(ł)f)(_n)=d_n - 1(g|) ∫_^n-1d_n-1   Λ(_n,_n-1;ł) f(_n-1)
with the kernel
Λ(_n,_n-1;ł)= e^ł(_n-_n-1)
(_n,_n-1)μ(_n-1),     x_j,y_j ∈and constant $d_{n - 1}(g|\bo)$ given by the formula (<ref>).

M. Hallnäs  and S. Ruijsenaars proved  in [6] that for real positive periods $\bo$ and the coupling constant
$g$,  $0<\Re g<\o_2$,  the function
Ψ__n(_n):=Ψ__n(_n; g|)=n(ł_n)   n-1(ł_n-1)⋯2(ł_2)   e^ł_1x_1
is a joint eigenfunction of Macdonald operators. Similar arguments show that the same statement holds for complex valued periods $\bo$ and the coupling constant $g$.
Under the conditions I0a,  I0b and I13a we have the equality
M_r(_n;g|) Ψ__n(_n)= e_r (e^2πλ_1ø_1, …, e^2πλ_nø_1) Ψ__n(_n).

Here $e_r(z_1,\ldots,z_n)$ is $r$-th elementary symmetric function,
e_r(z_1,…,z_n)=∑_1≤i_1<i_2<…< i_r≤nz_i_1⋯z_i_r,

The analytical properties of the function $\Psi_{\bl_n}(\bx_n)$ with respect to $x_j, \lambda_j$ (in the case of positive periods) are investigated in detail in [6]. Here we only remark that from definition (<ref>) and formula for the kernel (<ref>) the variables $x_j$ determine the positions of poles of functions $K(x_j-y_i)$, see (<ref>). In the case $x_j \in \R$ the sequences of poles are separated by the integration contours. The closest poles to the contours are $y_i = x_j \pm \imath g^*/2$.

Therefore, the wave function $\Psi_{\bl_n}(\bx_n)$ extends to the analytic function of $x_j \in \mathbb{C}$ in the domain

\begin{equation}
| \Im x_j | < \frac{g^*}{2}, \qquad j=1, \dots,n.
\end{equation}

From the convergence of integral representation (<ref>) we also have the  analyticity of the wave function with respect to the variables $\lambda_j \in \mathbb{C}$ in the domain

\begin{equation}
| \Im(\lambda_j - \lambda_k) | <\theta, \qquad j,k =1, \dots,n
\end{equation}

with some positive $\theta$, see Proposition <ref> in Section <ref>.

The short proof of Theorem <ref>, combining the arguments by M. Hallnäs and S. Ruijsenaars [6] and by I. Macdonald [13] is given in Section <ref>.  See also [10] for an analogous proof in the case of hyperbolic Calogero-Sutherland model, which represents a non-relativistic limit of the Ruijsenaars system.

Note that the double sine function is invariant under permutation of $\omega_1, \omega_2$, therefore the wave function $\Psi_{\bm{\lambda}_n}(\bx_n)$ is invariant too. This implies that under condition

\begin{equation}
0 < \Re g < \min(\Re \omega_1, \Re \omega_2)
\end{equation}

it is also an eigenfunction of difference Macdonald operators $M_r(\bm{x}_n; g | \omega_2, \omega_1)$ related to the shift of variables by $\imath \omega_2$, and thus satisfies two pairs of compatible difference equations.

The formula (<ref>) also can be written recursively

\begin{equation}\label{Psi-rec}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = \Lambda_{n}(\lambda_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}), \qquad \Psi_{\lambda_1}(x_1) = e^{{2\pi \imath}{} \lambda_1 x_1}.
\end{equation}

Hence, we call $\Lambda_{n}(\lambda)$ raising operator.

Our approach to study the wave functions in Ruijseenars system is based on the detailed investigation of the relations between the various products of $Q$-operators I14 and raising operators I11 which follow from commutativity relation of $Q$-operators QQcomm.

Firstly, an immediate consequence of QQcomm is the exchange relation for the raising operators, its short proof is given in Section <ref>.
Under assumptions I0a, I0b raising operators satisfy relation

\begin{equation}\label{ll}
\Lambda_{n}(\lambda) \, \Lambda_{n - 1}(\rho) = \Lambda_{n }(\rho) \, \Lambda_{n - 1}(\lambda).
\end{equation}

The kernels of the operators in both sides of ll are analytic functions of $\l,\rho$ in the strip
|(λ- ρ)| < ν_g.

The wave function $\Psi_{\bm{\lambda}_n}(\bm{x}_n)$ is invariant under permutations of $\bm{x}_n$ and $\bm{\lambda}_n$ components.

The symmetry with respect to permutations of $x_j$ follows immediately from the defining formula (<ref>). The symmetry over spectral variables $\l_j$ is a consequence of Theorem <ref>. To apply it we use absolute convergence of the corresponding multiple integral (<ref>) (see Proposition <ref> in Section <ref>) that allows us to change integration order in this integral.

The next exchange relation appears in the limit of the $Q$-operators commutation relation (<ref>) as one of the external kernel parameters $z_n$ tends to infinity. Set
â=a/ø_1ø_2for any $a\in\C$, so that
=(1/ø_2,1/ø_1) ,
ĝ=g/ø_1ø_2,    ≫=ø̂_1 + ø̂_2 - ĝ = g^∗/ø_1ø_2.
Also introduce the function
The operator identity

\begin{equation}\label{QLcomm}
Q_{n}(\lambda) \, \Lambda_{n}(\rho) = \KK(\lambda - \rho) \; \Lambda_{n}(\rho) \, Q_{n-1}(\lambda)
\end{equation}

holds true for $\lambda, \rho \in \mathbb{C}$ such that
| (λ- ρ) | < ν_g/2.
The  proof is given in Section <ref>.  Note that the strip for $\lambda,\rho$ (<ref>) is twice narrower than in the $Q$-commutativity relation (<ref>). Using the iterative representation of the wave function (<ref>) and the relation QLcomm we arrive at the spectral description of the Baxter $Q$-operator. Its proof with the necessary convergence arguments is given in Section <ref>.
The wave function $\Psi_{ \bm{\lambda}_n }(\bm{x}_n)$ is a joint eigenfunction of the commuting family of operators $Q_{n}(\lambda)$

\begin{equation}\label{QPsi}
Q_{n}(\lambda) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n) = \prod_{j = 1}^n \KK(\lambda-\lambda_j) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n).
\end{equation}

The integrals in both sides of QPsi converge if

\begin{equation}
| \Im(\lambda - \lambda_n) | < \frac{\nu_g}{2}(1-\ve), \qquad | \Im (\lambda_k - \lambda_j) | <\theta(\ve), \qquad k, j = 1, \dots, n
\end{equation}

for any $\ve \in[0,1)$ and
θ() = ν_g /4 (n - 1)! e.

The kernels I14 and Lker of the $Q$-operator and raising $\Lambda$-operator are related in two ways. Firstly, the kernel  $\Lambda(\bx_n,\by_{n-1};\l)$ appears in the certain limit of the kernel $Q(\bx_n,\by_{n};\l)$. Besides, the kernel $\Lambda(\bx_n,\by_{n-1};\l)$ equals the kernel $Q(\bx_{n-1},\by_{n-1};\l)$ multiplied by the function

\begin{equation}
e^{\const \l x_n} \prod_{j=1}^{n - 1}K(x_n - y_j)
\end{equation}

regular on the integration cycle. Then the  inductive use of Theorem <ref> allows to relate two integral representations for the wave function: the representation iterative with respect to variables $x_j$ and another one iterative with respect to spectral variables $\l_j$. This can be formulated as duality relation of the wave function conjectured by M. Hallnäs and S. Ruijsenaars in [6].

For this purpose we introduce the counterparts of basic functions I5, I6a and operators I14, I11, but in the space of spectral variables $\l_j$. Namely, similarly to the function (<ref>) we define the function

(ł):=μ_ĝ^*(ł|)=S_2(ł|) S_2^-1(ł+≫|) and their products

\begin{equation}\label{I6d}
\mm(\bl_n)=\prod_{\substack{i,j = 1\\i\not=j}}^n \mm(\l_i-\l_j),\qquad	\KK(\bl_n,\bg_m)=\prod_{i=1}^n \prod_{j = 1}^m \KK(\l_i-\g_j).
\end{equation}

Next we introduce counterparts of the $Q$-operator and raising $\Lambda$-operator that act on functions of spectral variables $\lambda_j$

\begin{equation}\label{I15d}\begin{aligned}
\bigl( \QQ_n(x) f\bigr) (\bm{\l}_n)& = d_n(\hat{g}^*|\hat{\bo})\, \int_{\mathbb{R}^n} d\bm{\gamma}_n \, \QQ(\bm{\l}_n, \bm{\gamma}_n; x) f(\bm{\gamma}_n),\\[7pt]
\bigl(\LLL_{n}(x)f\bigr)(\bl_n)&=d_{n - 1}(\hat{g}^*|\hat{\bo}) \int_{\R^{n-1}}d\bg_{n-1} \, \LLL(\bl_n,\bg_{n-1};x) f(\bg_{n-1})
\end{aligned}
\end{equation}

with the kernels
(_n,_n;x)   = e^x (_n-_n)	(_n,_n)μ̂(_n),

(_n,_n-1;x)   = e^x (_n-_n-1 ) (_n,_n-1) μ̂  (_n-1).

Additionally, in what follows we assume the condition
ν_g^*:= ĝ^* = g^*/ø_1 ø_2 > 0.Then together with the conditions (<ref>), (<ref>) we have two pairs of double inequalities for the coupling constant

\begin{equation}\label{gg-ineq}
0<\Re g < \Re\o_1 + \Re\o_2, \qquad 0< \Re \hat{g}^* < \Re\hat{\o}_1 + \Re\hat{\o}_2.
\end{equation}

Note that for the real periods they are equivalent.  Clearly, with these conditions introduced operators enjoy relations analogous to ll, (<ref>) and can be used for iterative construction of the function

\begin{equation}\label{Psi-rec2}
\hat{\Psi}_{\bx_n}(\bl_n) := \Psi_{\bx_n}(\bl_n;\gg|\hat{\bo}) = \LLL_{n}(x_n) \, \hat{\Psi}_{\bx_{n - 1}}(\bl_{n - 1}), \qquad \hat{\Psi}_{x_1}(\l_1) = e^{{2\pi \imath}{} \lambda_1 x_1}.
\end{equation}

The following statement was conjectured in [6].
The wave function $\Psi_{\bm{\lambda}_n}(\bm{x}_n)$ satisfies duality relation

\begin{equation}\label{duality}
\Psi_{\bm{\lambda}_n}(\bm{x}_n; g|\bo) = \Psi_{\bm{x}_n}(\bm{\lambda}_n; \gg|\hat{\bo}).
\end{equation}

The proof of this theorem is given in Section <ref>. It has several direct corollaries.
The wave function  $	\Psi_{\bm{\lambda}_n}(\bm{x}_n)$ is an eigenfunction of both operators $Q_{n}(\lambda)$ and $\QQ_{n}(x)$

\begin{equation}\label{QQPsi}
\begin{aligned}
Q_{n}(\lambda) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n) &= \prod_{j = 1}^n \KK(\lambda-\lambda_j) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n),\\%[5pt]
\QQ_{n}(x) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n) &= \prod_{j = 1}^n \K(x-x_j) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n).
\end{aligned}
\end{equation}

Under the condition $\Re g < \Re\o_2$ I13a and analogous one $\Re \hat{g}^* < \Re \hat{\o}_2$  the wave function  $\Psi_{\bm{\lambda}_n}(\bm{x}_n)$ is an eigenfunction of both families of Macdonald operators $M_r(\bx_n;g|\bo)$ and $M_s(\bl_n;\hat{g}^*|\hat{\bo})$

M_r(_n;g|) Ψ__n(_n)   =e_r (e^2πλ_1ø_1, …, e^2πλ_nø_1 )Ψ__n(_n),

M_s(_n;ĝ^*|)  Ψ__n(_n)   =e_s (e^2πx_1/ø_2, …, e^2πx_n/ø_2 )Ψ__n(_n).

In other words, the function $\Psi_{\bm{\lambda}_n}(\bm{x}_n)$ solves bispectral problems for pairs of dual $Q$-operators and for a pair of dual Macdonald operators.  In non-relativistic limit the bispectrality property of the wave function of Sutherland model was established in [9, 10]. The  duality MM provides us with two integral representations (<ref>), (<ref>), recursive with respect to $x_j$ and $\lambda_j$ variables correspondingly.

\begin{equation}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = \Lambda_{n}(\lambda_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}) = \LLL_{n}(x_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}).
\end{equation}

For real periods these two recursive procedures can be combined into the third one, symmetric with respect to coordinates and spectral parameters, see Section <ref> for details.

For real positive periods $\o_1, \o_2$ the wave function admits the iterative integral representation

\begin{equation}
\begin{aligned}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) &= e^{\const \lambda_n x_n} Q_{n - 1}(\lambda_n) \, \QQ_{n - 1}(x_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1} ) \\[6pt]
&= Q_{n - 1}(\lambda_n) \, \QQ_{n - 1}(x_n) \cdots Q_{1}(\lambda_2) \, \QQ_{1}(x_2) \, e^{\const (\lambda_1 x_1 + \ldots + \lambda_n x_n)}.
\end{aligned}
\end{equation}

An analogous representation is presented in [5] for the open Toda wave function. The wave function of open Toda chain also admits two integral representations iterative with respect to coordinates $x_j$ and spectral variables $\l_j$ correspondingly. In [5] their equivalence is established in a similar way with the help of two dual $Q$-operators and their spectral descriptions analogous to (<ref>). But these spectral properties are proved differently.

Finally, we note that in the works of M. Hallnäs and S. Ruijsenaars [6], [7], [8], where only the case of real periods $\o_1, \o_2$ is considered, the wave function $\Phi_{\bl_n}(\bx_n)$ differs from ours by renormalization of spectral parameters, that is

\begin{equation}
\Phi_{\bl_n}(\bx_n) = \Psi_{ \bl_n / \o_1 \o_2 }(\bx_n).
\end{equation}

Therefore, due to homogeneity of the double sine function (<ref>) (with $\gamma = \o_1 \o_2$), the duality relation (<ref>) for this function looks as

\begin{equation}\label{HR-duality}
\Phi_{\bl_n}(\bx_n; g | \bo) = \Phi_{\bx_n}(\bl_n; g^*|\bo ).
\end{equation}

However, for the complex periods the scaling by $\o_1 \o_2$ rotates integration contours in the dual integral representation. Note that for $n=2,3$ the relation (<ref>)  was proved  by different method in [7] .

§ EXCHANGE RELATIONS

In this section we derive exchange relations between $Q$-operators and raising operators,
as well as between raising operators themselves, starting from basic commutativity property of $Q$-operators. We start with exchange relation between raising operators.

In terms of kernels (<ref>) the identity ll is equivalent to the relation

\begin{equation}
\Lambda(\bm{x}_n, \bm{z}_{n - 2}; \lambda, \rho) = e^{\const (\rho - \lambda) (\bbx_n + \bbz_{n - 2} ) } \, \Lambda(\bm{x}_n, \bm{z}_{n - 2}; \rho, \lambda),
\end{equation}

\begin{equation}
\Lambda(\bm{x}_n, \bm{z}_{n - 2}; \lambda, \rho) = \int_{\mathbb{R}^{n - 1}} d\bm{y}_{n - 1} \, \mu  (\bm{y}_{n - 1}) \, e^{ \const (\rho - \lambda) \bby_{n-1} } \, \K(\bm{x}_n, \bm{y}_{n - 1}) \, \K(\bm{y}_{n - 1}, \bm{z}_{n - 2}).
\end{equation}

Denoting $x_n := z_{n - 1}$ we arrive at $Q$-operators commutation relation (<ref>)

\begin{equation}
Q_{n - 1}(\lambda) \, Q_{n - 1}(\rho) = Q_{n - 1}(\rho) \, Q_{n - 1}(\lambda)
\end{equation}

written in terms of kernels.

The following theorem plays the crucial role in further derivations of the duality properties of the wave functions studied in this paper.

Recall $Q$-operators commutativity relation (<ref>)

\begin{equation}\label{QQcomm2}
Q_{n}(\lambda) \, Q_n(\rho) = Q_n(\rho) \, Q_n(\lambda).
\end{equation}

The main idea of the current proof is to take the limit of the identity (<ref>) written in terms of kernels

\begin{equation}\label{QQker}
Q(\bm{x}_n, \bm{z}_n; \lambda, \rho) = Q_{}(\bm{x}_n, \bm{z}_n; \rho, \lambda)
\end{equation}

as $z_n \rightarrow \infty$, where

\begin{equation}
Q(\bm{x}_n, \bm{z}_n; \lambda, \rho) = \int_{\mathbb{R}^n} d\bm{y}_n \, Q(\bm{x}_n, \bm{y}_n; \lambda) \,Q(\bm{y}_n, \bm{z}_n; \rho), \qquad x_j, z_j \in \mathbb{R}.
\end{equation}

Here we canceled constants $d_n^2(g|\bo)$ from both sides. However, to obtain an appropriate limit we should first multiply (<ref>) by the function

\begin{equation}
r(\bm{z}_n; \rho) = \exp\Bigl( \pi \hat{g} \bigl[\bbz_{n - 1}+(2-n)z_n\bigr] + {2 \pi}\imath \rho  z_n \Bigr),
\end{equation}

where we recall the notation (<ref>)

$$\hat{g} = \frac{g}{\o_1 \o_2}.$$

We note that the integral (<ref>) is absolutely convergent for $\lambda, \rho \in \mathbb{C}$ such that
| (λ- ρ) | < ν_g = ĝ, see <cit.>. So, in what follows we assume this condition for $\l,\rho$.

Consider the left-hand side of (<ref>). Multiplying by the function $r(\bm{z}_n; \rho)$, we can rewrite it as

\begin{equation}\label{IV13}
r(\bm{z}_n; \rho) \, Q_{}(\bm{x}_n, \bm{z}_n; \lambda, \rho) = \int_{\mathbb{R}^n} d\bm{y}_n F(\bm{x}_n, \bm{y}_n, \bm{z}_n; \lambda, \rho),
\end{equation}

\begin{equation}\label{F}
\begin{aligned}
F &= e^{{2 \pi \imath}  \lambda (\bbx_n - \bby_n) + 2\pi\imath\left( \rho - \frac{\imath \hat{g}}{2} \right) ( \bby_n - \bbz_{n - 1} )  } \\[10pt]
&\times \K (\bm{x}_n, \bm{y}_n) \, \mu  (\bm{y}_n) \, \K (\bm{y}_n, \bm{z}_{n - 1}) \, \mu  (\bm{z}_{n - 1}) \\[5pt]
&\times \prod_{j = 1}^n e^{ \pi\hat{g} (z_n - y_j) } \K (y_j - z_n) \; \prod_{j = 1}^{n - 1} e^{ 2\pi \hat{g} (z_j - z_n) } \mu  (z_n - z_j) \mu(z_j - z_n).
\end{aligned}
\end{equation}

The variable $z_n$ is contained only in the last line. The function $r(\bm{z}_n; \rho)$ was chosen such that products in the last line have a pointwise limit

\begin{equation}
\lim_{z_n \rightarrow \infty} \, \prod_{j = 1}^n e^{ \pi\hat{g} (z_n - y_j) } \K (y_j - z_n) \; \prod_{j = 1}^{n - 1} e^{ 2\pi \hat{g} (z_j - z_n) } \mu  (z_n - z_j) \mu(z_j - z_n) = 1.
\end{equation}

To evaluate this limit we use asymptotic formulas (<ref>)

\begin{equation}\label{IV8}
\mu(x) \sim e^{\pi \hat{g} | x | \pm \imath \frac{\pi \hat{g} g^*}{2} }, \qquad K(x) \sim e^{- \pi \hat{g} |x|}, \qquad x\rightarrow \pm \infty.
\end{equation}

And for sufficiently large $z_n $, such that $z_n > z_j$ for all $j = 1, \dots, n - 1$, it is bounded from above

\begin{equation}\label{IV15}\left|
\prod_{j = 1}^n e^{ \pi\hat{g} (z_n - y_j) } \K (y_j - z_n) \; \prod_{j = 1}^{n - 1} e^{ 2\pi \hat{g} (z_j - z_n) } \mu  (z_n - z_j) \mu(z_j - z_n) \right| \leq C_1
\end{equation}

with some $C_1$ which depends on $g$ and $\bo$. For this we use bounds (<ref>)

\begin{equation}\label{IV16}
| \mu  (x) | \leq C e^{\pi\nu_g |x| }, \qquad |\K (x)| \leq C e^{- \pi\nu_g |x| }, \qquad x \in \mathbb{R}
\end{equation}

where again we denoted
ν_g=ĝ >0 and $C$ is a constant which depends on $g,\bo$.

The next step is to derive a bound for the whole function $F$ in order to use dominated convergence theorem and evaluate $z_n \rightarrow \infty$ limit of the integral (<ref>). For all functions in the second line (<ref>) use bounds (<ref>) together with triangle inequalities

\begin{equation}\label{Kmu-bounds2}
\begin{aligned}
|\K (x_i - y_j)| &\leq C \, e^{-\pi\nu_g |x_i - y_j|} \leq C \, e^{-\pi\nu_g (|y_j| - |x_i|)}, \\[6pt]
|\mu  (y_i - y_j)| &\leq C \, e^{\pi\nu_g |y_i - y_j|} \leq C \, e^{\pi\nu_g(|y_i| + |y_j|)}.
\end{aligned}
\end{equation}

Hence, using also bound (<ref>) for the whole function $F$ we obtain

\begin{equation}
| F | \leq C_2(g, \bm{\omega}, \bm{x}_n, \bm{z}_{n - 1}) \; \exp \Bigl( \bigl[ \, | 2 \pi \Im(\lambda - \rho) + \pi\nu_g | - \pi\nu_g \bigr] \| \bm{y}_n \| \Bigr)
\end{equation}

with some $C_2$, where by $\|\bm{y}_n\|$ we mean $L^1$-norm

\begin{equation}
\| \bm{y}_n \| = \sum_{j = 1}^n |y_j|.
\end{equation}

Function from the right is integrable when
(λ- ρ) ∈(-ν_g, 0). Note that this condition is stronger than the one assumed before (<ref>). So, assuming it we use dominated convergence theorem to evaluate the limit of the integral (<ref>)

\begin{equation}
\begin{aligned}
\lim_{z_n \rightarrow \infty} \, r(\bm{z}_n; \rho) \, Q_{}(\bm{x}_n, \bm{z}_n; \lambda, \rho) &= \int_{\mathbb{R}^n} d\bm{y}_n \; e^{2 \pi \imath \left[ \lambda (\bbx_n - \bby_n) + \left( \rho - \frac{\imath \hat{g}}{2} \right) ( \bby_n - \bbz_{n - 1} ) \right] } \\[10pt]
&\times \K (\bm{x}_n, \bm{y}_n) \, \mu  (\bm{y}_n) \, \K (\bm{y}_n, \bm{z}_{n - 1}) \, \mu  (\bm{z}_{n - 1}).
\end{aligned}
\end{equation}

The result coincides with the kernel of operator $Q_n(\lambda) \, \Lambda_{n}(\rho -  \imath \hat{g}/2)$ modulo multiplication by constant $d_n(g|\bo) d_{n - 1}(g|\bo)$.

The next step is to take the limit of the right-hand side of equality (<ref>) multiplied by the function $r(\bm{z}_n; \rho)$

\begin{equation}
r(\bm{z}_n; \rho) \, Q_{}(\bm{x}_n, \bm{z}_n; \rho, \lambda) = \int_{\mathbb{R}^n} d\bm{y}_n \; G(\bm{x}_n, \bm{y}_n, \bm{z}_n; \rho, \lambda).
\end{equation}

In contrast with the function $F$ appeared on the left-hand side, the function $G$ doesn't have a pointwise limit as $z_n \rightarrow \infty$. The trick is to perform a certain shift of integration variable $y_n$ and then take the limit.

But before we also need to slightly rewrite this integral. Notice that the function $G$ is symmetric with respect to variables $y_j$. Define the domain $D_j$ as

\begin{equation}
D_j = \{ \bm{y}_n \in \mathbb{R}^n \colon y_j \geq y_k, \; \forall k \in [n] \setminus \{j\} \}.
\end{equation}

We have the equality

\begin{equation}
\bm{1}_{D_1} + \bm{1}_{D_2} + \ldots + \bm{1}_{D_n} = 1,
\end{equation}

where $\bm{1}_{D_j}$ is the indicator function of the domain $D_j$. Therefore using symmetry of $G$ we can rewrite the last integral

\begin{equation}
r(\bm{z}_n; \rho) \, Q_{}(\bm{x}_n, \bm{z}_n; \rho, \lambda) = n \int_{\mathbb{R}^n} d\bm{y}_n \; \bm{1}_{D_n} \, G(\bm{x}_n, \bm{y}_n, \bm{z}_n; \rho, \lambda),
\end{equation}

where we gained factor of $n$ due to the equality above.

Now we should shift the integration variable $y_n \rightarrow y_n + z_n$. Let us write the result after the shift

\begin{equation}\label{IV24}
\begin{aligned}
r(& \bm{z}_n; \rho) \, Q_{}(\bm{x}_n, \bm{z}_n; \rho, \lambda) = n \int_{\mathbb{R}^n} d\bm{y}_n \; e^{2 \pi \imath \left[ \left( \rho - \frac{\imath \hat{g}}{2} \right) (\bbx_n - \bby_n) + \lambda ( \bby_n - \bbz_{n - 1} ) \right] } \\[10pt]
&\times \K (\bm{x}_n, \bm{y}_{n - 1}) \, \mu  (\bm{y}_{n - 1}) \, \K (\bm{y}_{n - 1}, \bm{z}_{n - 1}) \, \K (y_n) \, \mu  (\bm{z}_{n - 1}) \, R(\bm{x}_n, \bm{y}_n, \bm{z}_n),
\end{aligned}
\end{equation}

where by $R$ we denoted all factors that contain $z_n$

\begin{equation}
\begin{aligned}
R&(\bm{x}_n, \bm{y}_n, \bm{z}_n) =  \prod_{j = 1}^n e^{ \pi \hat{g} (y_n + z_n - x_j) } \K (x_j - y_n - z_n) \prod_{j = 1}^{n - 1} e^{ 2\pi \hat{g} (y_{jn} - z_n) } \mu   (y_{nj} + z_n) \mu(y_{jn} - z_n) \\[3pt]
&\times\prod_{j = 1}^{n - 1} e^{ \pi\hat{g}(z_{nj} + y_{nj} + z_n) } \K (y_n + z_{nj}) \, \K (y_j - z_n) \; \prod_{j = 1}^{n - 1} e^{ 2 \pi \hat{g}z_{jn}}  \mu  (z_n - z_j) \mu(z_j - z_n) \cdot \bm{1}_{D_n'}
\end{aligned}
\end{equation}

Here, for brevity, we introduced notation $y_{nj} = y_n - y_j$. Domain of the indicator function has changed after the shift

\begin{equation}
D'_n = \{ \bm{y}_n \in \mathbb{R}^n \colon y_n + z_n \geq y_k, \; \forall k \in [n-1] \}.
\end{equation}

Due to the asymptotics (<ref>) it has a pointwise limit

\begin{equation}
\lim_{z_n \rightarrow \infty} R(\bm{x}_n, \bm{y}_n, \bm{z}_n) = 1.
\end{equation}

Hence, to use the dominated convergence theorem we need to write a bound independent of $z_n$ for it. Note that in the presence of the indicator function we have equality

\begin{equation}
y_{nj} + z_n = |y_{nj} + z_n|, \qquad j = 1, \ldots, n - 1.
\end{equation}

Therefore, for a factor from the second product in $R$ we can write

\begin{equation}
\left|e^{ 2\pi \hat{g} (y_{jn} - z_n) } \mu   (y_{nj} + z_n) \mu(y_{jn} - z_n)\right| \leq C^2,
\end{equation}

where we used the bound for $\mu $ (<ref>). Factors from three other products can be estimated similarly using bounds on $\K $ and $\mu$ (<ref>) and assuming  that $z_n$ is sufficiently large, such that $z_n > z_j$ for all $j = 1, \dots, n - 1$. We conclude that $R$ is bounded

\begin{equation}
|R(\bm{x}_n, \bm{y}_n, \bm{z}_n)| \leq C_3(g, \bm{\omega})
\end{equation}

with some $C_3$.

Now using this fact together with bounds (<ref>) we estimate the integrand in (<ref>)

\begin{equation}
\begin{aligned}
\bigl| \bm{1}_{D'_n} \, &G(\bm{x}_n, \bm{y}_{n - 1}, y_n + z_n, \bm{z}_n; \rho, \lambda) \bigr| \\[6pt]
&\leq C_4 \exp \Bigl( \bigl[ \, \left| 2\pi \Im(\rho - \lambda) - \pi\nu_g \right| - 3\pi\nu_g \bigr] \| \bm{y}_{n - 1} \| \\[6pt]
&\hspace{1.5cm}+ \bigl[ \, \left| 2 \pi \Im(\rho - \lambda) - \pi\nu_g \right| - \pi\nu_g \bigr] | y_n |  \Bigr)
\end{aligned}
\end{equation}

with some $C_4(g, \bm{\omega}, \bm{x}_n, \bm{z}_{n - 1})$. Clearly, function from the right doesn't depend on $z_n$ and is integrable for $\Im(\rho - \lambda) \in (0, \nu_g)$. The same condition appeared when we considered the left-hand side of the equality (<ref>). Thus, assuming it we use dominated convergence theorem to evaluate the limit $z_n \rightarrow \infty$ of the integral (<ref>)

\begin{equation}
\begin{aligned}
\lim_{z_n \rightarrow \infty} \, r( \bm{z}_n; \rho) \, Q_{}(&\bm{x}_n, \bm{z}_n; \rho, \lambda) = n \int_{\mathbb{R}^n} d\bm{y}_n \; e^{{2 \pi \imath}{} \left[ \left( \rho - \frac{\imath \hat{g}}{2} \right) (\bbx_n - \bby_n) + \lambda ( \bby_n - \bbz_{n - 1} ) \right] } \\[10pt]
&\times \K (\bm{x}_n, \bm{y}_{n - 1}) \, \mu  (\bm{y}_{n - 1}) \, \K (\bm{y}_{n - 1}, \bm{z}_{n - 1}) \, \K (y_n) \, \mu  (\bm{z}_{n - 1}).
\end{aligned}
\end{equation}

The integral over $y_n$ has separated and represents Fourier transform of function $\K $ (<ref>)

\begin{equation}
\int_{\mathbb{R}} dy_n \; e^{{2 \pi \imath}{} \left( \lambda - \rho + \frac{\imath \hat{g}}{2} \right) y_n }  \K (y_n) = \sqrt{\omega_1 \omega_2} \, S(g) \, \KK \Bigl(\lambda - \rho + \frac{\imath \hat{g}}{2} \Bigr).
\end{equation}

The rest part of the integral coincides with the kernel of operator $\Lambda_{n}(\rho - \imath \hat{g}/2) \, Q_{n - 1}(\lambda)$ up to constant $d^2_{n - 1}(g|\bo)$.

So, starting with the commutativity of $Q$-operators in the limit $z_n \rightarrow \infty$ we obtain the equality

\begin{equation}
Q_n(\lambda) \, \Lambda_{n}\Bigl(\rho - \frac{\imath \hat{g}}{2} \Bigr) = \KK \Bigl(\lambda - \rho + \frac{\imath \hat{g}}{2} \Bigr) \Lambda_{n}\Bigl(\rho - \frac{\imath \hat{g}}{2} \Bigr) Q_{n-1}(\lambda),
\end{equation}

assuming $\Im(\rho - \lambda) \in (0, \nu_g) = (0, \Re \hat{g})$. Shifting the parameter $\rho \rightarrow \rho +  \imath \hat{g}/2$ we arrive at the identity stated in the theorem (<ref>).

Exchange relation (<ref>) at the edge of convergence strip $\Im(\lambda - \rho) = \pm \nu_g/2$ could be also obtained, although it is beyond the scope of the present paper. In this case one should think of the corresponding kernels (that won't be represented by absolutely convergent integrals anymore) as generalized functions.

§ EIGENFUNCTIONS

§.§ Eigenfunctions of Macdonald operators

Here we  prove that the wave functions I12 are eigenfunctions of Macdonald operators for complex valued periods $\bo$ and coupling constant $g$ assuming the conditions I0a,  I0b and I13a. The proof differs from that of [6] by complex valued periods and the use of bilinear scalar product BM3a.
We prove I13 by induction over $n$. It is clearly correct for $n=1$. For $n>1$ under the conditions  I0a and  I0b the function $\Psi_{\bl_n}(\bx_n)$ with
| x_j | < g^*/2,     | (ł_j - ł_k) | < θ,     j,k=1,…,nwith some positive $\theta$  is defined by absolutely converging multiple integral I12 (see Proposition <ref> in Section <ref>), which can be written via recursive procedure

\begin{equation}\label{Psi-rec2a}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = \Lambda_{n}(\lambda_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}), \qquad \Psi_{\lambda_1}(x_1) = e^{{2\pi \imath}{} \lambda_1 x_1}.
\end{equation}

In the space of functions $\vf(\by_{n-1})$ analytical in a small strip around the real plane
introduce the symmetric bilinear pairing
(,ψ)= ∫_^nd_n   μ(_n)   (_n)ψ(-_n)assuming that it is defined when the corresponding integral converges.
Denote by $\tau_y$ the operator that changes the sign of argument in a function
τ_y (_n)=(-_n). Then we can rewrite this pairing as
(,ψ)= ∫_^nd_n   μ(_n)   (_n)  τ_y [ψ(_n) ].In this notation, the inductive procedure Psi-rec2a can be written as
Ψ_λ_n(x_n)=e^ł_n_n(K(_n,_n-1),e^ł_n_n-1 τ_y Ψ_λ_n-1(y_n-1)).Note the important property of the integration contour $C= \R^n$ in the integral BM3b:
it separates two series of poles of the kernel function:

y_i=x_j+g^∗/2 +mø_1+kø_2,    and    y_i= x_j-g^∗/2 -mø_1-kø_2,      m,k≥0,
and two series of poles of the measure function
y_i=y_j+g  +mø_1+kø_2,     and    y_i=y_j-g  -mø_1-kø_2,     m,k≥0.
The kernel function $K(\bx_n,\by_n)$ satisfies the relation [17]
(M_r(_n;g|)- M_r(-_n;g|))K(_n,_n)=0which is the corollary of the trigonometric version of the kernel function identity [17, 11], valid for any tuples $\bx_n$ and $\by_n$ of $n$ complex variables and arbitrary parameter $\a$:

In the limit $y_n$ tends to $\imath\infty$ the identity   I9 degenerates to

which implies the relation

M_r(_n;g|   )K(_n,_n-1)

= [M_r(-_n-1;g|)+M_r-1(-_n-1;g|) ] K(_n,_n-1).

Besides, Macdonald operators $M_r(\by_{n-1};g|\bo)$ are symmetric with respect to the bilinear pairing
BM3a,	 see <cit.>,
Acting with Macdonald operator $M_r(\bx_n;g|\bo)$ on the wave function (<ref>) we have
M_r(_n;g|)   e^ł_n_n(K(_n,_n-1),e^ł_n_n-1 τ_y Ψ_λ_n-1(y_n-1))

= e^2πł_n_n + 2πr ł_n ø_1 (M_r(_n;g|) K(_n,_n-1), e^2πł_n_n-1   τ_y Ψ__n-1(_n-1)).	Using I8a we rewrite the right hand side of BM7 as

e^2πł_n_n + 2πr ł_n ø_1 [(τ_y M_r(_n-1;g|)τ_y K(_n,_n-1), e^2πł_n_n-1   τ_y Ψ__n-1(_n-1) )

+(τ_y M_r-1(_n-1;g|)τ_y K(_n,_n-1), e^2πł_n_n-1   τ_y Ψ__n-1(_n-1) )]

=     e^2πł_n_n + 2πr ł_n ø_1 [(M_r(_n-1;g|)τ_y K(_n,_n-1), e^-2πł_n_n-1   Ψ__n-1(_n-1) )

+(M_r-1(_n-1;g|)τ_y K(_n,_n-1), e^-2πł_n_n-1    Ψ__n-1(_n-1) )].
Then using the symmetricity of Macdonald operators BM6 we rewrite BM7a as

e^2πł_n_n + 2πr ł_n ø_1 [(τ_y K(_n,_n-1),  M_r(_n-1;g|)e^-2πł_n_n-1   Ψ__n-1(_n-1)  )
+ ( τ_y K(_n,_n-1), M_r-1(_n-1;g|)e^-2πł_n_n-1   Ψ__n-1(_n-1)  )]
under the assumption that the shifts of the integration contour implemented in the equality BM6 do not spoil separation of the poles t04 and t04a property. For t04 this happens when
ø_1<g^∗,which is equivalent to (<ref>). The conditions t04a could be violated in the right hand side of the adjointness relation. However, in this case the nearest poles $\imath y_a=\imath y_b+g$ are canceled by zeroes of the coefficients $\sh\frac{\pi}{\o_2} (y_a-y_b-\imath g)$ of Macdonald operators, see <cit.>, and we do not get here additional restrictions.
Finally we rewrite  BM5 as

e^2πł_n_n [ ( τ_y K(_n,_n-1),  e^-2πł_n_n-1  M_r(_n-1;g|) Ψ__n-1(_n-1)  )

+( τ_y K(_n,_n-1),e^-2πł_n_n-1+2πł_nø_1 M_r-1(_n-1;g|)   Ψ__n-1(_n-1) )]

=    e^2πł_n_n [ (  K(_n,_n-1),  e^2πł_n_n-1  τ_y M_r(_n-1;g|) Ψ__n-1(_n-1)  )

+( K(_n,_n-1),e^2πł_n_n-1+2πł_nø_1 τ_y M_r-1(_n-1;g|)   Ψ__n-1(_n-1) )].
By induction assumption,
M_s(_n-1;g|)Ψ__n-1(_n-1) =e_s(e^2πł_1ø_1,…, e^2πł_n-1ø_1) Ψ__n-1(_n-1).Then BM5c implies the relation
M_r(_n;g|)Ψ__n(_n) =

[e_r(e^2πł_1ø_1,…, e^2πł_n-1ø_1)+ e^ł_nø_1e_r-1(e^2πł_1ø_1,…, e^2πł_n-1ø_1)]Ψ__n(_n)
which in turn implies the theorem statement I13 due to the relation
e_r(z_1,…, z_n)=e_r(z_1,…, z_n-1)+z_ne_r-1(z_1,…, z_n-1).

§.§ Eigenfunctions of $Q$-operators

In this section we show that the iterated integrals I12 represent as well joint eigenfunctions of the commutative family of operators $Q_n(\l)$. Modulo convergence of the appearing integrals, it is a formal consequence of Theorem <ref>. However, these convergence details are not trivial.
Proof of this theorem follows from the iterative representation of the wave function

\begin{equation}\label{Psi-iter}
\Psi_{\bl_n}(\bx_n; g)=\LL{n}(\l_n) \, \LL{n-1}(\l_{n-1})\cdots \LL{2}(\l_2) \, e^{\const \l_1x_1}
\end{equation}

and the exchange relation given by Theorem <ref>

\begin{equation}
Q_n(\lambda) \, \Lambda_{n}(\rho) = \KK(\lambda - \rho) \; \Lambda_{n}(\rho) \, Q_{n-1}(\lambda).
\end{equation}

The subtle point is the possibility to switch integration order inside of multiple integrals appearing during this calculation. It is indeed possible due to the absolute convergence of these integrals, which we prove in the following proposition. After it we give a proof of Theorem <ref>.

Multiple integrals

\begin{align}\label{int1}
I_{\lambda, \bm{\lambda}_n} &= Q_n(\lambda) \, \LL{n}(\l_n) \, \LL{n-1}(\l_{n-1})\cdots \LL{2}(\l_2) \, e^{\const \l_1x_1}, \\[7pt]
\label{int2}
J_{\lambda, \bm{\lambda}_n} &= \Lambda_{n}(\lambda_n) \, Q_{n - 1}(\lambda) \, \LL{n-1}(\l_{n-1})\cdots \LL{2}(\l_2) \, e^{\const \l_1x_1}, \\[7pt]
\label{int3}
\Psi_{\bl_n} &= \LL{n}(\l_n) \, \LL{n-1}(\l_{n-1})\cdots \LL{2}(\l_2) \, e^{\const \l_1x_1}.
\end{align}

are absolutely convergent for $\lambda, \lambda_j \in \mathbb{C}$ such that

\begin{equation}
|\Im (\l-\l_n)|<\frac{\nu_g}{2}(1-\ve),\qquad
|  \Im(\lambda_k - \lambda_j)| < \theta(\ve), \qquad k,j = 1, \dots, n
\end{equation}

for any $\ve \in [0,1)$ and
θ() = ν_g/4 (n - 1)! e   .
In particular, the wave function $\Psi_{\bl_n}(\bx_n)$ is given by a multiple integral, converging when all the variables are in small strips around the real lines.
We present the proof for $n\geq 2$. For $n=1$ it is an elementary check.

Consider the first integral (<ref>). It has $n$ levels of integration variables. Denoting integration variables of $k$-th level by

\begin{equation}
\bm{y}_{k} = \bigl( y_1^{(k)} , \dots, y_k^{(k)} \bigl),
\end{equation}

we write integral in its full form

\begin{equation}
\begin{aligned}
I_{\lambda, \bm{\lambda}_n} (\bm{x}_n) = C_I \int_{\mathbb{R}^{n(n+1)/2} } d\bm{y}_n & d\bm{y}_{n - 1} \dots d\bm{y}_1 \; e^{ {2\pi \imath}{} \left[ \lambda \bbx_n + (\lambda_n - \lambda) \bby_{n} \right] } \, \K  (\bm{x}_n, \bm{y}_{n}) \\[5pt]
& \hspace{0.6cm} \times  \prod_{k = 2}^n e^{{2\pi \imath}{}(\lambda_{k - 1} - \lambda_{k}) \bby_{k - 1}} \, \mu ( \bm{y}_{k} ) \, \K  (\bm{y}_{k}, \bm{y}_{k - 1} ).
\end{aligned}
\end{equation}

Here $C_I$ contains all constants $d_k$ from kernels of operators. Denote the integrand as $F$.
Assume that
|(ł-ł_n)| ≤δ_Qν_g/2,    |(ł_k - 1-ł_k)| ≤δ_Λν_g/2,    2<k≤nwhere $\delta_Q$ and $\delta_{\Lambda}$ are some positive constants. Using bounds (<ref>)

\begin{equation}\label{Kmu-bounds}
| \mu  (x) | \leq C e^{ \pi\nu_g |x| }, \qquad |\K (x)| \leq C e^{ - \pi\nu_g |x| }, \qquad x \in \mathbb{R},
\end{equation}

together with triangle inequalities $|x_i - y_j| \geq |y_j| - |x_i|$ we get

\begin{equation}\label{int1b}
| F | \leq C_1 \exp \pi\nu_g\biggl( \bigl[\delta_Q - n   \bigr]  \|\bm{y}_n \| +  S_n (\bm{y}_1, \dots,  \bm{y}_n ) +\delta_{\Lambda}\sum_{k=1}^{n-1}\|\bm{y}_k\|\biggr),
\end{equation}

with some $C_1(g, \bm{\omega}, \bm{x}_n)$ and the function $S_n$ in exponent is defined by recurrence relation

\begin{equation}\label{Srec}
\begin{aligned}
S_n (\bm{y}_1, \dots, \bm{y}_n ) = \sum_{\substack{i, j = 1 \\ i \not= j}}^n \bigl| y_i^{(n)} - y_j^{(n)} \bigr| &- \sum_{i = 1}^{n} \sum_{j = 1}^{n - 1} \, \bigl| y_i^{(n)} - y_j^{(n - 1)} \bigr| \\
&+ S_{n - 1} (\bm{y}_1, \dots,  \bm{y}_{n - 1} )
\end{aligned}
\end{equation}

with $S_1 = 0$. As before, $\| \bm{y}_k \|$ denotes $L^1$-norm.  In appendix we prove the inequality (<ref>)

\begin{equation}
S_n \leq (n - 1 +  \epsilon) \| \bm{y}_n \| - \frac{\ve}{c_n} \, \sum_{k = 1}^{n - 1} \| \bm{y}_k \|
\end{equation}

for any $\epsilon \in \left[ 0,2(n - 1) \right]$, where $c_n$ is defined by the recurrence relation
c_n = (n - 1)(c_n-1 + 1)    with    c_1 = 0.Clearly,
(n-1)! ≤c_n<(n-1)!e,     n≥2.Inserting it in (<ref>) we obtain

\begin{equation}\label{kh3}
|F| \leq C_1 \exp \pi\nu_g \biggl( \Bigl[ ( -1+\delta_Q + \ve)  \Bigr]  \|\bm{y}_n \|  +\left(\delta_{\Lambda}-\frac{\ve}{c_n}\right)  \sum_{k = 1}^{n - 1} \| \bm{y}_k \| \biggr).
\end{equation}

We have an exponential decay in the right hand side of kh3, if
δ_Q<1-,    δ_Λ</(n-1)!e</c_n.
This proves the statement of the proposition for the integral int1.

Next, consider the second integral (<ref>). It also has $n$ levels of integration variables. We will denote by $\bm{t}_{n - 1} = (t_1, \dots, t_{n - 1})$ variables associated with the operator $\Lambda_{n}$ and the rest levels of variables, as before, will be denoted by $\bm{y}_{k}$ ($k = 1, \dots, n - 1$). Full form of the integral

\begin{align} \nonumber
J_{\lambda, \bm{\lambda}_n} (\bm{x}_n) &= C_J \int  d\bm{t}_{n - 1} d\bm{y}_{n - 1} \ldots d\bm{y}_1 \, e^{ {2\pi \imath}{} \left[ \lambda_n \bbx_n + (\lambda - \lambda_n) \bbt_{n - 1} + (\lambda_{n - 1} - \lambda) \bby_{n - 1} \right] } \, \K (\bm{x}_n, \bm{t}_{n - 1})  \\[4pt]
& \times \mu  (\bm{t}_{n - 1}) \, \K (\bm{t}_{n - 1}, \bm{y}_{n - 1} ) \prod_{k = 2}^{n - 1} e^{{2\pi \imath}{}  (\lambda_{k - 1} - \lambda_{k}) \bby_{k - 1}} \, \mu   ( \bm{y}_{k}) \, \K  (\bm{y}_{k}, \bm{y}_{k - 1}).
\end{align}

Again, $C_J$ contains all unimportant constants. The integration goes over $\mathbb{R}^{(n - 1)(n + 2)/2}$. Denote the integrand as $G$. Assuming the conditions kh1, using the bounds (<ref>) and triangle inequalities we arrive at

\begin{equation}\label{int2b}
| G | \leq  C_2 \exp \pi\nu_g\biggl(\delta_Q \, \left| \sum\limits_{j = 1}^{n - 1} \bigl( t_j - y_j^{(n-1)} \bigr) \right|+\delta_\Lambda \sum_{k = 1}^{n - 1} \| \bm{y}_k \|- n   \| \bm{t}_{n - 1} \| +  T_{n - 1} \biggr) ,
\end{equation}

with some $C_2(g, \bm{\omega}, \bm{x}_n)$, where the function $T_{n - 1}$ is defined as

\begin{equation}
\begin{aligned}
T_{n - 1}(\bm{y}_1, \dots, \bm{y}_{n - 1}, \bm{t}_{n - 1}) = \sum_{\substack{i, j = 1 \\ i \not= j}}^{n - 1} | t_i - t_j | &- \sum_{i, j = 1}^{n - 1} \bigl| t_i - y_j^{(n - 1)} \bigr| \\[4pt]
&+ S_{n - 1}(\bm{y}_1, \dots, \bm{y}_{n - 1}).
\end{aligned}
\end{equation}

In appendix we prove the inequality (<ref>) (here we replace $r$ by $1-\ve$)

\begin{equation}\label{kh6}
\begin{aligned}
T_{n - 1} \leq ( n -  \ve ) \| \bm{t}_{n - 1} \| %&
- \frac{ \ve}{2 (n - 1) c_{n - 1}} \sum_{k = 1}^{n - 1} \| \bm{y}_k \|
- (1-\ve) \, \biggl| \sum_{j = 1}^{n - 1} \bigl( t_j - y_j^{(n - 1)} \bigr) \biggr|
\end{aligned}
\end{equation}

for any $\ve \in [0, 1]$. Using it we have the bound

| G | ≤C_2 expπν_g(   - t_n - 1     +(δ_Λ-  /2 (n - 1) c_n - 1) ∑_k = 1^n - 1  y_k

+ (δ_Q-1+ ) | ∑_j = 1^n - 1 ( t_j - y_j^(n-1) ) |   ) .
We thus see that the integral (<ref>) converges once

δ_Q<1-,     δ_Λ<  /2 (n-1)! e   <  /2 (n - 1) c_n - 1.

Finally, the third integral (<ref>) in its full form

\begin{equation}
\begin{aligned}
\Psi_{\bm{\lambda}_n} (\bm{x}_n) = C_\Psi \int_{\mathbb{R}^{n(n-1)/2} } & d\bm{y}_{n - 1} \dots d\bm{y}_1 \; e^{ {2\pi \imath}{} \left[ \lambda_n \bbx_n + (\lambda_{n - 1} - \lambda_{n}) \bby_{n - 1} \right] } \K (\bm{x}_n, \bm{y}_{n - 1}) \\[5pt]
& \hspace{0.6cm} \times  \prod_{k = 2}^{n - 1} e^{{2\pi \imath}{}(\lambda_{k - 1} - \lambda_{k}) \bby_{k - 1}} \, \mu   ( \bm{y}_{k} ) \, \K  (\bm{y}_{k}, \bm{y}_{k - 1} )
\end{aligned}
\end{equation}

almost coincides with the first integral $I_{\lambda_n, \bm{\lambda}_{n - 1}}$ up to additional functions

\begin{equation}
e^{{2\pi \imath}{}  \lambda_n x_n} \,  \prod_{j = 1}^{n - 1} \K (x_n - y_j)
\end{equation}

that only improve convergence.
Its integrand, denote it by $H$, can be bounded in a similar way as

\begin{equation}\label{kh8}
|H| \leq C_3 \exp \pi\nu_g \biggl( \Bigl[ ( -2+\delta_\Lambda + \ve)  \Bigr]  \|\bm{y}_n \|  +\left(\delta_{\Lambda}-\frac{\ve}{c_n}\right)  \sum_{k = 1}^{n - 1} \| \bm{y}_k \| \biggr).
\end{equation}

It converges if
δ_Λ<2-,    δ_Λ</(n-1)!e
for some $\ve >0$. Assuming $\ve \in [0, 1)$ the first condition follows from the second one.

So, all three integrals converge under assumptions

\begin{equation}
|\Im(\l - \l_n) | < \frac{\nu_g}{2}(1 - \ve), \qquad |\Im(\l_j - \l_k)|< \frac{\nu_g}{2} \cdot  \frac{\ve}{2(n - 1)! e}
\end{equation}

with any $\ve \in [0,1)$.

Since all constants that appear in the bounds are uniform on compact sets of parameters $\bo, g$ preserving the conditions (<ref>), (<ref>), the considered multiple integrals are analytic functions of the parameters $\bo, g$, as well as of the variables $x_i$ and $\l_j$. In the case of real periods $\o_1, \o_2$ convergence of the third integral (<ref>), that is the wave function representation, was first proved in [6].

In the notation of the previous proposition the theorem states that

\begin{equation}
I_{\lambda, \bm{\lambda}_n} = \prod_{j = 1}^n \KK (\lambda - \lambda_j) \, \Psi_{\bm{\lambda}_n}.
\end{equation}

The proof goes by induction. In the case $n = 1$ this integral

\begin{equation}
I_{\lambda, \lambda_1} = d_1(g|\bo) \int_{\mathbb{R}} dy_1 \, e^{ {2\pi \imath}{} \left[ \lambda x_1 + (\lambda_1 - \lambda) y_1  \right] } \K (x_1 - y_1) = \KK(\lambda - \lambda_1) \, e^{ {2\pi \imath}{} \lambda_1 x_1 }
\end{equation}

coincides with the Fourier transform of the function $\K $ (<ref>). Next, consider the case of arbitrary $n$. The main ingredient of the proof is the local relation (<ref>)

\begin{equation}\label{QL}
Q_n(\lambda) \, \Lambda_{n}(\rho) = \KK(\lambda - \rho) \;\Lambda_{n}(\rho) \, Q_{n-1}(\lambda).
\end{equation}

Consider the function in question

\begin{equation}
I_{\lambda, \bm{\lambda}_n} = Q_n (\lambda) \, \Psi_{\bm{\lambda}_n} = Q_n(\lambda) \,\Lambda_{n} (\lambda_n) \cdots \Lambda_{2}(\lambda_2) \,  \Psi_{\lambda_1}.
\end{equation}

This integral has $n$ levels of integration variables that we denoted by

\begin{equation}
\bm{y}_k = \bigl(y_1^{(k)}, \dots, y_k^{(k)} \bigr), \qquad k = 1, \dots, n.
\end{equation}

By Proposition <ref> the integral is absolutely convergent. Therefore, we can use Fubini's theorem and evaluate integrals in any order. Take integral over $\bm{y}_n$ to be the first in this order. Then we use the local relation (<ref>) and obtain

\begin{equation}
\begin{aligned}
I_{\lambda, \bm{\lambda}_n} &= \KK(\lambda - \lambda_n) \;\Lambda_{n}(\lambda_n) \, Q_{n - 1}(\lambda) \, \Lambda_{n - 1}(\lambda_{n - 1}) \cdots  \Lambda_{2}(\lambda_2) \, \Psi_{\lambda_1} \\[4pt]
&= \KK(\lambda - \lambda_n) \; J_{\lambda, \bm{\lambda}_n}.
\end{aligned}
\end{equation}

Obtained integral coincides with the function $J_{\lambda, \bm{\lambda}_n}$ from Proposition <ref>. This integral has $n$ levels of integration variables that we denoted by $\bm{t}_{n - 1}$ (variables associated with the operator $\Lambda_{n}$ from the left) and $\bm{y}_{k}$ (the rest ones, $k = 1, \dots, n - 1$). Since by proposition it is also absolutely convergent we can change order of integrals and make the integral over $\bm{t}_{n - 1}$ to be the last one. Then the multiple integral we should take before it

\begin{equation}
I_{\lambda, \bm{\lambda}_{n - 1}} = Q_{n - 1}(\lambda) \, \Lambda_{n - 1}(\lambda_{n - 1}) \cdots  \Lambda_{2}(\lambda_2) \, \Psi_{\lambda_1}
\end{equation}

represents the function on the previous step of induction. So, to conclude the proof we use induction assumption and obtain desired formula.

§ DUALITY

In this section, we prove the equivalence of two integral representations of the wave function
$\Psi_{\bm{\lambda}_{n}}(\bm{x}_{n})$ encoded in duality relation duality. This will be proved in two steps. First, using the spectral properties of operators $Q_n(\l)$ and $\QQ_n(x)$ proved in the previous section and the bound on the wave function from [8], we establish the relation duality for the case of real positive periods $\o_1,\o_2$. Then using analiticity of the wave function with respect to the periods (see Remark <ref> in the previous section) we extend this result to complex periods.

We note that functions from both sides exist under assumptions (<ref>)

\begin{equation}
0< \Re g < \Re \o_1 + \Re \o_2, \qquad 0 < \Re \hat{g}^* < \Re \hat{\o}_1 + \Re \hat{\o}_2.
\end{equation}
We also assume $x_j, \l_j \in \R$; the duality relation for complex variables then follows from the analyticity arguments.  To prove this theorem, firstly, note that kernels of the $Q$-operator (<ref>) and raising $\Lambda$-operator (<ref>) are connected by the formula

\begin{equation}
\Lambda(\bm{x}_n, \bm{y}_{n - 1}; \lambda) = e^{ {2 \pi \imath}{} \lambda x_n } \, Q(\bm{x}_{n - 1}, \bm{y}_{n - 1}; \lambda) \, \prod_{j = 1}^{n - 1} \K (x_n - y_j).
\end{equation}

Therefore, the operators themselves are connected in the following way

\begin{equation}\label{LQ}
\Lambda_{n}(\lambda) = e^{{2 \pi \imath}{} \lambda x_n } \, Q_{n-1}(\lambda) \, \prod_{j = 1}^{n - 1} \K (x_n - x_j).
\end{equation}

Secondly, the function 	$\hat{\Psi}_{\bm{x}_n}(\bm{\lambda}_n)=\Psi_{\bm{x}_n}(\bm{\lambda}_n; \gg|\hat{\bo})	$ from the right-hand side of duality relation (<ref>) admits an analogous integral representation

\begin{equation}
\hat{\Psi}_{\bm{x}_n}(\bm{\lambda}_n) = \LLL_n(x_n) \, \hat{\Psi}_{\bm{x}_{n - 1}}(\bm{\lambda}_{n - 1}),
\end{equation}

iterative with respect to $\lambda_j$ variables

\begin{equation}\label{d1}
\bigl(\LLL_n(x) f\bigr)(\bm{\lambda}_n) = d_{n - 1}(\hat{g}^*|\hat{\bo}) \int_{\mathbb{R}^{n - 1}} d\bm{\gamma}_{n - 1} \, \LLL(\bm{\lambda}_n, \bm{\gamma}_{n - 1}; x) f(\bm{\gamma}_{n - 1})
\end{equation}

where the kernel

\begin{equation}
\LLL(\bm{\lambda}_n, \bm{\gamma}_{n - 1}; x) = e^{\const x (\bbl_n - \bbg_{n - 1})} \KK(\bl_n, \bg_{n - 1}) \hat{\mu}(\bg_{n - 1}).
\end{equation}

It is also an eigenfunction of the $Q$-operator acting on functions of $\lambda_j$

\begin{equation}\label{d2}
\bigl(\QQ_n(x) f\bigr)(\bm{\lambda}_n) = d_{n}(\hat{g}^*|\hat{\bo}) \int_{\mathbb{R}^n} d\bm{\gamma}_{n} \, \QQ(\bm{\lambda}_n, \bm{\gamma}_{n}; x) f(\bm{\gamma}_{n})
\end{equation}

with the kernel

\begin{equation}
\hat{Q}(\bm{\lambda}_n, \bm{\gamma}_{n}; x) = e^{\const x (\bbl_n - \bbg_{n})} \KK(\bl_n, \bg_{n}) \hat{\mu}(\bg_{n})
\end{equation}

and the corresponding eigenvalue

\begin{equation}\label{Q'Psi}
\QQ_{n}(x) \, \hat{\Psi}_{ \bm{x}_n }(\bm{\lambda}_n) = \prod_{j = 1}^n \K(x-x_j) \, \hat{\Psi}_{ \bm{x}_n }(\bm{\lambda}_n ).
\end{equation}

All of this is just a matter of renaming the variables $\lambda_j \leftrightarrow x_j$ and parameters $\bo\leftrightarrow \hat{\bo}$, $g \leftrightarrow \gg$, see the definitions of the functions $\KK$ (<ref>) and $\hat{\mu}$ (<ref>). The  operators d1 and d2 are also connected by the formula analogous to (<ref>)

\begin{equation}\label{LQ'}
\LLL_{n}(x) = e^{ {2 \pi \imath}{} x \lambda_n } \, \QQ_{n - 1}(x) \, \prod_{j = 1}^{n - 1} \KK (\lambda_n - \lambda_j).
\end{equation}

Thirdly, in the proof we use the bound derived in <cit.>. This bound holds when $x_j \in \mathbb{C}$ and $\Re g \in (0, \max(\omega_1, \omega_2)]$, but we only need its specialization for the case of real $x_j$. Denote

\begin{equation}\label{mu'}
\mu'(\bx_n) = \prod_{\substack{i,j=1 \\ i< j}}^n \mu(x_i - x_j), \qquad \hat{\mu}'(\bl_n) = \prod_{\substack{i,j=1 \\ i< j}}^n \hat{\mu}(\l_i - \l_j),
\end{equation}

so that comparing with (<ref>) and (<ref>) we have

\begin{equation}\label{mumu'}
\mu(\bm{x}_n) = \mu'(\bx_n) \mu'(-\bx_n), \qquad  \hat{\mu}(\bl_n) = \hat{\mu}'(\bl_n) \hat{\mu}'(-\bl_n).
\end{equation}

<cit.> Assuming $\Re g \in (0, \max(\omega_1, \omega_2)]$, the wave function admits the bound

\begin{equation}\label{Psi-HRbound}
\Bigl| \mu'(\bm{x}_n)\, \hat{\mu}'(\bm{\lambda}_n) \, \Psi_{ \bm{\lambda}_n }(\bm{x}_n) \Bigr| \leq P(\bm{x}_n; g| \bo),
\end{equation}

where $P$ is a polynomial with $\deg P \leq n(n - 1)/2$.

At first, assume $\o_1, \o_2$ are real positive. We prove the theorem by induction. The case $n = 1$

\begin{equation}
\Psi_{\lambda_1}(x_1) = \hat{\Psi}_{x_1}(\lambda_1) = e^{ {2 \pi \imath}{} \lambda_1 x_1 }
\end{equation}

is obvious. Suppose we proved the relation for $n - 1$ case

\begin{equation}\label{Psi-ind}
\Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}) = \hat{\Psi}_{\bm{x}_{n - 1}}(\bm{\lambda}_{n - 1}).
\end{equation}

Using formula (<ref>) we write the wave function recursively with the help of $Q$-operator

\begin{equation}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = e^{ {2 \pi \imath}{} \lambda_n x_n } \, Q_{n-1}(\lambda_n) \, \prod_{j = 1}^{n - 1} \K (x_n - x_j) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}).
\end{equation}

By induction assumption (<ref>) we can represent $\K $ functions here as action of the operator $\QQ_{n-1}(x_n)$, see eigenvalue identity (<ref>),

\begin{equation}\label{PsiQQ'}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = e^{ {2 \pi \imath}{} \lambda_n x_n } \, Q_{n-1}(\lambda_n) \, \QQ_{n - 1}(x_n) \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}).
\end{equation}

Since operators $Q_{n-1}(\lambda_n)$ and $\QQ_{n-1}(x_n)$ act on different types of variables ($\bm{x}_{n - 1}$ and $\bm{\lambda}_{n - 1}$ correspondingly), they commute: switching their order is equivalent to switching order of the corresponding integrals. Let us postpone convergence arguments to the end of this proof and assume that we can switch their order. Then we arrive at

\begin{equation}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = e^{ {2 \pi \imath}{} \lambda_n x_n } \, \QQ_{n - 1}(x_n) \, Q_{n-1}(\lambda_n)  \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}).
\end{equation}

Next step is to use the eigenvalue identity (<ref>) for the operator $Q_{n-1}(\lambda_n)$

\begin{equation}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = e^{ {2 \pi \imath}{} \lambda_n x_n } \, \QQ_{n - 1}(x_n) \, \prod_{j =1}^{n - 1} \KK (\lambda_n - \lambda_j)  \, \Psi_{\bm{\lambda}_{n - 1}}(\bm{x}_{n - 1}).
\end{equation}

Finally, using induction assumption (<ref>) and formula (<ref>) we obtain function $\hat{\Psi}_{\bm{x}_n}(\bm{\lambda}_n)$ from the right.

It is left to prove that we can switch integration order in representation (<ref>). Let us write it explicitly

\begin{equation}
\begin{aligned}
\Psi_{\bm{\lambda}_n}(\bm{x}_n) = e^{ {2 \pi \imath}{} \lambda_n x_n } \int_{\mathbb{R}^{n - 1}} d\bm{y}_{n - 1} \, \mu  (\bm{y}_{n - 1}) \, e^{ {2 \pi \imath}{} \lambda_n (\bm{x}_{n - 1} - \bm{y}_{n - 1} ) } \, \K (\bm{x}_{n - 1}, \bm{y}_{n - 1}) \\[5pt]
\times \int_{\mathbb{R}^{n - 1}} d\bm{\gamma}_{n - 1} \, \hat{\mu}_{}(\bm{\gamma}_{n - 1}) \, e^{ {2 \pi \imath}{} x_n (\bm{\lambda}_{n - 1} - \bm{\gamma}_{n - 1} ) } \, \KK (\bm{\lambda}_{n - 1}, \bm{\gamma}_{n - 1}) \; \Psi_{\bm{\gamma}_{n - 1}} (\bm{y}_{n - 1}).
\end{aligned}
\end{equation}

To switch order of integrals over $\bm{y}_{n - 1}$ and $\bm{\gamma}_{n - 1}$ we bound integrand and use Fubini's theorem. First, factor measure functions

\begin{equation}
\mu(\by_{n - 1}) = \mu'(\by_{n - 1}) \mu'(-\by_{n - 1}), \qquad \hat{\mu}(\bg_{n - 1}) = \hat{\mu}'(\bg_{n - 1}) \hat{\mu}'(-\bg_{n - 1}),
\end{equation}

where $\mu'$ and $\hat{\mu}'$ are defined in (<ref>). Then using bounds (<ref>) and analogous ones for the dual functions

\begin{align}
&| \mu  (y) | \leq C e^{ \pi\nu_g |y| }, \qquad \; \, |\K (y) |\leq C e^{ - \pi\nu_g|y| }, \qquad \;\, y \in \mathbb{R},\\[5pt]
&| \hat{\mu}  (\g) | \leq C e^{ \pi \nu_{\hat{g}^*} |\g| },  \qquad |\KK (\g) |\leq C e^{ - \pi\nu_{\hat{g}^*} |\g| }, \qquad \g \in \mathbb{R},
\end{align}

ν_g = ĝ>0,     ν_ĝ^* = g^*>0,
together with triangle inequalities

\begin{equation}
|x_i - y_j | \geq |y_j| - |x_i|, \qquad |\lambda_i - \gamma_j| \geq |\gamma_j| - |\lambda_i|
\end{equation}

we bound part of the integrand by exponentially decreasing function

\begin{equation}\label{QQ'bound}
\begin{aligned}
\Bigl| \mu' (-\bm{y}_{n - 1}) \, \hat{\mu}'(-\bm{\gamma}_{n - 1}) &\, \K (\bm{x}_{n - 1}, \bm{y}_{n - 1}) \, \KK (\bm{\lambda}_{n - 1}, \bm{\gamma}_{n - 1}) \Bigr| \\[5pt]
&\leq C_q \exp \Bigl( - \pi\nu_g \| \bm{y}_{n - 1} \| -  \pi\nu_{\hat{g}^*} \| \bm{\gamma}_{n - 1} \| \Bigr)
\end{aligned}
\end{equation}

with some $C_q(g, \bo, \bm{x}_{n - 1}, \bm{\lambda}_{n - 1})$.

For the rest part of the integrand (besides surely bounded oscillating exponents) we use estimate (<ref>). If $\Re g\in (0, \max(\omega_1, \omega_2)]$ we use it, as it is written. If $\Re g \in (\max(\omega_1, \omega_2), \omega_1 + \omega_2)$, then $\gg \in (0, \min(\hat{\o}_1,\hat{\o}_2))$. In this case we use an analogous to (<ref>) estimate for the dual wave function $\hat{\Psi}_{ \bm{y}_{n - 1} }(\bm{\gamma}_{n - 1})$

\begin{equation}
\Bigl| \mu'  (\bm{y}_{n - 1})\, \hat{\mu}'(\bm{\gamma}_{n - 1})  \, \hat{\Psi}_{ \bm{y}_{n - 1} }(\bm{\gamma}_{n - 1}) \Bigr| \leq P(\bm{\gamma}_{n - 1}; \hat{g}^*|\hat{\bo}),
\end{equation}

since wave functions are equal by induction assumption (<ref>). In both cases exponent in (<ref>) suppresses the polynomial growth. Thus, the whole integrand is absolutely integrable. This gives the proof of Theorem <ref> for real positive periods.

Since the estimate for multiple integral int3 presenting the wave function in Proposition <ref> is uniform on compact subsets of the periods, the wave function $\Psi_{\bl_n}(\bx_n)$ depends analitically on the periods $\bo$. The same is true for the function  $\hat{\Psi}_{\bx_n}(\bl_n)$. Hence, their coincidence for real periods implies their coincidence for complex periods. This ends the proof of Theorem <ref>.

Note that during the proof  we established the third, mixed integral representation of the wave function.

§ ACKNOWLEDGMENTS

The work of N. Belousov and S. Derkachov (Sections 1, 3, Appendices A, B) was supported by the Theoretical Physics and Mathematics Advancement Foundation BASIS.
The  work of S. Kharchev (Section 2) was supported by the Russian Science Foundation (Grant No. 23-41-00049). The work of S. Khoroshkin (Section 4) was supported by the International Laboratory of Cluster Geometry of National Research University Higher School of Economics, Russian Federation Government grant, ag. No. 075-15-2021-608 dated 08.06.2021.

§ APPENDIX

§ THE DOUBLE SINE FUNCTION

The  double sine  function $S_2(z):=S_2(z|\bo)$, see [12] and references therein, is a meromorphic function that satisfies two functional relations
S_2(z)/S_2(z+ø_1)=2sinπz/ø_2,    S_2(z)/S_2(z+ø_2)=2sinπz/ø_1
and inversion relation
S_2(z)S_2(-z)=-4sinπz/ø_1sinπz/ø_2,or equivalently
S_2(z)S_2(ø_1+ø_2-z)=1. The function $S_2(z)$ has poles at the points

z = m ø_1 + kø_2,     m,k≥1
and zeros at

z=-mø_1-kø_2,    m,k≥0.
For $\o_1 / \o_2 \not\in \mathbb{Q}$ all poles and zeros are simple.
In the analytic region $ \Re z \in ( 0, \Re(\omega_1 + \omega_2) )$ we have the following integral representation for the logarithm of $S_2(z)$

\begin{equation}\label{S2-int}
\ln S_2 (z) = \int_0^\infty \frac{dt}{2t} \left( \frac{\sh \left[ (2z - \omega_1 - \omega_2)t \right]}{ \sh (\omega_1 t) \sh (\omega_2 t) } - \frac{ 2z - \omega_1 - \omega_2 }{ \omega_1 \omega_ 2 t } \right).
\end{equation}

It is clear from this representation that the double sine function is homogeneous

S_2( γz | γø_1, γø_2 ) = S_2(z|ø_1, ø_2),     γ∈(0, ∞)
and invariant under permutation of periods

S_2(z| ø_1, ø_2) = S_2(z | ø_2, ø_1).
The double sine function can be expressed through the Barnes double Gamma function $\Gamma_2(z|\bo)$ [1],
and its properties follow from the corresponding properties of the double Gamma function.
It is also connected to the Ruijsenaars hyperbolic Gamma function $G(z|\bo)$ [16]

G(z|) = S_2(z + ø_1 + ø_2/2  |  )
and to the Faddeev quantum dilogarithm $\gamma(z|\bo)$ [3]
γ(z|) = S_2(-z + ø_1+ø_2/2  |  ) exp( π/2ø_1 ø_2 [z^2 + ø_1^2+ø_2^2/12 ]).
Both $G(z|\bo)$ and $\gamma(z|\bo)$ were investigated independently.

In the paper we deal only with ratios of double sine functions denoted by $\mu(x)$ (<ref>) and $K(x)$ (<ref>)

μ(x)    =S_2(x)S_2^-1 (x+g),

K(x)    =  S_2(x+ø_1+ø_2/2+g/2)S_2^-1(x+ø_1+ø_2/2-g/2).

Now we will give the key asymptotic formulas and bounds for them, which were derived in <cit.> from the known results for the double Gamma function. In what follows we assume conditions (<ref>), (<ref>)

\begin{equation}
\Re \o_j > 0, \qquad 0 < \Re g < \Re\o_1 + \Re \o_2, \qquad \Re \hat{g} > 0,
\end{equation}

where we denoted

\begin{equation}
\hat{g} = \frac{g}{\o_1 \o_2}.
\end{equation}

Let $\sigma_i$ be the arguments of the periods $\o_i$, $|\sigma_i|<\pi/2$. Because of the symmetry (<ref>), we may assume that $\sigma_1 \geq \sigma_2$. Let $D_+$ and $D_-$ be the cones of poles (<ref>) and zeros (<ref>) of $S_2(z)$:
D_+={ zσ_2< z<σ_1},    D_-={ zπ+σ_2< z<π+σ_1},    D=D_+∪D_-.Denote by $d(z,D)$ the distance between a point $z$ and the cones $D$. Using Barnes' Stirling formula for the asymptotic of double Gamma function, for the ratio of double sines one obtains
S_2(z)/S_2(z+g) =  e^∓πĝ(z-g^∗/2)(1+ O(d^-1(z,D))), where $z \in \mathbb{C} \setminus D$ and the sign $-$ (or $+$) is taken for $\Re z > 0$ (or $\Re z < 0$), see <cit.>.
Then from (<ref>) for the functions $\mu(x)$ and $K(x)$ (<ref>) with $x \in \R$ we have

\begin{equation}\label{Kmu-asymp}
\mu(x) \sim e^{\pi \hat{g} | x | \pm \imath \frac{\pi \hat{g} g^*}{2} }, \qquad K(x) \sim e^{- \pi \hat{g} |x|}, \qquad x\rightarrow \pm \infty.
\end{equation}

Denote also
ν_g= ĝ.Under the conditions (<ref>), (<ref>) by using the same Stirling formula we also have bounds

|μ(x)| ≤C e^πν_g |x|,     |K(x)| ≤C e^-πν_g |x|,      x ∈where $C$ is a positive constant uniform for compact subsets of parameters $\bo, g$ preserving the mentioned conditions, see <cit.>.

Another key result that we need in the paper is the following Fourier transform formula given in <cit.>, which we rewrite in terms of the double sine function using connection formula (<ref>).  This Fourier transform can be already found in [4, 14].
[18]
For real positive periods $\o_1, \o_2$ we have

\begin{equation}
\begin{aligned}
&\int_{\mathbb{R}} dx \, e^{\frac{2\pi\imath}{\o_1 \o_2} y x} S_2\Bigl(\imath x - \imath \nu + \frac{\o_1 + \o_2}{2} \Bigr) S_2^{-1} \Bigl( \imath x - \imath \rho + \frac{\o_1 + \o_2}{2} \Bigr) \\[6pt]
&= \sqrt{\o_1 \o_2} \, e^{ \frac{\pi\imath}{\o_1 \o_2} y (\nu + \rho) } S_2(\imath \rho - \imath \nu) \, S_2^{-1}\Bigl(\imath y + \frac{\imath(\rho - \nu)}{2} \Bigr) \, S_2^{-1}\Bigl(-\imath y + \frac{\imath(\rho - \nu)}{2} \Bigr),
\end{aligned}
\end{equation}

while the parameters $\nu, \rho, y$ satisfy the conditions

\begin{equation}\label{A19}
-\frac{\o_1 + \o_2}{2} < \Im \rho < \Im \nu < \frac{\o_1 + \o_2}{2}, \qquad |\Im y | < \Im \frac{\nu - \rho}{2}.
\end{equation}

In the special case

\begin{equation}
\nu = \frac{\imath g}{2}, \qquad \rho = -\frac{\imath g}{2}
\end{equation}

taking $y = \o_1 \o_2 \l$ and using homogeneity of the double sine (<ref>) (with $\gamma = \o_1 \o_2$) we arrive at the Fourier transform formula for the function $K(x)$ (<ref>)

\begin{equation}\label{K-fourier}
\int_{\mathbb{R}} dx \; e^{2 \pi \imath  \lambda x }  \K (x) = \sqrt{\omega_1 \omega_2} \, S_2(g) \, \KK (\lambda),
\end{equation}

where $| \Im \l | < \Re \hat{g}/2$ and conditions (<ref>) are satisfied due to the inequalities on the coupling constant $g$ (<ref>), (<ref>). Here we recall the notations (<ref>), (<ref>)

\begin{equation}
\hat{K}(\l) = K_{\hat{g}^*}(\l|\hat{\o}), \qquad \hat{g}^* = \frac{g^*}{\o_1\o_2}, \qquad \hat{\bo} = \Bigl( \frac{1}{\o_2}, \frac{1}{\o_1} \Bigr).
\end{equation}

Note that the right hand side of (<ref>) is analytic function of $\o_1, \o_2$ in the domain $\Re \o_j > 0$. The integral from the left is also analytic with respect to periods. Indeed, due to the bound (<ref>) it is absolutely convergent uniformly on compact sets of parameters $\bo, g$ preserving the conditions (<ref>), (<ref>). Hence, the formula (<ref>) also holds for complex periods under the mentioned conditions.

§ INEQUALITIES WITH ABSOLUTE VALUE

First, we will prove a little lemma that will be crucial for all other estimates.

For any $\epsilon \in [0, 2]$, $y_1, y_2, y \in \mathbb{R}$ we have

\begin{equation}
|y_1 - y_2| - |y_1 - y| - |y_2 - y| \leq \epsilon \left( |y_1| + |y_2| - |y| \right).
\end{equation}

Consider two cases. First, assume that $|y_1| + |y_2| \geq |y|$. Then using triangle inequality we obtain

\begin{equation}
|y_1 - y_2| - |y_1 - y| - |y_2 - y| \leq 0 \leq \epsilon \left( |y_1| + |y_2| - |y| \right)
\end{equation}

since $\epsilon \geq 0$.

Second, assume $|y_1| + |y_2| \leq |y|$. Now using $|y_j - y| \geq |y| - |y_j|$ and $|y_1 - y_2| \leq |y_1| + |y_2|$ we arrive at

\begin{equation}
|y_1 - y_2| - |y_1 - y| - |y_2 - y| \leq 2 \left( |y_1| + |y_2| - |y| \right) \leq \epsilon \left( |y_1| + |y_2| - |y| \right)
\end{equation}

because $\epsilon \leq 2$.

Note that Lemma <ref> holds for any norm and points $y_j$ in the corresponding normed space.

Denote by $\bm{y}_k$ a vector with $k$ components $\bm{y}_k = \bigl( y^{(k)}_{1}, \dots, y^{(k)}_{k} \bigr)$, $y_j^{(k)}\in \R$. Define a function $S_n$ by a recurrence relation

\begin{equation}\label{Srec2}
\begin{aligned}
S_n ( \bm{y}_1, \dots, \bm{y}_n ) = \sum_{\substack{i, j = 1 \\ i \not= j}}^n \bigl| y^{(n)}_i - y^{(n)}_j \bigr| &- \sum_{i = 1}^{n} \sum_{j = 1}^{n - 1} \, \bigl| y^{(n)}_i - y^{(n - 1)}_j \bigr| \\
&+S_{n - 1} ( \bm{y}_1, \dots, \bm{y}_{n - 1} )
\end{aligned}
\end{equation}

with $S_1 = 0$. Also by $\| \bm{y}_k \|$ denote $L^1$-norm

\begin{equation}
\| \bm{y}_k \| = \sum_{j = 1}^k \bigl| y^{(k)}_j \bigr|.
\end{equation}

The inequality

\begin{equation}\label{Sest}
S_n \leq \frac{1}{2}\sum_{\substack{i, j = 1 \\ i \not= j}}^n \bigl| y^{(n)}_i - y^{(n)}_j \bigr| + c_n \epsilon \| \bm{y}_n \|  - \epsilon \sum_{k = 1}^{n - 1} \| \bm{y}_k \|,
\end{equation}

holds for any $\epsilon \in \left[0, \frac{2(n - 1)}{c_n} \right]$, where the number $c_n$ is defined by recurrence relation $c_n = (n - 1) (1 + c_{n - 1})$ with $c_1 = 0$.

The proof goes by induction. Consider the case $n = 2$. Using Lemma <ref> we have

\begin{equation}
\begin{aligned}
S_2 = 2 \bigl| y_1^{(2)} - y_2^{(2)} \bigr| - \bigl| y_1^{(2)} - & y_1^{(1)} \bigr| - \bigl| y_2^{(2)} - y_1^{(1)} \bigr| \\[6pt]
&\leq \bigl| y_1^{(2)} - y_2^{(2)} \bigr| + \epsilon \left( \bigl| y_1^{(2)} \bigr| + \bigl| y_2^{(2)} \bigr| \right) - \epsilon \bigl| y_1^{(1)} \bigr|
\end{aligned}
\end{equation}

for any $\epsilon \in [0, 2]$. Thus, we proved the base case.

Now suppose we proved the statement for $S_{n - 1}$. Both sides of the stated inequality are symmetric with respect to components of $\bm{y}_k$ for all $k = 1, \dots, n$. Therefore, without loss of generality we can assume the ordering

\begin{equation}\label{ord}
y^{(k)}_1 \geq \ldots \geq y^{(k)}_k, \qquad k = 1, \dots, n.
\end{equation}

For the vector $\bm{y}_n$ with ordered components we can write

\begin{equation}\label{ord-comp}
\sum_{\substack{i, j = 1 \\ i \not= j}}^n \bigl| y^{(n)}_i - y^{(n)}_j \bigr| = 2 \sum_{m = 1}^{\lfloor n/2 \rfloor} (n - 2m + 1) \bigl| y_m^{(n)} - y_{n - m + 1}^{(n)} \bigr|.
\end{equation}

Consequently, due to the recurrence relation (<ref>) and induction assumption

\begin{equation}\label{Sineq}
\begin{aligned}
S_n &\leq 2 \sum_{m = 1}^{\lfloor n/2 \rfloor} (n - 2m + 1) \bigl| y_m^{(n)} - y_{n - m + 1}^{(n)} \bigr| - \sum_{i = 1}^n \sum_{j = 1}^{n - 1} \, \bigl| y^{(n)}_i - y^{(n - 1)}_j \bigr| \\[5pt]
&+ \sum_{m = 1}^{\lfloor (n - 1)/2 \rfloor} (n - 2m) \bigl| y_m^{(n - 1)} - y_{n - m}^{(n - 1)} \bigr|  + c_{n - 1} \epsilon \| \bm{y}_{n - 1} \| - \epsilon \sum_{k = 1}^{n - 2} \| \bm{y}_k \|.
\end{aligned}
\end{equation}

Next step it to regroup and estimate terms from the first three sums. For the term with $m = 1$ from the first sum we use the Lemma <ref> and for all other terms we simply use triangle inequalities.

Consider term with $m = 1$ from the first sum and terms with $i = 1, n$ from the second double sum and write the estimate

\begin{equation}\label{ineq1}
\begin{aligned}
(n - 1) \bigl| y_1^{(n)} - y_{n}^{(n)} \bigr| - \sum_{j = 1}^{n - 1} & \left( \bigl| y_1^{(n)} - y^{(n - 1)}_j \bigr| + \bigl| y_{n}^{(n)} - y^{(n - 1)}_j \bigr| \right) \\[4pt]
& \leq (n - 1) \, \epsilon' \left( \bigl| y_1^{(n)} \bigr| + \bigl| y_n^{(n)} \bigr| \right) - \epsilon' \, \| \bm{y}_{n - 1} \|,
\end{aligned}
\end{equation}

where we used Lemma <ref> with parameter $\epsilon' \in [0, 2]$ multiple times. Similarly let us estimate the term with $m > 1$ from the first sum together with the corresponding terms from the second double sum

\begin{equation}\label{ineq2}
\begin{aligned}
(n - 2m + 1) &\bigl| y_m^{(n)} - y_{n - m + 1}^{(n)} \bigr| \\[4pt]
&- \sum_{j = m}^{n - m} \left( \bigl| y_m^{(n)} - y^{(n - 1)}_j \bigr| + \bigl| y_{n - m + 1}^{(n)} - y^{(n - 1)}_j \bigr| \right) \leq 0,
\end{aligned}
\end{equation}

where we used triangle inequality multiple times. Remaining from the second double sum terms can be grouped with terms from the third sum

\begin{equation}\label{ineq3}
\begin{aligned}
(n - 2m) & \bigl| y_m^{(n - 1)} - y_{n - m}^{(n - 1)} \bigr| \\[4pt]
&- \sum_{i = m + 1}^{n - m} \left( \bigl| y_i^{(n)} - y^{(n - 1)}_m \bigr| + \bigl| y_i^{(n)} - y^{(n - 1)}_{n - m} \bigr| \right) \leq 0,
\end{aligned}
\end{equation}

where we again used triangle inequalities.

So, starting from (<ref>) and using inequalities (<ref>), (<ref>), (<ref>) we obtain

\begin{equation}
\begin{aligned}
S_n \leq \sum_{m = 1}^{\lfloor n/2 \rfloor} (n - 2m + 1) \bigl| y_m^{(n)} - y_{n - m + 1}^{(n)} \bigr| + (n - 1) \epsilon' \left( \bigl| y_1^{(n)} \bigr| + \bigl| y_n^{(n)} \bigr| \right) \\
+ (c_{n - 1} \epsilon - \epsilon') \| \bm{y}_{n - 1} \| - \epsilon \sum_{k = 1}^{n - 2} \| \bm{y}_k \|.
\end{aligned}
\end{equation}

Finally, for the first sum we again use equality (<ref>). For the second sum we have

\begin{equation}
\bigl| y_1^{(n)} \bigr| +  \bigl| y_{n}^{(n)} \bigr| \leq \| \bm{y}_{n} \|.
\end{equation}

Choosing $\epsilon' = (1 + c_{n - 1} )\epsilon$ we get bound stated in the lemma. Note that since $\epsilon' \leq 2$ then $\epsilon \leq 2/(1 + c_{n - 1}) = 2(n - 1)/c_n$.

The inequality

\begin{equation}\label{Sest2}
S_n \leq (n - 1 +  \epsilon) \| \bm{y}_n \| - \frac{\epsilon}{c_n} \sum_{k = 1}^{n - 1} \| \bm{y}_k \|
\end{equation}

holds for any $\epsilon \in \left[ 0, 2(n - 1) \right]$.

Use Lemma <ref> rescaling the parameter $\ve \rightarrow \ve/c_n$ and

\begin{equation}
\frac{1}{2} \sum_{\substack{i, j = 1 \\ i \not= j}}^n \bigl| y^{(n)}_i - y^{(n)}_j \bigr| \leq (n - 1) \| \bm{y}_n \|
\end{equation}

to arrive at stated inequality.

In fact, Lemma <ref> gives an estimate for the integrand in the integral representation of the wave function $\Psi_{\bl_n}(\bx_n)$ (<ref>) and, as a consequence, the exponential bound for the wave function. Denote, as before (<ref>),
μ'(_n) = ∏_i,j=1
i<j^n μ(x_i - x_j).
The wave function admits the bound
| μ'(_n) Ψ__n(_n) | ≤C expπν_g (  _n  -2/ν_g _n ł_n   )
with some $C(g, \bo)$ for any $\ve \in(0, 2(n - 1)]$, assuming $x_j \in \R$ and

\begin{equation}
| \Im(\l_k - \l_j) |< \frac{\nu_g}{2(n - 1)!e} \, \ve, \qquad k,j=1, \dots,n.
\end{equation}

The function $\mu'(\bx_n)$ is estimated as

\begin{equation}\label{B43}
|\mu'(\bx_n)| \leq C_1 \exp \pi \nu_g \Biggl( \; \sum_{\substack{i,j=1 \\ i<j}}^n |x_i-x_j| \, \Biggr)
\end{equation}

with some $C_1(g, \bo)$, where we used bound from (<ref>).
The integral representation (<ref>) in full form

\begin{equation}
\begin{aligned}
\Psi_{\bm{\lambda}_n} (\bm{x}_n) = C_\Psi \int_{\mathbb{R}^{n(n-1)/2} } & d\bm{y}_{n - 1} \dots d\bm{y}_1 \; e^{ {2\pi \imath}{} \left[ \lambda_n \bbx_n + (\lambda_{n - 1} - \lambda_{n}) \bby_{n - 1} \right] } \K (\bm{x}_n, \bm{y}_{n - 1}) \\[5pt]
& \hspace{0.6cm} \times  \prod_{k = 2}^{n - 1} e^{{2\pi \imath}{}(\lambda_{k - 1} - \lambda_{k}) \bby_{k - 1}} \, \mu   ( \bm{y}_{k} ) \, \K  (\bm{y}_{k}, \bm{y}_{k - 1} ),
\end{aligned}
\end{equation}

where $C_\Psi$ contains all constants $d_k$. Denote the integrand by $H$. Assuming

\begin{equation}
|\Im(\l_{k - 1} - \l_{k})| \leq \delta_\Lambda \frac{\nu_g}{2}
\end{equation}

and using bounds (<ref>) we arrive at

\begin{equation}
\begin{aligned}
|H| \leq C_2 \exp \pi\nu_g \Biggl(  -\frac{2}{\nu_g} \bbx_n \Im\l_n  &+\delta_\Lambda \sum_{k = 1}^{n - 1}\|\by_k\| \\
&+ S_n(\by_1, \dots,\by_{n - 1}, \bx_n) - \sum_{\substack{i,j=1\\i\not= j}}^n|x_i - x_j| \Biggr)
\end{aligned}
\end{equation}

with some $C_2(g,\bo)$. Using Lemma <ref> with rescaling $\ve \rightarrow \ve/c_n$ we have

\begin{equation}\label{B47}
|H| \leq C_2 \exp \pi\nu_g \Biggl(  -\frac{2}{\nu_g} \bbx_n \Im\l_n  + \Bigl(\delta_\Lambda - \frac{\ve}{c_n} \Bigr) \sum_{k = 1}^{n - 1}\|\by_k\| + \ve \|\bx_n\| - \frac{1}{2} \sum_{\substack{i,j=1\\i\not= j}}^n|x_i - x_j| \Biggr)
\end{equation}

for any $\ve \in [0, 2(n - 1)]$. For

\begin{equation}
\delta_\Lambda < \frac{\ve}{(n - 1)! e} < \frac{\ve}{c_n},
\end{equation}

see (<ref>), the bound (<ref>) represents integrable function. Combining it with (<ref>) we arrive at the bound stated in the corollary.

For the next inequality define another function that depends on vectors $\bm{y}_k$ and additional vector $\bm{t}_n = (t_1, \dots, t_n)$

\begin{equation}\label{Tdef}
T_n(\bm{y}_1, \dots, \bm{y}_n, \bm{t}_n) = \sum_{\substack{i, j = 1 \\ i \not= j}}^n | t_i - t_j | - \sum_{i, j = 1}^n \bigl| t_i - y_j^{(n)} \bigr| + S_n(\bm{y}_1, \dots, \bm{y}_n).
\end{equation}

The inequality

\begin{equation}\label{Test}
\begin{aligned}
T_n \leq ( n + r ) \| \bm{t}_n \| - \frac{1 - r}{2 n \, c_n} \sum_{k = 1}^{n} \| \bm{y}_k \| - r \, \biggl| \sum_{j = 1}^n \bigl( t_j - y_j^{(n)} \bigr) \biggr|
\end{aligned}
\end{equation}

holds for any $r \in [0, 1]$.

Both sides of the stated inequality are symmetric with respect to components of the vectors $\bm{y}_k, \bm{t}_n$ (separately). Therefore, without loss of generality we assume ordering

\begin{equation}\label{ord2}
y^{(k)}_1 \geq \ldots \geq y_k^{(k)}, \qquad t_1 \geq \ldots \geq t_n, \qquad k = 1, \dots, n.
\end{equation}

Then, as in the proof of the previous lemma, we can write

\begin{equation}\label{tord}
\sum_{\substack{i, j = 1 \\ i \not= j}}^n | t_i - t_j | = 2 \sum_{m = 1}^{\lfloor n / 2 \rfloor} (n - 2m + 1) | t_m - t_{n - m + 1} |.
\end{equation}

Next, take arbitrary $r \in [0, 1]$ and use triangle inequalities for the second double sum in (<ref>) to obtain

\begin{equation}\label{tr}
- \sum_{i, j = 1}^n \bigl| t_i - y_j^{(n)} \bigr| \leq - \left(1 - \frac{r}{n} \right) \sum_{i, j = 1}^n \bigl| t_i - y_j^{(n)} \bigr| - r \, \biggl| \sum_{j = 1}^n \bigl( t_j - y_j^{(n)} \bigr) \biggr|.
\end{equation}

Combining two previous formulas (<ref>), (<ref>) and Lemma <ref> with the parameter

\begin{equation}
\epsilon = \frac{1 - r}{2n \, c_n} \, \in \, \left[ 0, \frac{2(n - 1)}{c_n} \right]
\end{equation}

we arrive at inequality

\begin{equation}\label{Tineq1}
\begin{aligned}
T_n &\leq \sum_{m = 1}^{\lfloor n / 2 \rfloor} (n - 2m + 1) \left( 2 \, | t_m - t_{n - m + 1} | + \bigl| y_m^{(n)} - y^{(n)}_{n - m + 1} \bigr| \right) \\[6pt]
&- \left(1 - \frac{r}{n} \right) \sum_{i, j = 1}^n \bigl| t_i - y_j^{(n)} \bigr| + \frac{1 - r}{2n} \| \bm{y}_n \| - \frac{1 - r}{2n \, c_n} \sum_{k = 1}^{n - 1} \| \bm{y}_k \| - r \, \biggl| \sum_{j = 1}^n \bigl( t_j - y_j^{(n)} \bigr) \biggr|.
\end{aligned}
\end{equation}

Notice that two last sums from the right coincide with the corresponding terms in the stated inequality (<ref>), so it is left to bound the rest ones. Denote them by $R_n$

\begin{equation}\label{Rn}
\begin{aligned}
R_n = \sum_{m = 1}^{\lfloor n / 2 \rfloor} (n - 2m & + 1) \left( 2 \, | t_m - t_{n - m + 1} | + \bigl| y_m^{(n)} - y^{(n)}_{n - m + 1} \bigr| \right) \\[6pt]
&- \left(1 - \frac{r}{n} \right) \sum_{i, j = 1}^n \bigl| t_i - y_j^{(n)} \bigr| + \frac{1 - r}{2n} \| \bm{y}_n \|.
\end{aligned}
\end{equation}

We wish to prove the inequality

\begin{equation}\label{Rineq}
R_n \leq (n + r) \| \bm{t}_n \| - \frac{1 - r}{2n \, c_n} \| \bm{y}_n \|.
\end{equation}

Now fix $m \in \{ 1, \dots, \lfloor n / 2 \rfloor \}$ and from all of the remaining terms $R_n$ consider the following ones

\begin{equation}\label{m}
\begin{aligned}
R_{n, m} &= 2(n - 2m + 1) |t_m - t_{n - m + 1}| + (n - 2m + 1) \bigl| y_m^{(n)} - y^{(n)}_{n - m + 1} \bigr| \\[5pt]
& -\left(1 - \frac{r}{n} \right) \sum_{i = m}^{n - m + 1} \left( \bigl| t_i - y_m^{(n)} \bigr| + \bigl| t_i - y_{n - m + 1}^{(n)} \bigr| \right) \\[5pt]
&- \left(1 - \frac{r}{n} \right) \sum_{j = m + 1}^{n - m} \left( \bigl| t_m - y_j^{(n)} \bigr| + \bigl| t_{n - m + 1} - y_j^{(n)} \bigr| \right) + \frac{1 - r}{2n} \left( \bigl| y_m^{(n)} \bigr| + \bigl| y_{n - m + 1}^{(n)} \bigr| \right).
\end{aligned}
\end{equation}

Then for even $n$ we have

\begin{equation}\label{Revenn}
R_n = \sum_{m = 1}^{\lfloor n / 2 \rfloor} R_{n, m}
\end{equation}

and for odd $n$

\begin{equation}\label{Roddn}
R_n = \sum_{m = 1}^{\lfloor n / 2 \rfloor} R_{n, m} - \left( 1 - \frac{r}{n} \right) \bigl| t_{\frac{n + 1}{2}} - y_{\frac{n + 1}{2}}^{(n)} \bigr| + \frac{1 - r}{2 n} \, \bigl| y_{\frac{n + 1}{2}}^{(n)} \bigr|,
\end{equation}

since in this case two variables $t_{(n + 1)/2}$, $y_{(n + 1)/2}^{(n)}$ don't enter the first sum in (<ref>).

Let us prove the following bound

\begin{equation}\label{Rnm}
R_{n, m} \leq ( n + r ) \left( |t_m| + |t_{n - m + 1}| \right) - \frac{1 - r}{2n \, c_n} \left( \bigl| y_m^{(n)} \bigr| + \bigl| y_{n - m + 1}^{(n)} \bigr| \right).
\end{equation}

Note that for even $n$ the inequality we wish to prove (<ref>) directly follows from it due to (<ref>). For odd $n$ we have two more terms in (<ref>), but they can be easily bounded using triangle inequality

\begin{equation}
\begin{aligned}
- \left( 1 - \frac{r}{n} \right) \bigl| t_{\frac{n + 1}{2}} - y_{\frac{n + 1}{2}}^{(n)} \bigr| + \frac{1 - r}{2 n}  \, \bigl| y_{\frac{n + 1}{2}}^{(n)} \bigr|  &\leq \left( 1 - \frac{r}{n} \right) \bigl| t_{\frac{n + 1}{2}} \bigr| + \frac{1 + r - 2n}{2n} \, \bigl| y_{\frac{n + 1}{2}}^{(n)} \bigr| \\[10pt]
& \leq (n + r ) \bigl| t_{\frac{n + 1}{2}} \bigr| - \frac{1 - r}{2 n \, c_n}  \, \bigl| y_{\frac{n + 1}{2}}^{(n)} \bigr|.
\end{aligned}
\end{equation}

The proof of (<ref>) requires several steps. First, using triangle inequalities we can obtain

\begin{equation}\label{Rnm1}
\begin{aligned}
(n - 2m + 1) & \bigl| y_m^{(n)} - y^{(n)}_{n - m + 1} \bigr| \\[5pt]
&- \frac{n - 2m + 1}{n - 2m + 2} \sum_{i = m}^{n - m + 1} \left( \bigl| t_i - y_m^{(n)} \bigr| + \bigl| t_i - y_{n - m + 1}^{(n)} \bigr| \right) \leq 0.
\end{aligned}
\end{equation}

For the rest of the first sum in (<ref>) we use obvious inequality

\begin{equation}\label{Rnm2}
\begin{aligned}
-\biggl(  \frac{1}{n - 2m + 2} - \frac{r}{n} \biggr)  & \sum_{i = m}^{n - m + 1} \left( \bigl| t_i - y_m^{(n)} \bigr| + \bigl| t_i - y_{n - m + 1}^{(n)} \bigr| \right)\\[7pt]
&\hspace{-0.5cm} \leq -\frac{1 - r}{n} \, \Bigl( \bigl| t_m - y_m^{(n)} \bigr| + \bigl| t_m - y_{n - m + 1}^{(n)} \bigr|  \\[5pt]
&\hspace{2cm} + \bigl| t_{n - m + 1} - y_m^{(n)} \bigr| + \bigl| t_{n - m + 1} - y_{n - m + 1}^{(n)} \bigr| \Bigr).
\end{aligned}
\end{equation}

Then, we take the right-hand side of this formula and estimate it together with a part of $|t_m - t_{n - m + 1}|$ term using Lemma <ref> with the parameter $\epsilon$ = 1

\begin{equation}\label{Rnm3}
\begin{aligned}
\frac{2(1 - r)}{n} \, |t_m & - t_{n - m + 1}| - \frac{1 - r}{n} \, \Bigl( \bigl| t_m - y_m^{(n)} \bigr|  \\[7pt]
& + \bigl| t_m - y_{n - m + 1}^{(n)} \bigr|  + \bigl| t_{n - m + 1} - y_m^{(n)} \bigr| + \bigl| t_{n - m + 1} - y_{n - m + 1}^{(n)} \bigr| \Bigr) \\[7pt]
&\hspace{1.5cm} \leq \frac{1 - r}{n} \left( 2 |t_m| + 2 |t_{n - m + 1}| - \bigl| y_{m}^{(n)} \bigr| - \bigl| y_{n - m + 1}^{(n)} \bigr| \right).
\end{aligned}
\end{equation}

To bound the second sum in (<ref>) we again use triangle inequalities

\begin{equation}\label{Rnm4}
\begin{aligned}
(n - 2m) \biggl(1 & - \frac{r}{n} \biggr) |t_m - t_{n - m + 1}| \\[5pt]
&- \left(1 - \frac{r}{n} \right) \sum_{j = m + 1}^{n - m} \left( \bigl| t_m - y_j^{(n)} \bigr| + \bigl| t_{n - m + 1} - y_j^{(n)} \bigr| \right) \leq 0.
\end{aligned}
\end{equation}

Collecting results of the four inequalities (<ref>), (<ref>), (<ref>), (<ref>) we obtain the bound

\begin{equation}
\begin{aligned}
R_{n, m} \leq \biggl( n + r & + 2 - 2m - \frac{2 m r}{n} - \frac{2(1 - r)}{n} \biggr) | t_m - t_{n - m + 1} | \\[7pt]
&+ \frac{2 (1 - r)}{n} \left( |t_m| + |t_{n - m + 1}| \right) - \frac{1 - r}{2 n} \left( \bigl| y_{m}^{(n)} \bigr| + \bigl| y_{n - m + 1}^{(n)} \bigr| \right).
\end{aligned}
\end{equation}

Finally, to arrive at the inequality (<ref>) we use $|t_m - t_{n - m + 1}| \leq |t_m| + |t_{n - m + 1}|$ together with

\begin{equation}
2 - 2m - \frac{2 m r}{n} \leq 0, \hspace{1.5cm} \frac{1 - r}{2 n} \geq \frac{1 - r}{2 n \, c_n} \, .
\end{equation}

Thus, we proved the key bound (<ref>) and, consequently, the lemma.

[1] E. W. Barnes, The theory of the double gamma function, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 196 (1901), 265–387.

[2] N. Belousov, S. Derkachov, S. Kharchev, S. Khoroshkin, Baxter operators in Ruijsenaars hyperbolic system I. Commutativity of $Q$-operators, arXiv:2303.06383 (2023).

[3] L. D. Faddeev, Discrete Heisenberg-Weyl Group and modular group, Letters in Mathematical Physics 34 (1995), 249–254.

[4] L. D. Faddeev, R. M. Kashaev, A. Yu. Volkov, Strongly Coupled Quantum Discrete Liouville Theory. I: Algebraic Approach and Duality, Communications in Mathematical Physics 219:1 (2001), 199–219.

[5] A. Gerasimov, D. Lebedev, S. Oblezin, Baxter operator and Archimedean Hecke algebra, Communications in mathematical physics 284:3 (2008), 867–896.

[6] M. Hallnäs, S. Ruijsenaars, Joint Eigenfunctions for the Relativistic Calogero–Moser Hamiltonians of Hyperbolic Type: I. First Steps, International Mathematics Research Notices 2014:16 (2014), 4400–4456.

[7] M. Hallnäs, S. Ruijsenaars, Joint Eigenfunctions for the Relativistic Calogero–Moser Hamiltonians of Hyperbolic Type II. The Two-and Three-Variable Cases, International Mathematics Research Notices 2018:14 (2018), 4404–4449.

[8] M. Hallnäs, S. Ruijsenaars, Joint eigenfunctions for the relativistic Calogero–Moser Hamiltonians of hyperbolic type. III. Factorized asymptotics, International Mathematics Research Notices 2021:6 (2021), 4679–4708.

[9] S. Kharchev, S. Khoroshkin, Wave function for  $GL(n,\R)$ hyperbolic Sutherland model, arXiv:2108.04895 (2021).

[10] S. Kharchev, S. Khoroshkin, Wave function for $ GL (n,\mathbb {R}) $ hyperbolic Sutherland model II. Dual Hamiltonians, arXiv:2108.05393 (2021).

[11] Y. Kajihara, M. Noumi, Multiple elliptic hypergeometric series. An approach from the Cauchy determinant, Indagationes Mathematicae 14:3-4 (2003), 395–421.

[12] N. Kurokawa, S-Y. Koyama, Multiple sine functions, Forum Mathematicum 15 (2003), 839–876.

[13] I. Macdonald, Symmetric function and Hall Polynomials, Second edition, Oxford, Oxford University Press (1995).

[14] B. Ponsot, J. Teschner, Clebsch–Gordan and Racah–Wigner Coefficients for a Continuous Series of Representations of $ U_q (sl(2, \R))$, Communications in Mathematical Physics 224:3 (2001), 613–655.

[15] S. N. M. Ruijsenaars, Complete integrability of relativistic Calogero-Moser systems and elliptic function identities, Communications in Mathematical Physics 110 (1987), 191–213.

[16] S. N. M. Ruijsenaars, First-order analytic difference equations and integrable quantum systems, Journal of Mathematical Physics 38 (1997), 1069–1146.

[17] S. N. M. Ruijsenaars, Zero-eigenvalue eigenfunctions for differences of elliptic relativistic Calogero-Moser Hamiltonians, Theoretical and mathematical physics 146:1 (2006), 25–33.

[18] S. N. M. Ruijsenaars, A relativistic conical function and its Whittaker limits, SIGMA. Symmetry, Integrability and Geometry: Methods and Applications 7 (2011), 101.


G. Fejos, A. Patkos,
A Renormalized large-n solution of the U(n) x U(n) linear sigma model in the broken symmetry phase
Phys. Rev. D 82, (2010) 045011.

[45]
A. Jakovac and Zs. Szep,
Strange mass dependence of the tricritical point in the $U(3)_L \times U(3)_R$ chiral sigma model.
Phys. Rev. D 82, (2010) 125038.

[46]
L. Ferroni, V. Koch, and M. B. Pinto,
Multiple Critical Points in Effective Quark Models.
Phys. Rev. C82, (2010) 055205.

[47]
G. Marko and Zs. Szep, Influence of the Polyakov loop on the chiral phase transition in the two flavor chiral quark model,
Phys. Rev. D 82, (2010) 065021.

[48]A. Mocsy, I. N. Mishustin, and P. J. Ellis,
Role of fluctuations in the linear sigma model with quarks.
Phys. Rev. C70, (2004) 015204.

[49]
B.-J. Schaefer and J. Wambach, The Phase diagram of the quark meson model,
Nucl. Phys. A 757, (2005) 479.

[50]
B.-J. Schaefer and J. Wambach,
Susceptibilities near the QCD (tri)critical point
Phys. Rev. D 75, (2007) 085015.

[51]
B. J. Schaefer, J. M. Pawlowski, and J. Wambach, The Phase Structure of the Polyakov–Quark-Meson Model,
Phys. Rev. D  76, (2007) 074023.

[52]
B. J. Schaefer, M. Wagner, and J. Wambach, Thermodynamics of (2+1)-flavor QCD: Confronting Models with Lattice Studies,
Phys. Rev. D 81, (2010) 074013 .

[53]
H. Mao, J. Jin, and M. Huang, Phase diagram and thermodynamics of the Polyakov linear sigma model with three quark flavors,
J. Phys. G 37, (2010) 035001.

[54]
U. S. Gupta and V. K. Tiwari, Meson Masses and Mixing Angles in 2+1 Flavor Polyakov Quark Meson Sigma Model and Symmetry Restoration Effects,
Phys. Rev. D 81, (2010) 054019.

[55]
R. D. Pisarski and F. Wilczek, Remarks on the Chiral Phase Transition in Chromodynamics,
Phys. Rev. D 29, (1984) 338.

[56]
M. A. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot, On the phase diagram of QCD,
Phys. Rev. D 58, (1998) 096007.

[57]
V. Skokov, B. Friman, E. Nakano, K. Redlich, and B.-J. Schaefer,
Vacuum fluctuations and the thermodynamics of chiral models
Phys. Rev. D 82, (2010) 034029.

[58]
R. Khan and L. T. Kyllingstad,
The chiral phase transition and the role of vacuum fluctuations
AIP  Conf. Proc. 1343, (2011) 504 .

[59]
U. S. Gupta, V. K. Tiwari,
Revisiting the Phase Structure of the Polyakov-quark-meson Model in the presence of Vacuum Fermion Fluctuation
Phys. Rev. D 85, (2012) 014010 .

[60]
B.-J. Schaefer and M. Wagner, QCD critical region and higher moments for three flavor models,
Phys. Rev. D 85, (2012) 034027.

[61]
S. Chatterjee and K. A. Mohan, Including the Fermion Vacuum Fluctuations in the $(2+1)$ flavor Polyakov Quark Meson Model,
Phys. Rev. D 85, (2012) 074018.

[62]
J. O. Andersen and A. Tranberg, The Chiral transition in a magnetic background: Finite density effects and the functional renormalization group,
J. High Energy Phys. 08 (2012) 002.

[63]
V. K. Tiwari, Exploring criticality in the QCD-like two quark flavour models,
Phys. Rev. D 86, (2012) 094032.

[64]
S. Chatterjee and K. A. Mohan, Fluctuations and Correlations of Conserved Charges in the $(2+1)$ Polyakov Quark Meson Model,
Phys. Rev. D 86, (2012) 114021.

[65]
V. K. Tiwari, Comparing symmetry restoration trends for meson masses and mixing angles in the QCD-like three quark flavor models,
Phys. Rev. D  88, (2013) 074017.

[66]
T. K. Herbst, J. M. Pawlowski, and B.-J. Schaefer, Phase structure and thermodynamics of QCD,
Phys. Rev. D 88, (2013) 014007.

[67]
J. Weyrich, N. Strodthoff, and L. von Smekal, Chiral mirror-baryon-meson model and nuclear matter beyond mean-field approximation,
Phys. Rev. C 92, (2015) 015214.

[68]
P. Kovacs, Zs Szep, Gy Wolf, Existence of the critical endpoint in the vector meson extended linear sigma model,
Phys. Rev. D 93, (2016) 114014.

[69]
Andreas Zacchi and Jürgen Schaffner-Bielich, Effects of Renormalizing the chiral SU(2) Quark-Meson-Model,
Phys. Rev. D 97, (2018) 074011.

[70]
Andreas Zacchi and Jürgen Schaffner-Bielich, Implications of the fermion vacuum term in the extended SU(3) quark meson model on compact star properties,
Phys. Rev. D 100, (2019) 123024.

[71]
S. K. Rai and V. K. Tiwari, Exploring axial $U(1)$ restoration in a modified 2+1 flavor Polyakov quark meson model,
Eur. Phys. J. Plus 135, (2020) 844.

[72]
K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov
Generic rules for high temperature dimensional reduction and their application to the standard model
Nucl. Phys. B 458, (1996) 90.

[73]
P. Adhikari, J. O. Andersen and P. Kneschke, On-shell parameter fixing in the quark-meson model,
Phys. Rev. D 95, (2017) 036017.

[74]
S. Carignano, M. Buballa and B-J Schaefer, Inhomogeneous phases in the quark-meson model with vacuum fluctuations,
Phys. Rev. D 90, (2014) 014033.

[75]
J. O. Andersen, W. R. Naylor, and A. Tranberg, Phase diagram of QCD in a magnetic field: A review,
Rev. Mod. Phys. 88, (2016) 025001.

[76]
S. Carignano, M. Buballa, and W. Elkamhawy,
Consistent parameter fixing in the quark-meson model with vacuum fluctuations
Phys. Rev. D 94, (2016) 034023.

[77]
P. Adhikari, J. O. Andersen and P. Kneschke, Inhomogeneous chiral condensate in the quark-meson model,
Phys.Rev.D 96, (2017) 016013.

[78]
P. Adhikari, J. O. Andersen and P. Kneschke, Pion condensation and phase diagram in the Polyakov-loop quark-meson model,
Phys.Rev.D 98, (2018) 074016.

[79]
A.Folkestad, J. O. Andersen, Thermodynamics and phase diagrams of Polyakov-loop extended chiral models,
Phys.Rev.D 99, (2019) 054006.

[80]
S. K. Rai and V. K. Tiwari, On-shell versus curvature mass parameter fixing schemes in the quark-meson model and its phase diagrams,
Phys.Rev.D 105, (2022) 094010.

[81]
J. I. Kapusta and C. Gale, Finite Temperature Field Theory Principles and Applications (Cambridge University Press,
Cambridge, England), (2006).

[82]
R. Kobes, G. Kunstatter, and A. Rebhan, QCD plasma parameters and the gauge dependent gluon propagator,
Phys. Rev. Lett. 64, (1990) 2992.

[83]
R. Kobes, G. Kunstatter, and A. Rebhan, Gauge dependence identities and their application at finite temperature,
Nucl. Phys. B 355, (1991) 1.

[84]
A. K. Rebhan, The NonAbelian Debye mass at next-to-leading order, Phys. Rev.  D 48, (1993) R3967.