\section{Introduction}
There are various types of integrable
models of statistical mechanics
\cite{Baxter:1982zz} where interacting ``spins'' are assigned to
different parts of the square lattice. These includes the
\emph{Ising-type\/} (edge-interaction), \emph{vertex\/} and
\emph{iteraction-round-a-face\/} (IRF) models, where the Boltzmann
weights of local spin configurations are attributed to the edges,
vertices or faces of the lattice, respectively. The integrability of
these models requires the local Boltzmann weights to satisfy a
relevant type of the Yang-Baxter equation
\cite{McGuire:1964,Yang:1967,Baxter:1972}, which for the 
Ising-type models usually takes the form of Onsager's ``star-triangle
relation'' \cite{Onsager:1944}. 

The vertex type solutions of the Yang-Baxter equation (YBE), commonly called
$R$-matrices, are well understood. 
Their classification is based on the theory of Quantum Groups 
\cite{Drinfeld:1987,Jimbo:1986a,Faddeev:1987ih}. 
In particular, there are
infinite series of the so-called trigonometric solutions 
\cite{Cherednik:1980ey,Perk:1981,Zamolodchikov:1980ku,
IK81,Belavin:1982,Bazhanov:1984gu,Jimbo:1986,
Leites:1984pt,Bazhanov:1986av,
Kirillov:1987,
Khoroshkin:1991cmp,Warnaar:1992,
Delius:1994,Boos:2014,Chicherin:2014fqa},
connected with the evaluation representations of the
the quantized affine Lie algebras. 
Similar considerations also apply to the corresponding IRF
models, which are closely related to the vertex models through 
the IRF-vertex correspondence \cite{Bax73a,Bax73b,Jimbo:1988}. 
Here we consider the case of the $U_q(\widehat{sl}(2))$ algebra,  
where $q$ is a deformation (or ``anisotropy'') parameter. 
The standard 
trigonometric $R$-matrix ${\cal R}(\lambda\,|\,s_1,s_2)$ in this case 
depends on a
spectral variable $\lambda$ and parameters $s_1$, $s_2$, which define two
highest weight representations of the algebra $U_q({sl(2)})$. The
associated lattice model is 
usually referred to as ``a higher spin generalization of the 6-vertex
model''. 

Contrary to this holistic ``quantum group'' picture, 
the algebraic nature of the star-triangle
relation and associated Ising-type models 
is much less clear. 
Of course, these models are connected to the vertex and IRF type models. 
For instance, there is a ``star-square'' transformation (see
\eqref{tworel} below) between the Ising-type and IRF models, which has
been applied to some specific models 
\cite{Jungling:1975,Baxter:1986df,Boos:1997}. 
However, these examples in no way suggest 
how to construct an Ising-type formulation for an arbitrary integrable
IRF model,
in particular, for the IRF version of the higher spin generalization of
the 6-vertex model,
mentioned above.

Fortunately, some missing pieces  
in understanding of relevant algebraic structures are coming from
the 3-dimensional interpretation
\cite{Bazhanov:1992jqa,Bazhanov:2005as}  
of the Yang-Baxter equation based on Zamolodchikov's tetrahedron equation 
\cite{Zamolodchikov:1980rus,Zamolodchikov:1981kf,Baxter:1983qc} and
properties of its solutions \cite{Baxter:1986phd,Sergeev:1995rt,
Sergeev:1999,Sergeev:2001,BMS08a,Mangazeev:2013spa}.
In fact, the present paper is
motivated by a remarkable result by Mangazeev
\cite{Mangazeev:2014gwa}, 
who used the above 3D
interpretation to construct a very simple representation for the
$U_q(\widehat{sl}(2))$ trigonometric 
$R$-matrix ${\cal R}(\lambda\,|\,s_1,s_2)$ expressing it via a
terminating basic hypergeometric series ${}_4\varphi_3$. We found that
this result allows one to construct the desired Ising-type formulation of the
higher spin 6-vertex model with arbitrary values of $\lambda$, $s_1$ and $s_2$.
In general, this corresponds to infinite-dimensional representations,
so that the lattice spins take arbitrary non-negative integer
values. However, when the parameters $2s_1$ and $2s_2$ take
non-negative integer values all equations admit
a finite-dimensional reduction, with the spins taking finite number of values. 
In particular, when $s_1=s_2=\frac{1}{2}$ the model reduces
to the usual 6-vertex model with two spin states.

To illustrate the last point 
consider the standard vertex configurations of the usual 6-vertex
model, which are shown below in two ways: (i) by the edge arrangements 
(with bold and thin edges) and (ii) by the corresponding integer ``heights''
arrangements, $a\in{\mathbb Z}$, as in the unrestricted 
solid-on-solid (SOS) model,\footnote{%
The rules are: (a) two faces separated by the thin edge have the same
height, (b) the face height increases by $+1$ when crossing a bold
edge from left to right or from bottom to top. Note, that this trivial
vertex-SOS correspondence is different from that used in \cite{Bax73a,Bax73b}.}
\begin{equation}\label{figsix}
\input figsix.tex
\end{equation}
where $\omega_1,\omega_2,\ldots,\omega_6$ denote the corresponding
Boltzmann weights. For all other configurations the weights are
assumed to vanish identically. Below we assume that 
$\omega_5=\omega_6$, noting that this does not really reduce generality, since
for a periodic boundary conditions the weights $\omega_5$, $\omega_6$
always appear in pairs.  Now we replace the
lattice vertices
by the faces of the dual square lattice and go over from the SOS to the IRF
formulation of the model, as shown below,  
\begin{equation}\label{star6v}
\input fig0.tex
\input starsquare1.tex   
\end{equation}
It is not difficult to check that the resulting IRF weights can
be represented 
in two equivalent forms
\begin{equation}\label{tworel}
\begin{array}{rcl}
\displaystyle {\cal W}(a,b,c,d)&=&\displaystyle
\omega_1\, \frac{v(b-d)}{v(a-c)}\sum_{n=\max(b,c)}^a\,w_1(a-n)\,w_2(n-b)\,
w_3(n-c)\,w_4(n-d)\\[.8cm]
&=&\;\displaystyle
\omega_1\,
\frac{v(a-b)}{v(c-d)}\ \ \sum_{n=d}^{\min(b,c)}\;w_4(a-n)\,w_3(b-n)\, 
w_2(c-n)\,w_1(n-d)\,,
\end{array}
\end{equation}
as Boltzmann weights of four-edge stars of an Ising-type model with four
different types of oriented diagonal edges.  The corresponding edge weights 
$w_i(a-b)$, where $i=1,2,3,4$, depend on the difference of
spins $a$ and $b$ at the ends of the edge
(with the arrow pointing from $b$ to $a$).
These weights vanish for negative arguments 
\begin{equation}\label{wzero}
w_i(n)=0\,,\qquad n<0\,,\qquad i=1,2,3,4\,,
\end{equation}
therefore both summations over the central spins in \eqref{star6v} 
are restricted to finite intervals. Moreover, there are weights
$v(a-b)$ and $v(a-b)^{-1}$ associated with horizontal and vertical
edges in \eqref{star6v}. These weights only enter into the equivalence
transformation factors in \eqref{tworel}, which  
cancel out the partition function for periodic boundary conditions.

The only edge weights
that are required to reproduce \eqref{figsix} are given by
\begin{equation}
\begin{array}{l}
w_1(0)=w_2(0)=w_3(0)=w_4(0)=1\,,\\[.5cm] 
w_1(1)=\epsilon \,
(\omega_5^2-\omega_3^{}\omega_4^{})/(\omega_1\omega_5)\,,\quad
w_2(1)=\epsilon\,\omega_3/\omega_5\,\quad w_3(1)=\epsilon\,\omega_4/\omega_5\,,\quad 
w_4(1)=\omega_5/(\epsilon \,\omega_1)\,,\\[.5cm]
w_4(2)= \omega_5^2\,
(\omega_1^{}\omega_2^{}+\omega_3^{}\omega_4^{}-\omega_5^2)/(\epsilon^{2}\,\omega_1^2\, \omega_3^{}\,\omega_4^{})\,,\qquad \epsilon=v(1)/v(0)\;.
\end{array}
\end{equation}
Evidently, the weights ${\cal W}(a,b,c,d)$ only depend on spin
differences. For instance, if one sets $d=0$ then 
the central spin $n$ 
in \eqref{tworel} could only take three values $0,1,2$. Finally note,
that Eq.\eqref{tworel} relating the Ising-type and IRF models is
usually called the ``star-square transformation'', while the second
line in \eqref{tworel} relating the Boltzmann weights of two different
four-edge stars is called the ``star-star relation''.


Let us now return to the general higher spin 6-vertex
model, discussed above. Its transformation to an 
Ising-type model follows exactly the same steps. First one replaces the
vertex model by a SOS and then by an IRF model, as explained
above. 
Note that the edge indices in the vertex model (as well as the face
spin differences across the edge) can now take arbitrary non-negative
integer values. 
Next, the resulting IRF weights are represented as
Boltzmann weights of the first four-edge star in \eqref{tworel} 
with certain edge weights $w_i(a-b)$ and $v(a-b)$ 
defined for all non-negative 
spin differences $a-b\in {\mathbb Z}_{\ge0}$ 
but having the same vanishing conditions \eqref{wzero}. As a result
the range of summation over the central spin is finite and 
determined by the same
inequalities as in \eqref{tworel}.
Essentially, this representation is equivalent to the Mangazeev's formula
for the $R$-matrix ${\cal R}(\lambda\,|\,s_1,s_2)$ for the higher spin
6-vertex model. 


The existence of the Ising-type formulation
for this case
rises questions about the primacy of the quantum group structures in
the description of vertex models. It suggests, that 
the Yang-Baxter equations for 
${\cal R}(\lambda\,|\,s_1,s_2)$ may now be derived from
simpler relations imposed on the edge weights. Indeed, we show that
it is a corollary of the star-star relation, which has the same form
\eqref{tworel} as in the two-state
case. It is presented in Sect.\ref{star-section}, see Fig.~\ref{sspic} therein.
The proof of this relation
is based on the second Sears's transformation formula
\cite{Sears:1951} for the basic
hypergeometric series ${}_4\varphi_3$. 


The star-star relation as a primary integrablity condition has
appeared previously \cite{Bazhanov:1992jqa} in the context of 
the $sl(n)$-generalized
chiral Potts model \cite{BKMS,Date:1990bs}. 
For all other known integrable Ising-type models the star-star
relation could also be derived (see e.g., \cite{Baxter:1986df}), 
but it is a corollary of a much simpler star-triangle relation.
Therefore, even though we already have
a proof of the star-star
relation, mentioned above, it is
important to understand whether there are simpler integrability
conditions in the 6-vertex model and its generalizations, 
in particular, is there a star-triangle relation? 

Starting with Mangazeev's
${}_4\varphi_3$ formula for the $R$-matrix it was natural to first conduct the
search in the area of the basic hypergeometric series, especially because
it is extremely well documented, see, e.g., \cite{Andrews:1999,
Gasper:2004,Gasper:1995}. A review of this literature has eventually 
revealed that the required star-triangle relation can be derived from
1910 Jackson's $q$-analog \cite{Jackson:1910} 
of the Pfaff-Saalsch\"utz summation formula, 
originally
discovered by Johann Friedrich
Pfaff in 1797 \cite{Pfaff:1797}\footnote{%
The paper by Pfaff was a part of 1793 collection of the ``Analytical
Observations to L.Euler's Institute of Integral Calculus'' which was
published under the heading ``History 1793'' 
in ``New Transactions of Imperial Academy of Sciences in
St.Petersburg'', volume {\bf 11}. This volume is listed as published in 1793,
but Pfaff's paper is marked as ``Presented to Academy 14 January
1797''. This is probably the reason of a frequent 
volume-year mismatch in referencing
this paper, though after two hundred years the exact publication date
might not be that important.}
 and rediscovered by 
Saalsch\"utz in 1890 \cite{Saalschutz:1890}. An interesting feature of
this star-triangle relation is that it does not have an apparent
``difference property'', i.e., some of its edge weight depend on two
spectral variables rather than their ratio (we use multiplicative 
spectral variables). Nevertheless, we show that this relation implies 
the star-star relation and the Yang-Baxter equation and therefore  
does, indeed, plays the role of a primary integrability condition
in the higher spin generalization of the 6-vertex model. 

The organization of the paper is as follows. 
In Sect.~2 we formulate
a new inhomogeneous Ising-type model on the square lattice. 
In Sect.~3,  we present an equivalent IRF model and discuss its
integrability conditions, including the Yang-Baxter equation,
the star-star and the star-triangle relations. Connection to the
6-vertex model and its higher spin generalization is discussed in
Sect.~4, where also study most important finite-dimensional
reductions of the model. Factorized $R$-matrices are considered in
Sect.~5. Next we derive all functional relations for transfer 
matrices and Baxter's ${\bf Q}$-operators (Sect.6) and present the
Bethe ansatz solution (Sect 7.) for a general
inhomogeneous model. The column homogeneous case is considered in
Sect.~8. Some generalizations and extensions of the results to other
models are discussed in Sect.~9. In Conclusion we briefly summarize
a few important aspects of our work and mention some possible applications. 






\section{An Ising-type model\label{isingtype}}
Here we introduce a new two-dimensional solvable edge-interaction  
model. The model can be formulated on a rather general planar
graphs, however, for the purposes of this presentation it is convenient to
take a regular square lattice. 
Consider an oriented square lattice
drawn diagonally as in
Fig.~\ref{fig-lattice}. The edges of the lattice 
are shown by thick lines and the sites 
are shown with either open or filled circles in a
checkerboard order. We will refer to the latter as to the 
``white'' or ``black'' 
sites, respectively.
The edges of the lattice are oriented as indicated by arrows, 
namely, all the SE-NW edges
have the same (SE-NW) direction, while the SW-NE edges are oriented 
in a checkerboard order, always pointing towards the black sites.   

\bigskip
\begin{figure}[ht]
\input{fig1c.tex}
\caption{\small 
An oriented square lattice and its medial lattice shown by thick and
thin lines, respectively. All edges and lines are oriented as
indicated by arrows. There are two types of edges distinguished by relative 
orientation of the edge and the associated thin lines, passing through
the edge, shown in \eqref{twotypes}.
\label{fig-lattice}}
\end{figure} 

The thin vertical and horizontal
lines in Fig.~\ref{fig-lattice} represent the {\em medial}\/ graph 
whose vertices lie on the edges of the original square lattice. 
The lines are directed as shown by arrows. Each of these lines
carry its own parameter (in general complex) which we call 
the ``rapidity'' or  
``spectral variable''. In general, these variables may all be
different for different lines as illustrated in Fig.~\ref{fig-lattice}.
The minimal level of
generality corresponds to the homogeneous case 
when one assigns the same variable $x$ to
all lines directed upwards (i.e., when 
$x_1=x_2=\ldots=x_N\equivx$) 
and the same variable $x'$ to all lines directed
downwards (i.e., $x'_1=x'_2=\ldots=x'_N\equivx'$). 
Similarly, one assigns the same variables $y$ and $y'$ to all
thin horizontal lines directed right and left, respectively. 

At each lattice site place an integer spin variable $a\in {\mathbb Z}$. 
Two spins interact only if they are connected with an edge. 
There are two types
of edges distinguished by relative orientation of the edge and the
associated thin lines, passing through the edge,
\begin{equation}
\input{fig2c.tex}
\label{twotypes}
\end{equation}
The corresponding Boltzmann weights $V_{x/y}(a-b)$ and
$\big(V_{x/y}(a-b)\big)^{-1}$ depend on the difference of
spins $a$ and $b$ at the ends of the edge 
and on the ratio of two thin line variables $x$ and $y$. 
The edge arrow pointing from $b$ to $a$ 
shows that the argument of the weight function is
$(a-b)$ rather than $(b-a)$. 
For an easier perception the edges are also 
distinguished by their color: the blue
color 
for the edges of type $(i)$ and the red color for those of the type $(ii)$.

The function $V_x(n)$ is defined as 
\begin{equation}
V_x(n)=\left(\frac{q}{x}\right)^n\,\frac{(x^2;q^2)_n}{(q^2;q^2)_n}\,,\label{Vdef}
\qquad n\in {\mathbb Z}\,,
\end{equation}
where $n$ is an integer, $n\in {\mathbb Z}$,
\begin{equation}\label{qprod0}
(x;q^2)_n=\prod_{j=0}^{n-1}(1-x\, q^{2j})\,,\qquad
(x;q^2)_{-n}=\frac{1}{(x q^{-2n};q^2)_n}\,,
\end{equation}
and $q$ is an arbitrary parameter of the model. Note that 
\begin{equation}
V_x(n)\equiv0\,,\qquad \mbox{for} \ n<0\,,\label{vanish}
\end{equation}
and 
\begin{equation}\label{inv}
\sum_{n=-\infty}^\infty V_{a-n}(x^{-1}) V_{n-b}(x) \;=\;\delta_{a,b}\,,
\end{equation}
where, due to 
\eqref{vanish}, the summation is
effectively restricted to a finite interval $a\ge n\ge b$. Obviously, the edge 
weights are ``chiral'', since 
\begin{equation}
V_x(n)\not=V_x(-n)\,,\qquad n\not=0\,.
\end{equation}
Moreover, the vanishing property \eqref{vanish} is potentially 
dangerous for the type $(ii)$ edges with negative spin difference arguments
in \eqref{twotypes}, since their Boltzmann weights then diverge. 
However, this does not create any real problems, since 
such spin configurations will be excluded, see below.   

Now we need to describe the boundary conditions, which play rather
important role in the definition of the model. 
First, we assume that the values of white spins at the bottom row
of the lattice (see Fig.~\ref{fig-lattice}) are ordered as
\begin{subequations}
\begin{equation}\label{bc1}
a_1 \le a_2 \le a_3 \le \ldots 
\end{equation}
and so on, starting from right to left
and, similarly, the values of white spins 
at the right vertical boundary of the
lattice are ordered as 
\label{bc}
\begin{equation}\label{bc2}
a_1 \le b_1 \le c_1 \le \ldots
\end{equation}
\end{subequations}
starting from bottom to top. 

The Boltzmann weight of a particular spin configuration of the 
whole lattice is equal to the product of all the edge weights. 
The lattice spin configuration will be called {\it admissible\/}  if all the 
spin differences appearing in the edge weight functions \eqref{twotypes} 
are non-negative (i.e., for admissible 
configurations the values of spins cannot decrease  
along the edge directions). Remarkably, for the boundary conditions
\eqref{bc} the admissible configurations exhaust all
lattice spin configurations with non-zero Boltzmann weight.
Of course, the exclusion of non-admissible configurations
containing only 
type $(i)$ edges with negative spin arguments is obvious, since their
weights vanish identically due to \eqref{vanish}. However, this is not
obvious for configurations containing type $(ii)$ edges with $a<b$ in
\eqref{twotypes}, since their weights then become infinite. Our statement
follows from the limiting procedure based on the replacement of
\eqref{vanish} by 
\begin{equation}
V_x(n)=O(\varepsilon)\,,\qquad \varepsilon\to0\,,\qquad n<0\,.\label{vanish2} 
\end{equation}
Then a careful analysis shows, that as a consequence of \eqref{bc}, 
all divergent $1/\varepsilon$ pole
factors
coming from the type $(ii)$ edges are always compensated by zeros in
$\varepsilon$ coming the type $(i)$ edges. Actually, in all cases the
degree of zero is always higher than the degree of pole, so that the Boltzmann
weights of all non-admissible spin configurations with the 
boundary conditions \eqref{bc} do indeed vanish. Some additional 
explanations of this point are given after \eqref{starbound} 
in the next section.  
 

Next, we impose 
(quasi) periodic boundary conditions in the horizontal direction.
The lattice in Fig.~\ref{fig-lattice} has $N$ sites per row. Clearly, 
for consistency with spin difference dependence of the Boltzmann weights, 
one needs to formulate the periodic boundary condition in terms of 
differences of spins on the left and right vertical boundaries. Here we 
choose the simplest option
\begin{equation}\label{quasi}
b_{N+1}-a_{N+1}=b_1-a_1\,,\qquad 
c_{N+1}-b_{N+1}=c_1-b_1\,,\qquad \mbox{etc.}
\end{equation}
by equating the jumps of boundary spins between successive rows on
both sides of the lattice. If required, the periodic boundary
conditions in the vertical direction should also be defined in terms
differences of boundary spins between successive columns of lattice.
Finally, introduce the ``boundary fields'' by assigning additional
weights to the spin jumps on the horizontal and vertical boundaries 
of the lattice
\begin{equation}\label{ffactor}
w_{\rm fields}=\Big( \omega_v^{a_2-a_1} \, \omega_v^{a_3-a_2} \,
\cdots\, \omega_v^{a_{N+1}-a_N} \Big)\times
\Big(\omega_h^{b_1-a_1} \, \omega_h^{c_1-b_1} \, \cdots \,\Big)
\,, 
\end{equation}
where $\omega_h$ and $\omega_v$ are arbitrary parameters. 




The partition function is defined as a sum of the 
products of edge weights \eqref{twotypes} and the field factors 
\eqref{ffactor} over all admissible spin configurations,
\begin{equation} 
Z=\sum_{\mbox{\scriptsize (admissible spins)}}
w_{\rm fields}\ \prod_{\mbox{\scriptsize (edges)}}\ 
 \mbox{ (edge weights \eqref{twotypes})}\,. \label{part1}
\end{equation}
Obviously, the weights 
remain unchanged if all
lattice spins are simultaneously incremented by the same 
integer. Therefore, one 
particular spin, say the spin $a_1$ at the bottom right corner of the
lattice in Fig.~\ref{fig-lattice}, should be set to some fixed value (we
chose $a_1\equiv0$) to 
avoid a repeated counting of equivalent spin configurations.

To summarize, we have defined an Ising-type model on the square
lattice, see Fig.~\ref{fig-lattice}, 
with integer-valued spins at the lattice sites. The spins interact
with their nearest neighbours (two spins interact only if they are
connected by an edge).  
The edge Boltzmann weights 
depends on the difference of spins at the ends of the edge and the
ratio of two spectral variables assigned to the medial graph lines
passing through the edge. There are two types of edges, which are
assigned with different Boltzmann weights defined in \eqref{twotypes},
\eqref{Vdef}. The admissible spin configurations are controlled by the
boundary conditions \eqref{bc}, which play an essential role in the
definition of the model.
In the homogeneous lattice case the model has six
arbitrary, in general complex, parameters: the ``anisotropy'' parameter
$q$, the boundary field parameters $\omega_h, \omega_v$ and three independent
ratios of the four spectral variables $x,x',y,y'$. Note, that this 
counting does not include the trivial overall normalization of the
Boltzmann weights.    
For a fully inhomogeneous lattice there are additional spectral variables: 
there are 
two such variables for each column and two variables for each row (these are
the variables $x_1,x_1',x_2,x_2',\ldots$ and
$y_1,y_1',y_2,y_2',\ldots$, respectively, shown in Fig.~\ref{fig-lattice}).  
Remind, the weights only depend on the ratios of these variables.

\section{Equivalent models and integrability conditions\label{integrability}}

\subsection{Interaction-round-a-face formulation\label{irf-form}}
The edge interaction model introduced above could be
reformulated either as a {\em vertex\/} model
or as an {\em interaction-round-a-face\/} (IRF) model\/. In fact, this can be
done in a 
variety of ways. Let us first describe an equivalent IRF model on the square
lattice formed by the white sites in Fig.~\ref{fig-lattice}. 
The lattice in Fig.~\ref{fig-lattice} 
can be divided into 
four-edge ``stars'' with white corner sites, 
such as the star 
shown Fig.~\ref{fig-star}.
\begin{figure}[ht]
\begin{center}
\input fig3.tex
\caption{An elementary four-edge star centred around a black site.
\label{fig-star}}
\end{center}
\end{figure}
Applying the rules \eqref{twotypes} one can write the IRF-type 
Boltzmann weights for this star as
\begin{equation}\label{irf1}
{\mathcal S}(a,b,c,d\,|\,x,x',y,y')=
\sum_{n=\max(b,c)}^{a}\frac{V_{x/y'}(a-n)V_{y'/x'}(n-b)
V_{y/x}(n-c)}{V_{y/x'}(n-d)}\,.
\end{equation}
The weights depend on the four corner spins $a,b,c,d\in {\mathbb Z}$, 
which are supposed to obey the relations
\begin{equation}\label{starbound}
a\ge b\ge d\,,\qquad a\ge c\ge d\,.
\end{equation}
Outside this domain the weights \eqref{irf1} are assumed to vanish
identically. 
It is important to understand the origin of these
restrictions as well as the range of summation over the central spin
in \eqref{irf1}. Consider the four-edge star at the bottom right
corner of Fig.~\ref{fig-lattice}. The boundary conditions \eqref{bc}
require the ordering of the values of spins $a_2\ge a_1$ 
at the bottom and $b_1\ge a_1$ at the right side of that star. In the
notation of Fig.~\ref{fig-star} this corresponds to 
\begin{equation}\label{starbound1}
c\ge d\,,\qquad b\ge d \,.
\end{equation}
Let us now examine the sum in \eqref{irf1}. To regularize the summand 
for the case when its denominator vanishes we replace the strict
vanishing condition \eqref{vanish} by the regularized version
\eqref{vanish2}. The function ${V_{y/x'}(n-d)}$ 
in the denominator of \eqref{irf1} vanishes if its argument is
negative, i.e., when $n<d$. However, due to \eqref{starbound1} this
means that in this case $n<b$ as well as $n<c$, therefore, at least 
two factors $V_{y'/x'}(n-b)$ and $V_{y/x}(n-c)$
in the numerator of \eqref{irf1} must vanish, i.e., the degree of zero
in $\varepsilon$ is higher than the degree of the pole. Therefore the
configurations with a vanishing denominator from the type $(ii)$ edges 
do not contribute. Thus, the non-vanishing
contributions to the sum \eqref{irf1} could only arise when
$n\ge\max(b,c)$. Next, note that 
the first factor $V_{x/y'}(a-n)$ in the numerator 
is non-zero only when $a\ge n$. This explains the range of summation
in  \eqref{irf1} as well as the domain 
\eqref{starbound} where the star weights \eqref{irf1} do not
vanish. 

Let us now return to the bottom right star in
Fig.~\ref{fig-lattice}. The conditions \eqref{starbound} imply an
ordering of spins at the top and left sides of that star,
\begin{equation}
b_2\ge b_1\,,\qquad b_2\ge a_2\,.
\end{equation}     
Combining these with \eqref{bc} we can now apply the same reasonings
to the stars which are either directly above the or the first one to
the left of the bottom right corner star, that we have started. 
In this way the ordering of white spins in the horizontal and vertical
directions, which is initially only implied by \eqref{bc} to the
boundary, will propagate to the whole lattice. \footnote{%
It is not difficult to
see that the this result is essentially equivalent to the statement 
about the admissible spin configurations with non-zero Boltzmann weights
for the edge-interaction formulation presented in the previous
section.
}   
Thus, we have reformulated our model as an IRF model on the square
lattice formed by the white sites, with Boltzmann weights given by
\eqref{irf1}, \eqref{starbound}.
The partition function \eqref{part1} of the edge-interaction model from the
previous Section can now be equivalently rewritten 
\begin{equation}\label{part2}
Z=\sum_{\mbox{\scriptsize(white\ spins)}}
w_{\rm fields}
\prod_{\mbox{\scriptsize(black-center\ stars)}}
{\mathcal S}(a,b,c,d\,|\,x,x',y,y')\,. 
\end{equation}

\subsection{Commuting transfer matrices}

Below it will be convenient to use the IRF weights 
\begin{equation}\label{W-def}
{\mathcal W}(a,b,c,d\,|\,\,x,x',y,y')=\frac{V_{y/y'}(b-d)}{V_{y/y'}(a-c)}\,
\ {\mathcal S}(a,b,c,d\,|\,x,x',y,y')\,,
\end{equation}
which only differ from \eqref{irf1} by simple equivalence transformation
factors. 
Consider the transfer matrix of the edge-interaction
model between the two bottom rows of white spins in
Fig.~\ref{fig-lattice}, with the quasi-periodic
boundary conditions \eqref{quasi},
\begin{equation}\label{quasi1}
b_{N+1}-a_{N+1}=b_1-a_1\,.
\end{equation}
Its matrix elements can be written using
the IRF weights \eqref{W-def}, 
\begin{equation}\label{tmat}
\big({{\mathbb T}}(y,y'\,|\,x_N,x_N',\ldots,x_1,x_1')\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}=\omega_h^{b_1-a_1}\,\omega_v^{a_{N+1}-a_1}\,\prod_{i=1}^N
{\mathcal W}(b_{i+1},b_i,a_{i+1},a_i\,|\,x_i,x'_i,y,y')
\end{equation}
with the $a$'s and $b$'s labeling matrix rows and columns,
respectively. 
The parameters $y_1,y_1'$ from 
Fig.~\ref{fig-lattice} have been replaced by $y,y'$. 
Note, that the equivalence transformation factors,
introduced in \eqref{W-def},
cancel out in the definition of the
transfer matrix (i.e., the expression \eqref{tmat} remains
unchanged, if the weights ${\mathcal W}$ therein are replaced by the weights ${\mathcal S}$).
Remind, that the weights ${\mathcal W}(a,b,c,d)$ take non-zero values only if
the corner spins $a,b,c,d$ lie in the domain \eqref{starbound}. 
As a result the matrix elements \eqref{tmat} are non-zero 
\begin{equation}\label{nonzero}
\big({{\mathbb T}}(y,y'\,|\,x_N,x_N',\ldots,x_1,x_1')\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}\not=0\,, 
\qquad \mbox{iff}\qquad a_i\le b_i\,,\qquad 
i=1,2,\ldots,N+1\,,
\end{equation}
only if $a_i\le b_i$ for all $i=1,2,\ldots,N+1$. 

The transfer matrix acts in the infinite-dimensional 
vector space ${\cal H}^{{\rm(IRF)}}_N$ spanned by the vectors $\psi_{\bf a}$ 
labeled by ordered spin sequences 
\begin{equation}\label{secdef}
\psi_{\bf a}\in {\cal H}^{\rm(IRF)}_N\,,\qquad 
{\bf a}=(a_{N+1},a_N,\ldots,a_{1})\,\qquad \,a_i\in {\mathbb
  Z},\,\qquad 
a_{N+1}\ge a_N\ge \cdots\ge a_{1}\,.
\end{equation} 
As a consequence of \eqref{quasi1}, 
the transfer matrix preserves an integer quantity 
\begin{equation}\label{down}
M=\sum_{i=1}^N (a_{i+1}-a_{i})
=a_{N+1}-a_1=b_{N+1}-b_1
=\sum_{i=1}^N (b_{i+1}-b_i)\,,
\end{equation}
which is similar to the conserved 
number of ``down arrows'' in the six-vertex model. For this reason the
dependence on the vertical field parameter $\omega_v$ in \eqref{tmat} 
becomes trivial. Therefore, in what follows we set $\omega_v=1$.


The lattice model defined above is integrable. Namely, the transfer matrices
\eqref{tmat} with different values of $y$ and $y'$ (but the same
sets of $\{x\}$  and $\{x'\}$) form a two-parameter commuting family
\begin{equation}\label{tcomm}
[{{\mathbb T}}(y,y'\,|\,x_N,x_N',\ldots,x_1,x_1'),
{{\mathbb T}}(z,z'\,|\,x_N,x_N',\ldots,x_1,x_1')]=0\,,\qquad \forall\, y,y',z,z'\,.
\end{equation}
This follows in a standard way from the IRF-type Yang-Baxter
equation for the weights \eqref{W-def},
\begin{equation}\label{ybe}
\begin{array}{l}
\displaystyle \sum_{a\in{\mathbb Z}} 
{\mathcal W}(g,a,b,f\,|\,{\boldsymbol x},{\boldsymbol y})\ 
{\mathcal W}(c,e,g,a\,|\,{\boldsymbol x},{\boldsymbol z})\ 
{\mathcal W}(e,d,a,f\,|\,{\boldsymbol y},{\boldsymbol z})\;=\;\\
[3mm]
\displaystyle \qquad\qquad\qquad\qquad\qquad=\;
\sum_{h\in{\mathbb Z}} 
{\mathcal W}(c,h,g,b\,|\,{\boldsymbol y},{\boldsymbol z})\ 
{\mathcal W}(h , d ,b , f\,|\,{\boldsymbol x},{\boldsymbol z})\ 
{\mathcal W}(c , e ,h , d\,|\,{\boldsymbol x},{\boldsymbol y})
\end{array}
\end{equation}
where we have introduced compact notations for pairs of spectral
variables 
\begin{equation}
\qquad {\boldsymbol x}=(x,x')\,,\qquad {\boldsymbol y}=(y,y')\,,\qquad {\boldsymbol z}=(z,z')\,,\qquad \mbox{etc.}\,,\label{xpair}
\end{equation}
since they 
always appear in pairs.  Note that due to the corner spin constrains 
\eqref{starbound} for the weights \eqref{W-def} 
both summations in \eqref{ybe} are
restricted to finite intervals. 


For future references mention here some elementary symmetry
transformations of the YBE \eqref{ybe}. In particular, it is not 
affected by the re-scaling of the weights 
\begin{equation}\label{W-tilde}
{\mathcal W}(a,b,c,d\,|\,{\boldsymbol x},{\boldsymbol y})\to
\frac{g_{{\boldsymbol x}}(c,d)\,
  g_{{\boldsymbol y}}(a,c)}{g_{\boldsymbol x}(a,b)\,g_{{\boldsymbol y}}(b,d)}\,
{\mathcal W}(a,b,c,d\,|\,{\boldsymbol x},{\boldsymbol y})\,,
\end{equation}
with arbitrary non-vanishing factors $g_{\boldsymbol x}(a,b)$. 




\subsection{The star-star relation\label{star-section}}
The proof of Yang-Baxter equation \eqref{ybe} 
is based on the so-called ``star-star'' relation 
for the edge weights \eqref{twotypes}, which we present below.
The Boltzmann weights \eqref{W-def}, \eqref{irf1} are associated with
the black-centered four-edge stars of the type shown in Fig.~\ref{fig-star}.
Writing them in full, one gets
\begin{equation}\label{W-def1}
{\mathcal W}(a,b,c,d\,|\,\,x,x',y,y')=\frac{V_{y/y'}(b-d)}{V_{y/y'}(a-c)}\,
\sum_{n=\max(b,c)}^{a}\frac{V_{x/y'}(a-n)V_{y'/x'}(n-b)
V_{y/x}(n-c)}{V_{y/x'}(n-d)}\,.
\end{equation}
Similarly, 
one 
can define IRF-type weights associated with the white-centered stars
in Fig.~\ref{fig-lattice} having black corner spins,
\begin{equation}\label{W-def2}
{\overline {\mathcal W}}(a,b,c,d\,|\,x,x',y,y')
=\frac{V_{x/x'}(a-b)}{V_{x/x'}(c-d)}
\,\sum_{n=d}^{\min(b,c)} \ \frac{V_{x/y'}(n-d)
V_{y'/x'}(c-n) 
V_{y/x}(b-n)}
{V_{y/x'}(a-n)}\,.
\end{equation}
It turns out that the two definitions (with suitably chosen pre-factors 
in \eqref{W-def1} and \eqref{W-def2}) lead to the same result,
\begin{equation}\label{ssrel}
{\mathcal W}(a,b,c,d\,|\,x,x',y,y')=\overline{{\mathcal W}}(a,b,c,d\,|\,x,x',y,y')\,.
\end{equation}
This relation is a
statement of equality of the Boltzmann weights of the two ``stars'' shown
graphically in Fig.~\ref{sspic}, 
so it can naturally be called the star-star relation.
\begin{figure}[ht]
\begin{center}
\input st-st_pic.tex
\end{center}
\caption{A graphical representation of the star-star relation
  \eqref{ssrel}.\label{sspic}} 
\end{figure}
A verification of this identity 
is based on a repeated application of the second Sears's 
transformation formula for basic hypergeometric series (see \cite{Sears:1951}).
Details of these calculations together with a subsequent proof of
the Yang-Baxter equation \eqref{ybe} are presented in the Appendix~\ref{appA}.

\subsection{The star-triangle relation}
The ``star-star'' relation as a primary 
integrability condition has previously appeared in
\cite{Baxter:1986phd,Bazhanov:1992jqa} 
for 3-dimensional solvable models
\cite{Zamolodchikov:1981kf,Bazhanov:1992jqa}, which also cover  
the $sl(n)$-generalized
chiral Potts model in 2D \cite{BKMS,Date:1990bs}. 
For all other known integrable 2D edge interaction models the star-star
relation could also be derived (see e.g., \cite{Baxter:1986df}), 
but it is a corollary a much simpler star-triangle relation.
Therefore, even though we already have
a proof of the star-star
relation \eqref{ssrel}, mentioned above (see Appendix~\ref{appA}), it is
important to understand whether there are simpler integrability
conditions in our model, in particular, the star-triangle relation. 
Taking into account that the topic of the hypergeometric series have
been intensively studied by many mathematician for more than two
centures, it appears that 
the relevant star-triangle relation --- if it exists at all
--- should have already been discovered. Indeed, a purposeful 
study of the literature shows that the desired relation is the 1910
Jackson's $q$-analog \cite{Jackson:1910} 
of the Pfaff-Saalsch\"utz formula, originally
discovered by Pfaff in 1797 \cite{Pfaff:1797} and rediscovered by 
Saalsch\"utz in 1890 \cite{Saalschutz:1890}.  For a brief introduction
see the review article \cite{Gasper:1995}.  

To present the star-triangle relation we need to introduce yet another
edge weight function $W_{x,y}(a-b)$, which depends on the difference 
of edge spins $a$ and $b$ 
and on two spectral variables $x$ and $y$ (but not just on their ratio),
\begin{equation}\label{Wxy}
W_{x,y}(n) 
\;=\; \left(\frac{y}{x}\right)^n \;
\frac{(x^2;q^2)_n}{(y^2;q^2)_n}\;,\qquad n\in{\mathbb Z}\,.
\end{equation}
It possesses the following symmetries
\begin{equation}
W_{x,y}(n) \;=\; W_{q/y,q/x}(-n)\;=\;\frac{1}{W_{y,x}(n)}\;=\;\frac{1}{W_{q/x,q/y}(-n)}\;,\qquad n\in{\mathbb Z}\,.
\end{equation}
Note, in particular, that
\begin{equation}
W_{x,q}(n)=V_x(n)\,,\qquad W_{q,x}(n)=\frac{1}{V_x(n)}\,.
\end{equation}
The weights $V$ and $W$ satisfy the following two star-triangle
relations, 
\begin{subequations}
\label{str}
\begin{equation}
\sum_d \,V_{y/z}(d-a)\, W_{z,x}(b-d)\, V_{x/y}(c-d)  \;=\;
W_{z,y}(b-c) \,V_{x/z}(c-a)  \,W_{y,x}(b-a)\,,\label{str1}
\end{equation}
and
\begin{equation}
\sum_d V_{y/z}(a-d)\,W_{z,x}(d-b)\,V_{x/y}(d-c)   \;=\;
W_{z,y}(c-b)\, V_{x/z}(a-c)\, W_{y,x}(a-b)\,.\label{str2}
\end{equation}
\end{subequations}
where due to \eqref{vanish} the summations are restricted to finite
intervals.  
The two relations are
corollaries of each other under the substitution
\begin{equation}\label{str-tr}
b\leftrightarrow c\,,\qquad x\to q/z\,,\qquad y\to q/y\,,\qquad z\to
q/x\,.
\end{equation}
As mathematical identities, 
they are equivalent to the most general form of the 
of the Pfaff-Saalshch\"utz-Jackson summation
formula for the terminating balanced series
$\phantom{|}_3\varphi_2$, see eq.(3.5.1) in \cite{Gasper:1995}. 
The details of the correspondence are given in Appendix~\ref{appB}, 
where we also
present a graphical representation of the star-triangle relation
\eqref{str} and define yet another integrable edge interaction model 
involving the edge weights \eqref{Wxy}. 
Moreover, in the same Appendix we show that 
this relation 
implies the star-star relation \eqref{ssrel} and, consequently, the
Yang-Baxter equation \eqref{ybe} and commutativity of the transfer
matrices \eqref{tcomm}. Therefore  the star-triangle relation
\eqref{str} does, indeed, plays the role of a primary integrability condition
in the model.


\section{Connection to the six-vertex model\label{connection}}
The reader might have already become perplexed 
on how exactly the above constructions are related to the six-vertex model
mentioned in the title of the paper. These connections are discussed below. 
\subsection{Higher spin generalizations of the six-vertex
  model\label{irf-vertex}}
In the previous section we have reformulated our edge interaction
(Ising-type) model as an IRF model on the square lattice. Its face
weights  
depend on differences of corner spins, 
therefore one can easily convert this IRF model 
into a vertex model by replacing the face weights by an equivalent
vertex $R$-matrix, 
\begin{equation}
{\mathcal R}(x,x',y,y')^{j_1,j_2}_{i_1,i_2}=\delta_{i_1+i_2,j_1+j_2} \,
{\mathcal W}(i_1+i_2,j_2,i_1,0\,|x,x',y,y')\,.\label{R-def}
\end{equation}
where the weights ${\mathcal W}$ are defined by \eqref{W-def1} (or, equivalently, by 
\eqref{W-def} and \eqref{starbound}).
\begin{figure}[ht!]
\begin{center}
\input fig5c.tex
\caption{A typical face of the square lattice
  with corner spins $a$, $b$, $c$, $d$ and the associated vertex of
  the dual lattice with the edge spin variables $i_1$, $j_1$, $i_2$, $j_2$. 
The spectral parameters $(\lambda_1,s_1)$ and $(\lambda_2,s_2)$ are related 
to $(x,x')$ and $(y,y')$ as in \eqref{pardef}.
\label{fig5c}}
\end{center}
\end{figure}
The matrix indices 
\begin{equation}\label{ij-def}
i_1=c-d\,,\qquad i_2=a-c\,,\qquad j_1=a-b\,,\qquad j_2=b-d\,,
\end{equation}
are defined as differences of corner spins, arranged as in Fig.~\ref{fig5c}. 
Due to \eqref{starbound} these indices take non-negative integer values
\begin{equation}\label{j-pos}
i_1,i_2,j_1,j_2\ge0\,,\qquad  i_1,i_2,j_1,j_2\in {\mathbb Z}\,.
\end{equation}
The emergence of the conservation law 
\begin{equation}
i_1+i_2=j_1+j_2\,,\label{ar-cons}
\end{equation}
appearing in the delta-function in \eqref{R-def} is a trivial
consequence of the definitions 
\eqref{ij-def}. Using the two representations for the IRF weights
\eqref{W-def1} and \eqref{W-def2} one can write \eqref{R-def} in the form
\begin{eqnarray}\label{R-def1}
{\cal R}(x,x',y,y')^{j_1,j_2}_{i_1,i_2}&=&\delta_{i_1+i_2,j_1+j_2}\ 
\frac{V_{y/y'}(j_2)}{V_{y/y'}(i_2)}\,
\sum_{n=0}^{\min(i_2,j_1)}\frac{V_{x/y'}(n)V_{y'/x'}(j_1-n)
V_{y/x}(i_2-n)}{V_{y/x'}(i_1+i_2-n)}\,\nonumber\\
&&\label{R-def2}\\
&=&\delta_{i_1+i_2,j_1+j_2}\ \frac{V_{x/x'}(j_1)}{V_{x/x'}(i_1)}
\,\sum_{n=0}^{\min(i_1,j_2)} \ \frac{V_{x/y'}(n)
V_{y'/x'}(i_1-n) 
V_{y/x}(j_2-n)}
{V_{y/x'}(i_1+i_2-n)}\,.\nonumber
\end{eqnarray}
With these
notations the Yang-Baxter equation \eqref{ybe} takes the form
\begin{equation}
\begin{array}{r}
\displaystyle\sum_{j_1,j_2,j_3}
{\cal R}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}\ \ 
{\cal R}(x,x',z,z')_{j_1,i_3}^{k_1,j_3}\ \
{\cal R}(y,y',z,z')_{j_2,j_3}^{k_2,k_3}{\qquad \qquad}
\\[.6cm]
\displaystyle =\sum_{j_1,j_2,j_3}
{\cal R}(y,y',z,z')_{i_2,i_3}^{j_2,j_3}\ \
{\cal R}(x,x',z,z')_{i_1,j_3}^{j_1,k_3}\ \
{\cal R}(x,x',y,y')_{j_1,j_2}^{k_1,k_2}
\end{array}\label{Rybe}
\end{equation}
Note that due to \eqref{j-pos} and \eqref{ar-cons} all summations
here are restricted to finite intervals. It is worth
mentioning some simple symmetries of this equations. 
Obviously, it is invariant under the
simultaneous transposition in all three vector spaces. 
Below we will use this symmetry 
in a combination with the rescaling of the weights as in  
\eqref{W-tilde}. It is easy to
check that \eqref{Rybe} is invariant under the transformation
\begin{equation}\label{R-tilde}
\begin{array}{rcl}
\displaystyle{{\cal R}}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}\to
\displaystyle\widetilde{{\cal R}}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}&=&\displaystyle
\Big(\frac{\displaystyle y^{\prime j_1}x^{j_2}}{\displaystyle y^{i_1}x^{\prime
    i_2}}\Big)^{-1}
 \;\mathcal{R}^T(x,x',y,y')_{i_1,i_2}^{j_1,j_2}\,
\end{array}
\end{equation}
for all three $R$-matrices (with corresponding substitutions of their
arguments).  Here the superscript $T$
denotes the matrix transposition of $\cal R$ in both spaces
\begin{equation} \label{T-trans}
{\mathcal R}^{T}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}=
{\mathcal R}(x,x',y,y')^{i_1,i_2}_{j_1,j_2}\,.
\end{equation}
The scaling factor in \eqref{R-tilde} exactly coincides to that in
\eqref{W-tilde} with $g_{{\boldsymbol x}}(a, b) = x^b/x'^a$ (with an account of
\eqref{R-def}, \eqref{ij-def}). 

Note, that that the $R$-matrix \eqref{R-def2} 
only depends on three independent  
ratios of the four spectral parameters $x,x',y,y'$. Therefore, 
it is useful to
introduce a new set of parameters  
\begin{equation}\label{pardef}
x=\lambda_1\,q^{-s_1}\,,\quad x'=\lambda_1\,q^{s_1}\,,
\quad y=\lambda_2\,q^{-s_2}\,,\quad y'=\lambda_2\,q^{s_2}\,,
\quad z=\lambda_3\,q^{-s_3}\,,\quad z'=\lambda_3\,q^{s_3}\,.
\end{equation}
Then it is easy to see that \eqref{R-def2} only depend on the
variables $\lambda_1/\lambda_2$, $s_1$ and $s_2$. Therefore we will
write it as 
\begin{equation} \label{RR-rel}
{\cal R}(\lambda_1/\lambda_2\,|\,s_1,s_2)={\cal R}(x,x',y,y')\,,
\end{equation}
provided the parameters are related as in \eqref{pardef}. 
With these
notations the Yang-Baxter equation \eqref{ybe} takes the form
\begin{equation}
\begin{array}{r}
\displaystyle\sum_{j_1,j_2,j_3}
{\cal R}(\lambda_1/\lambda_2\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}\ \ 
{\cal R}(\lambda_1/\lambda_3\,|\,s_1,s_3)_{j_1,i_3}^{k_1,j_3}\ \
{\cal R}(\lambda_2/\lambda_3\,|\,s_2,s_3)_{j_2,j_3}^{k_2,k_3}\phantom{\qquad
\qquad }
\\[.6cm]
\displaystyle =\sum_{j_1,j_2,j_3}
{\cal R}(\lambda_2/\lambda_3\,|\,s_2,s_3)_{i_2,i_3}^{j_2,j_3}\ \
{\cal R}(\lambda_1/\lambda_3\,|\,s_1,s_3)_{i_1,j_3}^{j_1,k_3}\ \
{\cal R}(\lambda_1/\lambda_2\,|\,s_1,s_2)_{j_1,j_2}^{k_1,k_2}
\end{array}\label{Rybe2}
\end{equation}

Mention also two symmetry relations 
\eqref{R-def2},
\begin{eqnarray}
{\cal R}(\lambda\,|\,s_1,s_2)_{j_1,j_2}^{i_1,i_2}&=&
\frac{V_{q^{-2s_1}}(i_1) V_{q^{-2s_2}}(i_2)}
{V_{q^{-2s_1}}(j_1) V_{q^{-2s_2}}(j_2)}\ 
{\cal R}(\lambda\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}\,,
\label{sym1}\\[.4cm]
{\cal R}(\lambda\,|\,s_2,s_1)_{i_2,i_1}^{j_2,j_1}&=&
{\cal R}(\lambda\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}\,,\label{sym2}
\end{eqnarray}
and the normalization conditions 
\begin{equation}
\begin{array}{rcl}\displaystyle
\lim_{\lambda\to0}
\,\mathcal{R}(\lambda\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}&=&\delta_{i_1,j_1}
\delta_{i_2,j_2} q^{-2(s_1-j_1)(s_2-j_2)+2s_1s_2}\;,
\\[.4cm]\displaystyle
\lim_{\lambda\to\infty}
\mathcal{R}(\lambda\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}&=&\delta_{i_1,j_1}
\delta_{i_2,j_2} q^{+2(s_1-j_1)(s_2-j_2)-2s_1s_2}\;. 
\end{array}\label{R-norm}
\end{equation}
which simply follow from \eqref{R-def2} for fixed values of $s_1$ and $s_2$.

We state that the $R$-matrix \eqref{R-def2} coincides with Mangazeev's
$R$-matrix \cite{Mangazeev:2014gwa} for the higher spin generalizations of the 
6-vertex model. It is the trigonometric
$R$-matrix associated with two infinite-dimensional evaluation
representations of $U_q(\widehat{sl_2})$ with the highest weights
$2s_1$ and $2s_2$. The details of the correspondence are presented in
Appendix~\ref{appC}. 

\subsection{Finite-dimensional reductions\label{red}}
As usual we consider the $R$-matrix \eqref{R-def2} as a matrix acting in the
tensor product of two vector spaces, assigning the matrix indices
$i_1,j_1$ and $i_2,j_2$ to the first and second spaces, respectively.
Following the terminology of the QISM we will call them as the ``quantum''
and ``auxiliary'' spaces. Having in mind further applications 
we will discuss various reductions associated with the second
(auxiliary) space, however, in view of the symmetry \eqref{sym2},
the same considerations also apply to the first
space.  
When $2s_2$ is a non-negative integer the corresponding evaluation 
representation of $U_q(\widehat{sl_2})$ becomes reducible and splits
into a semi-direct 
sum of the $(2s_2+1)$-dimensional and an infinite-dimensional
representations. The $R$-matrix then takes the block-triangular form in the
second space\footnote{%
Note that in  
the considered case the formula \eqref{R-def2} requires a limiting
procedure where a non-integer 
$2s_2$ is continuously
approaching the required integer value. In this way all the matrix
elements \eqref{R-def2} remain finite and well defined.}, 
\begin{equation}\label{twos2}
2s_2\in{\mathbb
    Z}_{\ge0}\,:\qquad \qquad
{\mathcal R}(\lambda\,|\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}=0\,,\qquad 0\le i_2\le 2s_2,
\qquad j_2> 2s_2\,,\qquad \forall\, i_1,j_1\,.
\end{equation} 
To select its finite-dimensional block one just needs
to restrict the indices to a finite set
\begin{equation}\label{fin}
i_2,j_2\in {\cal I}_{s_2}=\{0,1,2,\ldots,2s_2\}\,.
\end{equation}
Similarly, for the infinite-dimensional block one needs
to select the complimentary set of non-negative integers, 
\begin{equation}\label{inf}
i_2,j_2\in \overline{\cal I}_{s_2}=\{2s_2+1,2s_2+2,\ldots,\infty\}\,.
\end{equation}
The Yang-Baxter equation \eqref{Rybe2} turns into two equations of the
same form, where 
the indices $i_2,j_2,k_2$ related to the second space are restricted to
ether the finite ${\cal I}_{s_2}$ or the infinite set 
$\overline{\cal I}_{s_2}$. Below we consider a few examples of this reduction. 
\subsubsection{Inversion relation, $s_2=0\,.$}
Setting $s_2=0$ in \eqref{pardef}, using the definition
\eqref{R-def2} and the inversion relation 
\eqref{inv} one obtains for the finite-dimensional part \eqref{fin},
\begin{equation}\label{spin0}
s_2=0:\qquad\qquad
{\mathcal R}(\lambda\,|\,s_1,0\,)_{i_1,0}^{j_1,0}=\delta_{i_1,j_1}\,,
\end{equation}
where $\lambda=\lambda_1/\lambda_2$.
For the infinite-dimensional part \eqref{inf} one get an important identity 
\begin{equation}\label{rel0}
s_2=0:\,\qquad\qquad
{\mathcal R}(\lambda\,|\,s_1,0\,)_{i_1,i_2}^{j_1,j_2}=\,\frac{[\lambda q^{-s_1}]\,
[q^{i_2}]}{[\lambda q^{+s_1}]\,[q^{j_2}]}\ 
{\mathcal R}(\lambda\,|\,s_1,-1\,)_{i_1,i_2-1}^{j_1,j_2-1}\,,\qquad i_2,j_2\ge1\,,
\end{equation}
relating the $R$-matrices with $s_2=0$ and $s_2=-1$. Here and below we
use the notation
\begin{equation}
[x]=x-x^{-1}\,.
\end{equation}

\subsubsection{The $L$-operator, $s_2=\frac{1}{2}\,.$\label{Lop-red}}
The next important case to
consider is when $s_2=\frac{1}{2}$ and $s_1$ is arbitrary. 
With the relations \eqref{pardef} this corresponds to $y'=q\,y$.
Introducing the notation
\begin{equation}\label{Lop-first}
{\mathcal R}(x,x',y,q\,y)_{i_1,i_2}^{j_1,j_2}
=\frac{1}{[x'/y]}{\cal L}(x,x',y)_{i_1,i_2}^{j_1,j_2}\,,
\end{equation}
and selecting the
2-dimensional subspace 
${\cal I}_{\frac{1}{2}}$ of the second space and 
one obtains from \eqref{R-def2},
\begin{equation}
\bm{\mathcal{L}}(x,x',y)_{i_1}^{j_1}=
\displaystyle\begin{pmatrix}
{\mathcal{L}}_{i_1,0}^{j_1,0}&
\mathcal{L}_{i_1,0}^{j_1,1}\\[.4cm]
\mathcal{L}_{i_1,1}^{j_1,0}&
\mathcal{L}_{i_1,1}^{j_1,1}
\end{pmatrix}=
\displaystyle\begin{pmatrix}
\ \delta_{i_1,j_1} [q^{-i_1}x'/y]&\delta_{i_1,j_1+1} {[q^{i_1}]}\\[.4cm]
\delta_{i_1+1,j_1}[q^{-i_1}x'/x]&
\ \ \delta_{i_1,j_1}[q^{i_1}x/y]\end{pmatrix}\,,
\label{Lmat}
\end{equation}
where $[x]=x-1/x$.  Notice, that here for the arguments of ${\cal
  L}(x,x',y)$ we use the same (original) variables as in ${\cal
  R}(x,x',y,qy)$.   
As usual we consider the above
expressions as matrix elements of the $L$-operator, which is regarded
as a two by two
matrix in the second space, whose elements are operators acting in the 
first space,
\begin{equation}
\boldsymbol{{\cal L}}(x,x',y)
=\begin{pmatrix}
\mu\, q^{\frac{H}{2}}-\mu^{-1} \,q^{-\frac{H}{2}}&
F\\[3mm]
E& \mu \,q^{-\frac{H}{2}}-\mu^{-1} \,q^{\frac{H}{2}}
\end{pmatrix}\,,\quad q^{2s_1}=x'/x\,,\quad \mu=(x x')^{\tfrac{1}{2}}/y\,,
\label{Lop}
\end{equation}
where the dependent variable $\mu$ and was used to bring the above
formula to the standard form.
Here $E,F,H$ stand for the generators of 
the quantum universal enveloping algebra $U_q(sl_2)$ 
\begin{equation}
[H,E]=2 E\,,\qquad [H,F]=-2F\,,\qquad
[E,F]={(q-q^{-1})}\,(q^H-q^{-H})\,.\label{uqsl2}
\end{equation}
The matrix elements \eqref{Lmat} correspond to the
infinite-dimensional highest weight 
representation $\pi^+_{s_1}$ of this algebra 
(the parameter $s_1$  is defined by $x'/x=q^{2 s_1}$),
\begin{equation}\label{piplus_s1}
\pi^+_{s_1}[H]\,|j\rangle=(2s_1-2j)\,|j\rangle\,,
\qquad \pi^+_{s_1}[E]\,|j\rangle=[q^{2s_1+1-j}]\,|j-1\rangle\,,\qquad 
\pi^+_{s_1}[F]\,|j\rangle=[q^{j+1}]\,|j+1\rangle\,,
\end{equation}
spanned on the basis vectors $|j\rangle\in {\mathbb C}^{(\infty)}$ 
with $j=0,1,2,\ldots,\infty$.

\bigskip
For the infinite-dimensional part of the $R$-matrix 
with the indices $i_2,j_2\ge2$ one gets the relation,
\begin{equation}\label{rel1}
s_2={\textstyle\frac{1}{2}}:\,\qquad
{\mathcal R}\big(\lambda\,|\,s_1,{\textstyle\frac{1}{2}}\,\big)_{i_1,i_2}^{j_1,j_2}
=\,\frac{[\lambda q^{-s_1-\frac{1}{2}}]\,
[\lambda q^{-s_1+\frac{1}{2}}]\,
[q^{i_2}]\,[q^{i_2-1}]}{[\lambda q^{+s_1-\frac{1}{2}}]\,[\lambda
    q^{+s_1+\frac{1}{2}}]\,
[q^{j_2}]\,[q^{j_2-1}]}\ 
{\mathcal R}\big(\lambda\,|\,s_1,-{\textstyle\frac{3}{2}}\,\big)_{i_1,i_2-2}^{j_1,j_2-2}\,,
\qquad i_2,j_2\ge2\,,
\end{equation}
which connect the $R$-matrices with $s_2=\frac{1}{2}$ and $s_2=-\frac{3}{2}$. 
Note, that here we use the notations for ${\cal
  R}(\lambda\,|\,s_1,s_2)={\mathcal R}(x,x',y,y')$,
as it is defined in \eqref{RR-rel} and \eqref{pardef}.




\subsubsection{Six-vertex $R$-matrix, $s_1=s_2=\frac{1}{2}\,.$ }
Consider the case when
$s_1=s_2=\frac{1}{2}$. Setting 
\begin{equation} \label{Rhh}
{\cal
  R}(\lambda\,|\,{\textstyle\frac{1}{2}},{\textstyle\frac{1}{2}})
=\frac{1}{[\lambda q]}\,
{\cal R}^{\rm (6v)}(\lambda)\,,
\end{equation}
and 
using \eqref{R-def2}, \eqref{pardef}, 
one gets the standard $R$-matrix of the 6-vertex model,
\begin{equation}
\begin{array}{ll}
{\mathcal R}^{\rm (6v)}(\lambda)_{00}^{00}={\mathcal R}^{\rm (6v)}(\lambda)_{11}^{11}=
\rho(\lambda)\,(\lambda q-\lambda^{-1} q^{-1})\,,\qquad\qquad 
&{\mathcal R}^{\rm (6v)}(\lambda)_{01}^{10}={\mathcal R}^{\rm (6v)}(\lambda)_{10}^{01}=
\rho(\lambda)\,(q-q^{-1})\,,\\[.6cm]
{\mathcal R}^{\rm (6v)}(\lambda)_{01}^{01}={\mathcal R}^{\rm (6v)}(\lambda)_{10}^{10}=
\rho(\lambda)\,(\lambda -\lambda^{-1})\,,
&
\end{array}\label{R6v}
\end{equation}
where $\lambda=\lambda_1/\lambda_2$ and $\rho(\lambda)\equiv1$.
Alternatively, one could get the same expression from \eqref{Lop}
by choosing the 2-dimensional representation of \eqref{uqsl2} with
$s_1=\frac{1}{2}$.  

\subsubsection{Infinite-dimensional $R$-matrices, $2s_2\in{\mathbb
    Z}_{\ge0}\,. $}
Here we present the generalization of the relations \eqref{rel0},
\eqref{rel1},
\begin{equation}
{\mathcal R}\big(\lambda\,|\,s_1,s_2\,\big)_{i_1,i_2}^{j_1,j_2}
=\,\frac{V_{y/x}(2s_2+1)\,V_{q^{-2i_2}}(2s_2+1)}
{V_{y/x'}(2s_2+1)\,V_{q^{-2j_2}}(2s_2+1)}\ 
{\mathcal R}\big(\lambda\,|\,s_1,-s_2-1\,\big)_{i_1,i_2-2s_2-1}^{j_1,j_2-2s_2-1}\,,
\qquad i_2,j_2\ge2s_2+1\,,\label{twos2inf}
\end{equation}
valid for $2s_2\in{\mathbb Z}_{\ge0}$ and arbitrary $s_1$. The
factor in the RHS is expressed via the weight function \eqref{Vdef}, 
the parameter relations \eqref{pardef} are assumed and
$\lambda=\lambda_1/\lambda_2$. 



\subsubsection{Symmetry transformation}
Below we will also use an alternative variant of the $L$-operator
originating from the $R$-matrix \eqref{R-tilde},
\begin{equation}
\widetilde{{\mathcal R}}(x,x',y,q\,y)_{i_1,i_2}^{j_1,j_2}
=\frac{1}{[x'/y]}\widetilde{{\cal
    L}}(x,x',y)_{i_1,i_2}^{j_1,j_2}\,,\qquad i_2,j_2=0,1\,.
\end{equation}
It is simply related to \eqref{Lop-first} (note the matrix transposition in both spaces)
\begin{equation}\label{Lop-second}
\displaystyle\widetilde{{\cal L}}(x,x',y)_{i_1,i_2}^{j_1,j_2}=\displaystyle
\Big(\frac{q^{j_1}y^{j_1} x^{j_2} }{y^{i_1}{x'}^{i_2} }\Big)^{-1}\,
 {\cal
  L}(x,x',y)_{j_1,j_2}^{i_1,i_2}\,.
\end{equation}
Writing $\widetilde{\bm{\mathcal{L}}}$ as a two by two matrix 
in the second space, but keeping the matrix indices of the first space, 
one gets
\begin{equation}
\widetilde{\bm{\mathcal{L}}}(x,x',y)_{i_1}^{j_1}=\begin{pmatrix}
\displaystyle
q^{-i_1}\,[\mu\,q^{s_1-j_1}]\,\delta_{i_1,j_1}
& \displaystyle  \mu^{-1}\,q^{s_1-j_1}\,[q^{2s_1-j_1}]
\,\delta_{i_1,j_1+1}\\[7mm]
\displaystyle\mu\,q^{s_1-j_1}\,[q^{j_1}]\,\delta_{i_1,i_1-1}
&
\displaystyle  q^{2s_1-i_1}\,[\mu\,q^{-s_1+i_1}]\,\delta_{i_1,j_1}
\end{pmatrix}\label{Lmat4}
\end{equation}
Here we use the same notations as in \eqref{pardef} and $\mu=\lambda_1
q^{\frac{1}{2}}/\lambda_2$. Similar to \eqref{Lop} we can rewrite 
\eqref{Lmat4} in the operator form
\begin{equation}
\widetilde{\bm{\mathcal{L}}}\;=\;
\left(\begin{array}{cc}
\displaystyle q^{-s_1}\,\big(\mu \,q^H -\mu^{-1}\big) & \displaystyle \mu^{-1}\,F\,q^{\frac{H}{2}} \\
[5mm]
\displaystyle \mu\,E\,q^{\frac{H}{2}} & \displaystyle q^{+s_1}\,\big(\mu - \mu^{-1}\, q^H\big)
\end{array}\right)
\end{equation}
where $E,F,H$ denote the generators
of the quantum universal enveloping algebra $U_q(sl(2))$, defined by
\eqref{uqsl2}. The matrix elements \eqref{Lmat4} corresponds to the
infinite-dimensional representation 
\begin{equation}
\pi^+_{s_1}[H]\,|j\rangle=(2s_1-2j)\,|j\rangle\,,
\qquad \pi^+_{s_1}[E]\,|j\rangle=[q^{j}]\,|j-1\rangle\,,\qquad 
\pi^+_{s_1}[F]\,|j\rangle=[q^{2s_1-j}]\,|j+1\rangle\,,
\end{equation}
spanned on the basis vectors $|j\rangle$ with $j=0,1,2,\ldots,\infty$.
This representation is equivalent to \eqref{piplus_s1}.

It is instructive to rewrite \eqref{Lmat4} in terms of generators of the 
Weyl algebra
\begin{equation}\label{weyl}
\boldsymbol{u}\,\boldsymbol{v}\;=\;q^2\,\boldsymbol{v}\,\boldsymbol{u}\;.
\end{equation}
Choosing the representation
\begin{equation}\label{uvj}
\boldsymbol{u}\,|\, j\, \rangle \;=\; |\, j\, \rangle\, q^{2s_1-2j}\;,\quad
\boldsymbol{v}\,|\, j\, \rangle \;=\; |\, j-1\, \rangle\;,
\end{equation}
one gets
\begin{equation}\label{L-form}
\widetilde{\bm{\mathcal{L}}}\;=\;\left(\begin{array}{cc}
\displaystyle \mu \, q^{-s_1} \, \boldsymbol{u} \,-\, \mu^{-1} \, q^{-s_1} & \displaystyle 
\mu^{-1} \, q^{s_1} \, \boldsymbol{v}^{-1} \, ( \boldsymbol{u}\,-\, q^{-2s_1} ) \\
\\
\displaystyle \mu \, q^{s_1} \, \boldsymbol{v} \, ( 1\,-\,q^{-2s_1} \, \boldsymbol{u} ) & 
\displaystyle \mu \, q^{s_1} \,-\, \mu^{-1} \, q^{s_1} \, \boldsymbol{u}
\end{array}\right)\;.
\end{equation}

 


\subsection{New solutions of the Yang-Baxter equation}
The existence of the star-triangle relation \eqref{str} 
allows one in the standard
way to construct new solutions of the Yang-Baxter equation.
For completeness we present them in the Appendix~\ref{sappb}, 
though they are not relevant to the main topic of this paper.
From the quantum group point of view the $R$-matrix \eqref{R-def1} 
intertwines two infinite-dimensional 
{\em evaluation highest weight\/} representations of
the $U_q(\widehat{sl}(2))$ algebra. The new $R$-matrices 
presented in Appendix~\ref{sappb} are related to infinite-dimensional
representations without highest and lowest weights.




\section{Six-vertex model as a descendant of the six-vertex 
model\label{descen}}

The Ising-type model presented above is inspired by a  
remarkable structure of the Mangazeev’s $R$-matrix 
\cite{Mangazeev:2014gwa} for the higher spin
generalizations of the 6-vertex model, defined in \eqref{R-def2} (see 
Appendix~C for connections with the notations of \cite{Mangazeev:2014gwa}). 
In this section we demonstrate how the observed Ising-type structure can
be discovered by elementary calculations, based on the idea of
factorized $R$-matrices and Baxter's ``propagation
through the vertex'' techniques \cite{Baxter:1982zz}. Previously these
methods were used in \cite{Bazhanov:1989nc}
for the chiral Potts model \cite{Au-Yang:1987syu,McCoy:1987pt,Baxter:1987eq}
which is related to cyclic representations of $U_q(\widehat{sl}(2))$ with $q$
being a root of unity. In this way the chiral Potts model (which is an
Ising-type model) has turned out 
to be  a ``descendant'' of the six-vertex model. 
Here we show that the same
methods are equally powerful for arbitrary values of $q$ in the context
of the six-vertex model itself.


\subsection{Factorized $R$-matrices}

In Sect.~\ref{irf-form} and \ref{irf-vertex},
we have reformulated our Ising-type model first as 
an IRF model and then as 
a vertex model. It is easy to see that to within boundary
effects the same model can also be reformulated as yet another
vertex model. Indeed, 
the square lattice in Fig.~\ref{fig-lattice} 
can be formed by periodic translations of
the ``box diagrams'' shown in Fig.~\ref{fig-box}. 
One can choose any
of these two diagrams, for definiteness we take the left diagram in
Fig.~\ref{fig-box}.    
\begin{figure}[ht]
\begin{center}
\input fig-box_new.tex
\end{center}
\caption{Graphical representation of the factorized $R$-matrces. 
Eq.\eqref{def3} corresponds to the left diagram.}
\label{fig-box}
\end{figure} 
The Boltzmann weight of this box can be conveniently associated with
an R-matrix
\begin{equation}\label{def3}
{\mathbb R}({\boldsymbol x},{\boldsymbol y})_{a\,,\,b}^{a',\,b'}\;=\;
\frac{V_{y'/x'}(b-a') V_{x/y'}(a'-b') V_{y/x}(b'-a)}{V_{y/x'}(b-a)}\;,
\end{equation}
where ${\boldsymbol x}=(x,x')$ and ${\boldsymbol y}=(y,y')$ stand for the pairs of spectral
variables in the vertical and 
horizontal directions. The spins $a,b,a',b'\in{\mathbb Z}$ are
supposed to satisfy the relations 
\begin{equation}\label{dom2}
b\ge a'\ge b'\ge a\,.
\end{equation}
Outside this domain the RHS of \eqref{def3} is 
assumed to vanish identically. 
The above inequalities ensure that all the spin differences appearing
in \eqref{def3} 
are non-negative, as it is required for admissible lattice spin 
configurations defined in Sect.~\ref{isingtype}.

Naturally, as the reader might have expected, the $R$-matrix
\eqref{def3} satisfy the vertex form Yang-Baxter equation 
\begin{equation}\label{ybe4}
\sum_{a',\,b',\,c'}{\mathbb R}({\boldsymbol x},{\boldsymbol y})_{a\,,\,b}^{a',\,b'} \,
{\mathbb R}({\boldsymbol x},{\boldsymbol z})_{a'\;,\,c}^{a'',\,c'}\,
{\mathbb R}({\boldsymbol y},{\boldsymbol z})_{b'\;,\,c'}^{b'',\,c''}
\;=\;
\sum_{a',\,b',\,c'}
{\mathbb R}({\boldsymbol y},{\boldsymbol z})_{b\,,\,c}^{b',\,c'}\,
{\mathbb R}({\boldsymbol x},{\boldsymbol z})_{a\,,\,c'}^{a',\,c''}\,
{\mathbb R}({\boldsymbol x},{\boldsymbol y})_{a'\;,\,b'}^{a'',\,b''}\;.
\end{equation}
A sketch of the proof, based on the star-star relation \eqref{ssrel}, 
is presented
in the Appendix~\ref{box-YBE}. Note also, that 
apparently related factorized $R$-matrices
in an operator form were previously 
constructed in \cite{Faddeev:1994,Derkachev:2006}.

At first sight the above equation \eqref{ybe4} looks identical to
\eqref{Rybe} (apart from the relabelling of spin indices). The
difference is in the range of allowed values of these indices.
For \eqref{Rybe} the indices
can take arbitrary non-negative integer values and the restrictions on
summations are determined by the conservation laws of the type 
\eqref{ar-cons} for
every vertex. Contrary to this the external indices in \eqref{ybe4} 
are taking arbitrary integer values, satisfying the inequalities 
\begin{equation}
a\le c''\le b'' \le a'' \le c\,,\qquad a\le b \le c\,,
\end{equation}
resulting from the the conditions
\eqref{dom2} for each $R$-matrix.
By the same reason the summation indices in  \eqref{ybe4}
are restricted by  
\begin{equation}
a\leq b'\leq c''\;,\quad b'' \leq c' \leq a''\;,\quad b'\leq a' 
\leq \min (c',b) 
\end{equation}
for the LHS and
\begin{equation}
c''\leq a'\leq b''\;,\quad a'' \leq b' \leq c\;,\quad \max (a',b) \leq
c' \leq b' 
\end{equation}
for the RHS. Note, that all the summations in \eqref{ybe4} go over
finite intervals.   

A remarkable feature of the $R$-matrix \eqref{def3} is that its 
matrix elements factorize into a product of factors depending 
only on two spins connected by an edge.
On the other hand the $R$-matrices \eqref{def3} and 
\eqref{R-def2} describe the same model and, it is natural to expect
that they should be related 
to each other by a linear transformation.\footnote{%
In case when $q$ is a root
of unity such relations was obtained in
\cite{Baxter:1986df,Boos:1997}.} The corresponding relation is given
by a discrete Fourier transform 
\begin{equation}\label{fourier}
\sum_{a',b'\in{\mathbb Z}}\mathbb{R}({\boldsymbol x},{\boldsymbol y})_{a\,,\,b}^{a',\,b'}\ 
\Big(\frac{x}{x'}q^{2j_1}\Big)^{\varepsilon a'} 
\Big(\frac{y}{y'}q^{2j_2}\Big)^{\varepsilon b'}\;=\;
\sum_{i_1,i_2\ge0}\Big(\frac{x}{x'}q^{2i_1}\Big)^{\varepsilon a}
 \Big(\frac{y}{y'}q^{2i_2}\Big)^{\varepsilon b}\, 
\Big(\frac{\displaystyle y^{\prime j_1}x^{j_2}}{\displaystyle y^{i_1}x^{\prime
    i_2}}\Big)^{\varepsilon}
 \;\mathcal{R}^T({\boldsymbol x},{\boldsymbol y})_{i_1,i_2}^{j_1,j_2}
\end{equation}
where ${\boldsymbol x}=(x,x')$, ${\boldsymbol y}=(y,y')$, $\varepsilon=\pm1$,  and the superscript $T$
denotes 
the matrix transposition in both spaces
\begin{equation}\label{Rtr}
{\cal R}^T({\boldsymbol x},{\boldsymbol y})_{i_1,i_2}^{j_1,j_2}=
{\cal R}({\boldsymbol x},{\boldsymbol y})_{j_1,j_2}^{i_1,i_2}\,.
\end{equation}
Note, that Eq.\eqref{fourier} contains 
two different (but equivalent) relations with $\varepsilon=+1$ and 
$\varepsilon=-1$. A proof of \eqref{fourier} is presented below. 
Separate considerations of the second diagram (on the right side of
Fig.~\ref{fig-box}) are really not required, since its Boltzmann weight is
expressed through \eqref{def3} as ${\mathbb R}({\boldsymbol x},{\boldsymbol y})_{-a',-b'}^{-a,\,\,-b}$ and all relevant formulae can be obtained by simple symmetry transformations. 



\subsection{The $RLL$-relations.}


As a preparation to the proof of \eqref{fourier} consider a few important  
cases of the Yang-Baxter equation \eqref{Rybe}, 
which are sufficient for an {\em ab initio} calculation of 
the $R$-matrix \eqref{R-def2}. 
First, quote the defining relation for the $L$-operator,  
\begin{equation}
\begin{array}{r}
\displaystyle\sum_{j_1,j_2,j_3}
{\cal L}(x,x',y)_{i_1,i_2}^{j_1,j_2}\ \ 
{\cal L}(x,x',z)_{j_1,i_3}^{k_1,j_3}\ \
{\cal R}^{\rm (6v)}(y/z)_{j_2,j_3}^{k_2,k_3}\phantom{\qquad
\qquad }
\\[.6cm]
\displaystyle =\sum_{j_1,j_2,j_3}
{\cal R}^{\rm (6v)}(y/z)_{i_2,i_3}^{j_2,j_3}\ \
{\cal L}(x,x',z)_{i_1,j_3}^{j_1,k_3}\ \
{\cal L}(x,x',y)_{j_1,j_2}^{k_1,k_2}
\end{array}\label{LLR1}
\end{equation}
where ${\cal R}^{\rm (6v)}$ is given by \eqref{R6v}. 
This equation is obtained from \eqref{Rybe} via the
$2$-dimensional reduction (see \eqref{Lop-first} in Sect.~\ref{Lop-red}) 
in the second and third spaces when $y'=q y$
and $z'=q z$. The indices $i_2,j_2,k_2,i_3,j_3,k_3=0,1$ take two values, while 
the indices $i_1,j_1,k_1\ge0$ take arbitrary non-negative integer values.
It is well known
that with the substitution \eqref{Lop} the above equation just reduces to 
the commutation relations \eqref{uqsl2} of the algebra $U_q(sl(2))$. 
Thus, any representation of this algebra gives a solution of
\eqref{LLR1}. 
Here we choose the representation \eqref{piplus_s1},
which leads to the expression \eqref{Lmat} for the matrix elements of 
${{\cal L}}(x,x',y)_{i_1,i_2}^{j_1,j_2}$. 

Next, consider the $2$-dimensional reduction of \eqref{Rybe} in the
third space, when $z'=q z$,
\begin{equation}
\begin{array}{r}
\displaystyle\sum_{j_1,j_2\in{\mathbb Z}_{\ge0}}\ 
\bm{{\cal L}}(y,y',z)_{i_2}^{j_2}\ 
\bm{{\cal L}}(x,x',z)_{i_1}^{j_1}\ 
{\cal R}(x,x',y,y')_{j_1,j_2}^{k_1,k_2}\phantom{\qquad
\qquad }
\\[.6cm]
\displaystyle =\sum_{j_1,j_2}
\ {\cal R}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}\ 
\bm{{\cal L}}(x,x',z)_{j_1}^{k_1}\ 
\bm{{\cal L}}(y,y',z)_{j_2}^{k_2}\,.
\end{array}\label{LLR2}
\end{equation}
Note, that the boldface
$
\bm{{\cal L}}(x,x',z)_{i}^{j}$ and  
$\bm{{\cal L}}(y,y',z)_{i}^{j}$ 
are defined by \eqref{Lmat}.
They are considered as
two-by-two matrices and the above equation
involves their matrix product.  
As usual we call Eq.\eqref{LLR2} as the ``$RLL$-relation''. 
Once the $L$-operators $
\bm{{\cal L}}(x,x',z)$ and  
$\bm{{\cal L}}(y,y',z)$ are fixed by \eqref{Lmat} 
this relation determines ${\cal R}(x,x',y,y')$ to within an overall
normalization (see, e.g., Appendix~C of ref.\cite{Mangazeev:2014gwa}).
Furthermore,  
by using \eqref{fourier} this relation can
be transformed to another $RLL$-relation, which uniquely determines the
factorized $R$-matrix \eqref{def3}. 

To proceed, first note the Yang-Baxter equations \eqref{LLR1} and \eqref{LLR2} are invariant under the substitution 
\begin{equation}\label{t-trans}
\bm{{\cal R}}\to 
\widetilde{\bm{{\cal R}}}\,,\qquad\qquad
\bm{{\cal L}}\to 
\widetilde{\bm{{\cal L}}}\,,\qquad\qquad
\bm{{\cal R}}^{(6v)}\to 
\widetilde{\bm{{\cal R}}}^{(6v)}\,,\qquad\qquad
\end{equation}
where the modified $R$-matrices are defined by \eqref{R-tilde},
\eqref{Lop-second} and 
\begin{equation}\label{R6v-tilde}
\widetilde{\mathcal R}^{(6v)}(\lambda)_{i_1,i_2}^{j_1,j_2}=
q^{i_2-j_1}\,\lambda^{i_2-j_2}\,
{\mathcal R}^{(6v)}(\lambda)_{i_1,i_2}^{j_1,j_2}\,.
\end{equation}
As before, we will consider $\widetilde{\bm{{\cal
      L}}}_i^j$, defined by \eqref{Lmat4},
as a two-by-two matrix in the second space,
with operator entries acting in the first space 
(the indices $i,j$ label their elements). 
Changing the basis in this space via the formula
\begin{equation}\label{L-trans}
\sum_{a'}\,\bm{{\mathbbm L}}(x,x',y)_{a}^{a'}\ \Big(\frac{x}{x'} \,q^{2
  j}\Big)^{-a'} =\sum_{i} \Big(\frac{x}{x'}\,q^{2i}\Big)^{-a}\ 
\widetilde{\bm{{\mathcal L}}}(x,x',y)_i^j
\end{equation}
and using \eqref{Lmat4}, 
one obtains the matrix elements of the transformed $L$-operator
\begin{equation}\label{Lop-new}
\bm{{\mathbbm L}}(x,x',y)_{a}^{a'}=
\begin{pmatrix}\displaystyle
\frac{x}{y}\, \delta_{a,a'-1}-\frac{y}{x'}\,\delta_{a,a'}&
\displaystyle q^{-2a}\Big(\,\frac{y}{x}\,\delta_{a,a'-1}-\frac{y}{x'}\,
\delta_{a,a'}\,\Big)\\[.5cm]
\displaystyle q^{2a}\Big(\,\frac{x'}{y}\,\delta_{a,a'}
- \frac{x}{y}\,\delta_{a,a'-1}\,\Big)& 
\displaystyle \frac{x'}{y}\,\delta_{a,a'}-\frac{y}{x}\,\delta_{a,a'-1} 
\end{pmatrix}\,.
\end{equation}
Evidently, it still has the form \eqref{L-form}, but the Weyl algebra \eqref{weyl} is now realized as
\begin{equation}\label{uva}
\boldsymbol{u}\,|\, a\, \rangle \;=\; |\, a-1\, \rangle \;,\quad
\boldsymbol{v}\,|\, a\, \rangle \;=\; |\, a\, \rangle\, q^{2a}\;.
\end{equation}

Next, taking the ``tilded'' version of \eqref{LLR2}, obtained by the symmetry transformation \eqref{t-trans}, and using \eqref{L-trans} and the variant of \eqref{fourier} with
$\varepsilon=-1$ one obtains
\begin{equation}\label{RLL-irf}
\begin{array}{r}\displaystyle
\sum_{a',b'} \mathbb{R}(x,x',y,y')_{a\,,\,b}^{a',\,b'} \ 
\bm{{\mathbbm L}}(x,x',z)_{a'}^{a''}\  \bm{{\mathbbm
    L}}(y,y',z)_{b'}^{b''}{z\qquad \qquad \qquad}\\[.5cm]
=\displaystyle\sum_{a',b'}
\bm{{\mathbbm{L}}}(y,y',z)_{b}^{b'}\  \bm{{\mathbbm{L}}}(x,x',z)_{a}^{a'}
\ 
\mathbb{R}(x,x',y,y')_{a'\;,\,b'}^{a'',\,b''}\,.
\end{array}
\end{equation}
It should be stressed that this relation is equivalent to
\eqref{LLR2}. It is obtained from the latter by the symmetry transformation 
\eqref{t-trans}, followed by linear
transformation, changing basis in the first and second vector spaces. 
At this stage we regard \eqref{fourier} 
as definition of   
$\mathbb{R}(x,x',y,y')$. 
However, remembering that the new $RLL$-relation \eqref{RLL-irf} (with
the $L$-operators specified by \eqref{Lop-new}) determines
$\mathbb{R}(x,x',y,y')$ up an overall normalization, one concludes that 
in order to prove \eqref{fourier} it suffices to show
that the factorized expression \eqref{def3} satisfies \eqref{RLL-irf}.





\subsection{Propagation through the vertex techniques}
The matrix elements of the $L$-operator \eqref{Lop-new} possess
important factorization properties, which we describe below. Let 
\begin{equation}
\bm{v}_0=\begin{pmatrix}1\\0\end{pmatrix}\,,
\qquad \bm{v}_1=\begin{pmatrix}0\\1\end{pmatrix}\,,\qquad
\overline{\bm{v}}_i=\bm{v}^t_i\,,\qquad i=0,1\,,\label{v-basis}
\end{equation}
be the two-dimensional basis vectors.
Introduce the following set of vectors labelled by $a,a'\in{\mathbb Z}$ 
\begin{equation}\label{sos-vertex}
\bm{\Phi}(x,y)_{a}^{a'}=\sum_{i=0,1}\Phi(x,y)_{a,i}^{a'}\, 
\bm{v}_i\,,\qquad
\overline{\bm{\Phi}}(x,y)_{a}^{a'}=
\sum_{i=0,1}\overline{\Phi}(x,y)_{a,i}^{a'}\, 
\overline{\bm{v}}_i\,,
\end{equation}
with the components 
\begin{equation}\label{Phi-def}
\begin{array}{rcl}
\displaystyle\Phi(x,y)_{a,i}^{a'}&=&\displaystyle\sigma_i\,q^{-\sigma_i a}\,
\big((x/y)^{-\sigma_i}\,\delta_{a,a'}+\delta_{a,a'-1}\big)\,,\\[.3cm]
\displaystyle\overline{\Phi}(x,y)_{a,i}^{a'}&=&\displaystyle\phantom{\sigma_i\,}q^{+\sigma_i a}\,
\big(\delta_{a,a'}-(x/y)^{\sigma_i}\,\delta_{a,a'-1}\big)\,,
\end{array}
\end{equation}
where $i=0,1$ and  $\sigma_0=+1$, $\sigma_1=-1$.
The above vectors satisfy important orthogonality relations
\begin{equation}\label{Phi-ort}
\sum_{a'} \Phi(x,y)_{a,i}^{a'}\,
\overline{\Phi}(x,y)_{a,j}^{a'}\;=\;[y/x]\,
\delta_{i,j}\;,\qquad
\sum_{i} \overline{\Phi}(x,y)_{a,i}^{a'}\,
\Phi(x,y)_{a,i}^{a''}\;=\;[y/x]\,\delta_{a',a''}\;.
\end{equation}
Using the explicit form of the two-by-two $L$-operator \eqref{Lop-new} 
it is easy to check that 
\begin{equation}\label{LL2}
\begin{tikzpicture}[baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.1, blue] (1.5,-1.5) -- (1.5,1.5) -- (2,1.7) -- (2.5,1.5) -- (2.5,-1.5) -- (2,-1.7) -- (1.5,-1.5);
\draw [-open triangle 45, thin] (1.5,1.5) -- (1.5,-1.5); \node [below] at (1.5,-1.5) {$x'$};
\draw [-open triangle 45, thin] (2.5,-1.5) -- (2.5,1.5); \node [above] at (2.5,1.5) {$x$};
\draw [thick] (1.5,0) -- (0.8,0); 
\draw [-latex, thick] (0,0) -- (0.8,0);
\node [left] at (0,0) {$i$}; \node [right] at (4,0) {$i'$};
\node [below] at (0.45,0) {$z$};
\node [above] at (3.3,0) {$z$};
\draw [thick, dashed] (2.5,0) -- (1.5,0); 
\draw [thick] (4,0) -- (3.5,0); 
\draw [-latex,thick] (2.5,0) -- (3.5,0);

\node [anchor=south east] at (1.5,0) {$\Phi$};
\node [anchor=north west] at (2.5,0) {$\overline\Phi$};
\node [below] at (2,-1) {$a$};
\node [above] at (2,1) {$a'$};
\node [left] at (-0.6,0)
      {$\big(\mathbb{L}(x,x',z)_{a}^{a'}\big)_{i,i'}
\;=\;\Phi(x',z)_{a,i}^{a'}\,
\overline{\Phi}(x,z)_{a,i'}^{a'}\;=\;$};
\end{tikzpicture}
\end{equation}
Remind, that the indices $(i,i')=(0,0)$ refer to the top left element 
in \eqref{Lop-new}.  
Here we have used the following graphical notations 
\begin{equation}\label{Psi2}
\begin{tikzpicture}[baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue,path fading=east] (1.5,0) -- (1.5,3) -- (3.0,3) -- (3.0,0) -- (1.5,0);
\draw [-open triangle 45, black, thin] (1.5,3) -- (1.5,0); \node [above] at (1.5,3) {$x$};
\draw [thick] (1.5,1.5) -- (0.80,1.5); \draw [-latex,thick] (0,1.5) --  
(0.8,1.5);
\node [left] at (0,1.5) {$i$};
\node [below] at (0.45,1.5) {$z$};
\node [anchor = south east] at (1.5,1.5) {$\Phi$};
\draw[dashed, thick] (1.5,1.5) -- (3.0,1.5); \node[right] at (2.0,2.5)
     {$a'$}; \node[right] at (2.0,0.5) {$a$};  
\node[left] at (-0.5,1.5) {$\displaystyle \Phi(x,z)^{a'}_{a,i}\;=$};
\end{tikzpicture}
\qquad\qquad \qquad
\begin{tikzpicture}[baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue, path fading=west] (1.5,0) -- (1.5,3) -- (0,3) -- (0,0) -- (1.5,0);
\draw [-open triangle 45, black, thin] (1.5,0) -- (1.5,3); \node [above] at (1.5,3) {$x$};
\draw [dashed, thick] (1.5,1.5) -- (0,1.5); \node[right] at (3,1.5) {$i$};
\node [above] at (2.3,1.5) {$z$};
\node [anchor = north west] at (1.5,1.5) {$\overline{\Phi}$};
\draw[thick] (3,1.5) -- (2.5,1.5); \draw [-latex, thick] 
(1.5,1.5) -- (2.5,1.5); 
\node[left] at (1,2.5) {$a'$}; \node[left] at (1,0.5) {$a$};
\node[left] at (0,1.5) {$\displaystyle
  \overline{\Phi}(x,z)_{a,i}^{a'}\;=\ $};
\end{tikzpicture}
\end{equation}
The are two types of spin variables in these pictures. 
The variables $a,a'\in{\mathbb Z}$ are the ``face spins'' 
assigned to the shaded faces, while $i,i'=\pm1$ are the edge
spins (like in vertex models) assigned to the directed solid lines in
the unshaded areas.  Moreover these lines (as well as their dashed
continuations into shaded areas) carry a spectral variable $z$. Another
spectral parameter $x$ is assigned to the directed vertical lines,
which separate the shaded and unshaded areas.

The vectors \eqref{sos-vertex} enter the so-called ``SOS-vertex''
correspondence and ``propagation through the vertex techniques''
invented by Baxter \cite{Baxter:1972,Bax73b} and further developed by
many authors including \cite{Andrews:1984,Date:1986,Krichever:1981}. 
The approach we use here is most similar to that of \cite{BKMS}, 
based on the factorized form of the $L$-operator \eqref{LL2}.   
The key relation of the SOS-vertex correspondence in our case reads
\begin{equation}\label{sos-duality}
\sum_{j_2,j_3=0,1} 
\; \widetilde{\mathcal{R}}^{(6v)}(y/z)_{i_2,i_3}^{j_2,j_3} \,
\Phi(x',z)_{d,j_3}^{b} \,
\Phi(x',y)_{b,j_2}^{a}\;=\;
\sum_{c=d}^a \; \Phi(x',y)_{d,i_2}^{c} \,
\Phi(x',z)_{c,i_3}^{a}\,
\mathcal{R}^{(6v)}(y/z)_{c-d,\;a-c}^{a-b,\;b-d}
\end{equation}
where $\mathcal{R}^{(6v)}$ and   $\widetilde{\mathcal{R}}^{(6v)}$ 
are given by \eqref{R6v} and \eqref{R6v-tilde}, respectively. 
This relation is presented graphically in
Fig.~\ref{fig-Rpsi-1}. Using the orthogonality \eqref{Phi-ort} one can
convert \eqref{sos-duality} into a similar relation containing the
$\overline{\Phi}_a^{a'}$ vectors instead of ${\Phi}_a^{a'}$,
\begin{equation}
\sum_{i_2,i_3=0,1} \; 
\overline{\Phi}(x,y)_{d,i_2}^{c}\,
\overline{\Phi}(x,z)_{c,i_3}^{a} \,
\widetilde{\mathcal{R}}^{(6v)}(y/z)_{i_2,i_3}^{j_2,j_3} 
\;=\;
\sum_{b=d}^a \; 
\mathcal{R}^{(6v)}(y/z)_{c-d,\;a-c}^{a-b,\;b-d} \,
\overline{\Phi}(x,z)_{d,j_3}^{b}\,
\overline{\Phi}(x,y)_{b,j_2}^{a}\,.
\end{equation}
\begin{figure}[ht]
\begin{center}
\begin{tikzpicture}[scale=.9,baseline=(current  bounding  box.center)]
\draw [-open triangle 45, thin] (2,3) -- (2,-1); \node [below] at
(2,-1) {$x'$}; 
\draw [thick] (2,0) -- (-2.,2); 
\draw [-latex, thick] (-2,2) -- (-1.5,1.75); 

\node[left] at (-2,2) {$i_3$};
\draw [dashed, thick] (2,0) -- (3,-0.5);
\node [above] at (-0,1.2) {$\widetilde{\cal R}^{(6v)}$};
\draw [black, ultra thick, fill] (0,1) circle [radius=0.09];
\node [anchor = south east] at (2.05,.2) {$\Phi$};
\draw [thick] (2,2) -- (-2,.0); 
\draw [-latex, thick] (-2,0) -- (-1.5,.25); 
\node[left] at (-2,0) {$i_2$};
\draw [dashed, thick] (2,2) -- (3,2.5);
\node [anchor = south east] at (2.05,1.95) {$\Phi$};
\node [above] at (1,1.5) {$j_2$};
\node [below] at (1,0.5) {$j_3$};
\node [above] at (-1.8,1.3) {$z$};
\node [below] at (-1.8,0.7) {$y$};
\node [right] at (2.25,-0.9) {$d$};
\node [right] at (2.25,1) {$b$};
\node [right] at (2.25,2.9) {$a$};
\draw [fill, opacity=0.25, blue,path fading=east] (2,-1) -- (2,3) -- (3,3.5) -- (3,-1.5) -- (2,-1);
\node [right] at (4,1) {$=$};
\end{tikzpicture}
\hskip 2mm
\begin{tikzpicture}[scale=.9,baseline=(current  bounding  box.center)]
\draw [-open triangle 45, thin] (0,3) -- (0,-1); \node [below] at (0,-1) {$x'$};
\draw [thick] (0,0) -- (-1,-0.5); 
\draw [-latex, thick] (-1,-0.5)--(-0.5,-0.25); 
\node[left] at (-1,2.5) {$i_3$};
\draw [dashed, thick] (0,0) -- (4,2);
\draw [thick] (0,2) -- (-1,2.5); 
\draw [-latex, thick]  (-1,2.5) -- (-0.5,2.25); 
\node[left] at (-1,-0.5) {$i_2$};
\draw [dashed, thick] (0,2) -- (4,0.);
\node [below] at (-0.7,2.3) {$z$};
\node [below] at (-0.7,-0.4) {$y$};
\node [anchor = south east] at (0,0) {$\Phi$};
\node [anchor = south east] at (0,2.2) {$\Phi$};
\node [right] at (0.4,1) {$c$};
\node [above] at (2,2) {$a$};
\node [below] at (2,0.) {$d$};
\node [right] at (3,1) {$b$};
\draw [fill, opacity=0.25, blue,path fading=east] (0,-1) -- (0,3) -- (4,3.5) -- (4,-1.5) -- (0,-1);
\node [above] at (2,1.2) {${\cal R}^{(6v)}$};
\draw [black, ultra thick, fill] (2,1) circle [radius=0.09];
\end{tikzpicture}
\end{center}
\caption{Graphical representation of the 
SOS-vertex duality \eqref{sos-duality} for the 6-vertex model.}
\label{fig-Rpsi-1}
\end{figure} 
 


Similarly to \eqref{sos-vertex} we need to define additional vectors
\begin{equation}\label{omega-def}
\begin{array}{rcl}
\displaystyle\Omega(x,y)_{a,i}^{a'}&=&\displaystyle\sigma_i\,q^{-\sigma_i a}\,
\big(\delta_{a,a'}+(y q/x)^{-\sigma_i}\,\delta_{a,a'-1}\big)\,,\\[.3cm]
\displaystyle\overline{\Omega}(x,y)_{a,i}^{a'}&=&\displaystyle\phantom{\sigma_i\,}q^{+\sigma_i a}\,
\big((y/x)^{\sigma_i}\delta_{a,a'}-q^{\sigma_i}\,\delta_{a,a'-1}\big)\,,
\end{array}
\end{equation}
satisfying the following orthogonality
relations
\begin{equation}\label{omega-ort}
\sum_{a} \Omega(x,y)_{a,i}^{a'}\,
\overline{\Omega}(x,y)_{a,i'}^{a'}\;=\;[y/x]\,
\delta_{i,i'}\;,\quad
\sum_{i} \overline{\Omega}(x,y)_{a',i}^{a}\,
\Omega(x,y)_{a'',i}^{a}\;=\;[y/x]\,\delta_{a',a''}\;.
\end{equation}
These vectors are represented graphically as 
\begin{equation}\label{Psi1}
\begin{tikzpicture}[baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue,path fading=east] (1.5,0) -- (1.5,3) -- (3,3) -- (3,0) -- (1.5,0);
\draw [-open triangle 45, black, thin] (1.5,0) -- (1.5,3); \node [above] at (1.5,3) {$x$};
\draw [thick] (1.5,1.5) -- (0.80,1.5); \draw [-latex,thick] (0,1.5) --  
(0.8,1.5);
\node [left] at (0,1.5) {$i$};
\node [below] at (0.45,1.5) {$z$};
\node [anchor = south east] at (1.5,1.5) {$\Omega$};
\draw[dotted, thick] (1.5,1.5) -- (3,1.5); \node[right] at (2,2.5) {$a'$}; \node[right] at (2,0.5) {$a$};
\node[left] at (-0.5,1.5) {$\displaystyle \Omega(x,z)_{a,i}^{a'}\;=$};
\end{tikzpicture}
\qquad\qquad \qquad
\begin{tikzpicture}[baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue,path fading=west] (1.5,0) -- (1.5,3) -- (0,3) -- (0,0) -- (1.5,0);
\draw [-open triangle 45, black, thin] (1.5,3) -- (1.5,0); \node
      [above] at (1.5,3) {$x$};
\draw [dotted, thick] (1.5,1.5) -- (0,1.5); \node[right] at (3,1.5) {$i$};
\node [above] at (2.3,1.5) {$z$};
\node [anchor = north west] at (1.5,1.5) {$\overline{\Omega}$};
\draw[thick] (3,1.5) -- (2.5,1.5); \draw [-latex, thick] 
(1.5,1.5) -- (2.5,1.5); 


\node[left] at (1,2.5) {$a'$}; \node[left] at (1,0.5) {$a$};
\node[left] at (0,1.5) {$\displaystyle
  \overline{\Omega}(x,z)^{a'}_{a,i}\;=\ $};
\end{tikzpicture}
\end{equation}
Note that here the directions of the
vertical lines, carrying the spectral variable $x$, 
are reversed in comparison with \eqref{Psi2}. 
 


Now, return to Fig.~\ref{fig-lattice} on page \pageref{fig-lattice}
and shade alternatively all faces of the 
medial lattice to cover those containing the white and black sites (the medial
lattice is a square lattice formed by thin rapidity lines). Then we
can regard the spins in our Ising-type model as face spins and update
the graphical notations \eqref{twotypes} for the edge Boltzmann weights
\begin{equation}
\input{fig2c_shade.tex}
\label{V-shade}
\end{equation}
where $V_x(n)$ is defined in \eqref{Vdef}.


Now we can present three easily verifiable Yang-Baxter relations
involving the SOS-vertex vectors \eqref{Psi2}, \eqref{Psi1} and the
edge weights \eqref{V-shade}. The first equation is 
\begin{equation}\label{psieq1}
V_{x/y}(a-b) \; \sum_{i} \overline{\Phi}(x,z)_{a,i}^{a'}\, \Phi(y,z)_{b,i}^{b'} \;=\;
\sum_{i} \overline{\Omega}(y,z)_{a,i}^{a'} \,\Omega(x,z)_{b,i}^{b'} \; V_{x/y}(a'-b')\;.
\end{equation}
It has the following graphical representation
\begin{equation}\label{peq1}
\begin{tikzpicture}[scale=0.65]
\draw [fill, opacity=0.2, blue, path fading=west] (-2,-1) -- (-0.5,-1) -- (0,0) -- (-2,4) -- (-4,4) -- (-2,-1);
\draw [fill, opacity=0.2, blue, path fading=east] (2,-1) -- (0.5,-1) -- (0,0) -- (2,4) -- (4,4) -- (2,-1);
\draw [-open triangle 45, black, thin] (2,4) -- (-0.5,-1); \node [above] at (2,4) {$y$};
\draw [-open triangle 45, black, thin] (0.5,-1) -- (-2,4); \node [above] at (-2,4) {$x$};
\draw [-latex, thick] (-1,2) -- (0,2); \draw [thick] (1,2) -- (-1,2);
\draw[dashed, thick] (-1,2) -- (-3.2,2);  \draw[dashed, thick] (1,2) -- (3.2,2); \node [below] at (0,2) {$z$};
\draw [blue!60!black, ultra thick, fill] (1.3,0) circle [radius=0.09]; \draw [blue!60!black, ultra thick, fill] (-1.3,0) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (1.3,0) -- (0,0); \draw[green!20!black, ultra thick] (0,0) -- (-1.3,0);
\node [left] at (-1.35,0) {$a$}; \node[right] at (1.35,0) {$b$};
\node [left] at (-2.1,3) {$a'$}; \node [right] at (2.1,3) {$b'$};
\node [anchor=south west] at (-1,2) {$\overline{\Phi}$};
\node [anchor=south east] at (1,2) {$\Phi$};
\node [right] at (4.5,1.5) {$=$};
\end{tikzpicture}
\hskip 10mm
\begin{tikzpicture}[scale=0.65,baseline=-40]
\draw [fill, opacity=0.2, blue,path fading=west] (-4,-2) -- (-2,-2) -- (0,2) -- (-0.5,3) -- (-2,3) -- (-4,-2);
\draw [fill, opacity=0.2, blue,path fading=east] (4,-2) -- (2,-2) -- (0,2) -- (0.5,3) -- (2,3) -- (4,-2);
\draw [-open triangle 45, black, thin] (0.5,3) -- (-2,-2) ; \node [below] at (-2,-2) {$y$};
\draw [-open triangle 45, black, thin] (2,-2) -- (-0.5,3); \node [below] at (2,-2) {$x$};
\draw [-latex, thick] (-1,0) -- (0,0); \draw [thick] (1,0) -- (-1,0);
\draw[dashed, thick] (-1,0) -- (-3.2,0);  \draw[dashed, thick] (1,0) -- (3.2,0); \node [above] at (0,0) {$z$};
\draw [blue!60!black, ultra thick, fill] (1.3,2) circle [radius=0.08]; \draw [blue!60!black, ultra thick, fill] (-1.3,2) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (1.3,2) -- (0,2); \draw[green!20!black, ultra thick] (0,2) -- (-1.3,2);
\node [left] at (-1.35,2) {$a'$}; \node[right] at (1.35,2) {$b'$};
\node [left] at (-2.1,-1) {$a$}; \node [right] at (2.1,-1) {$b$};
\node [anchor=north west] at (-1,0) {$\overline{\Omega}$};
\node [anchor=north east] at (1,0) {$\Omega$};
\end{tikzpicture}
\end{equation}
Note, that after a substitution of the definitions \eqref{Psi2}, \eqref{Psi1}
this relation uniquely determines the function $V_x(n)$, defined by
\eqref{Vdef}-\eqref{vanish}, up to an overall normalization.  The next
relation is 
\begin{equation}
\label{psieq2}
\sum_{b} V_{y/x}(a-b)\  \Phi(x,z)_{b,i}^{c} \  
\overline{\Omega}(y,z)_{b,i'}^{c}\;=\;
\sum_{b} \Phi(y,z)_{a,i}^{b} \ 
\overline{\Omega}(x,z)_{a,i'}^{b}\  
V_{y/x}(b-c)\;,
\end{equation}
It is represented graphically as
\begin{equation}\label{peq2}
\begin{tikzpicture}[scale=0.7,baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.25, blue, path fading=north] (0,1) -- (-1.5,4) -- (0,4.5) -- (1.5,4) -- (0,1);
\draw [fill, opacity=0.2, blue,path fading=south] (0,1) -- (1,-1) -- (0,-1.5) -- (-1,-1) -- (0,1);
\draw [-open triangle 45, black, thin] (1.5,4) -- (-1,-1); \node [above] at (1.5,4) {$y$};
\draw [-open triangle 45, black, thin] (-1.5,4) -- (1,-1); \node [above] at (-1.5,4) {$x$};
\draw [dashed, thick] (1,3) -- (0,3); \draw [dashed, thick] (0,3) -- (-1,3);
\draw[-latex, thick] (-2,3) -- (-1.3,3); \draw [thick] (-1.,3) -- (-2,3);  
\draw[-latex, thick] (1,3) -- (1.7,3); \draw [thick] (2.,3) -- (1,3); 
\node [left] at (-2,3) {$i$}; \node [right] at (2,3) {$i'$};
\node [above] at (-1.75,3) {$z$};
\node [above] at (1.6,3) {$z$};
\draw [blue!60!black, ultra thick, fill] (0,-0.3) circle [radius=0.09]; \draw [blue!60!black, ultra thick, fill] (0,2.3) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (0,2.3) -- (0,1); \draw [blue!60!black, ultra thick] (0,1) -- (0,-0.3);
\node [below] at (0,-0.35) {$a$}; \node[above] at (0,2.34) {$b$};
\node [above] at (0,3.35) {$c$};
\node [anchor=north east] at (-1,3) {$\Phi$};
\node [anchor=north west] at (1,3) {$\overline{\Omega}$};
\node [right] at (3.5,1.5) {$=$};
\end{tikzpicture}
\hskip 10mm
\begin{tikzpicture}[scale=0.7,baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue,path fading=north] (0,2) -- (-1,4) -- (0,4.5) -- (1,4) -- (0,2);
\draw [fill, opacity=0.25, blue,path fading=south] (0,2) -- (1.5,-1) -- (0,-1.5) -- (-1.5,-1) -- (0,2);
\draw [-open triangle 45, black, thin] (1,4) -- (-1.5,-1); \node [above] at (1,4) {$y$};
\draw [-open triangle 45, black, thin] (-1,4) -- (1.5,-1); \node [above] at (-1,4) {$x$};
\draw [dashed, thick] (1,0) -- (0,0); \draw [dashed, thick] (0,0) -- (-1,0);
\draw[-latex, thick] (-2,0) -- (-1.3,0); \draw [thick] (-1.,0) -- (-2,0);  
\draw[-latex, thick] (1,0) -- (1.7,0); \draw [thick] (2.,0) -- (1,0); 
\node [left] at (-2,0) {$i$}; \node [right] at (2,0) {$i'$};
\node [below] at (-1.75,0) {$z$};
\node [below] at (1.6,0) {$z$};
\draw [blue!60!black, ultra thick, fill] (0,0.7) circle [radius=0.09]; \draw [blue!60!black, ultra thick, fill] (0,3.3) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (0,3.3) -- (0,2); \draw [blue!60!black, ultra thick] (0,2) -- (0,0.7);
\node [below] at (0,0.66) {$b$}; \node[above] at (0,3.35) {$c$};
\node [below] at (0,-0.3) {$a$};
\node [anchor=south east] at (-1,0) {$\Phi$};
\node [anchor=south west] at (1,0) {$\overline{\Omega}$};
\end{tikzpicture}
\end{equation}
Note the summation over the spin $b$ corresponding to the internal
face. Finally, the third relation reads
\begin{equation}
\label{psieq3}
\sum_{b} V_{y/x}(b-a)\  \Omega(x,z)_{b,i}^{c} \  
\overline{\Phi}(y,z)_{b,i'}^{c}\;=\;
\sum_{b} \Omega(y,z)_{a,i}^{b} \ 
\overline{\Phi}(x,z)_{a,i'}^{b}\  
V_{y/x}(c-b)\;.
\end{equation}
It is depicted below
\begin{equation}\label{peq3}
\begin{tikzpicture}[scale=0.7,baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.25, blue,path fading=north] (0,1) -- (-1.5,4) -- (0,4.5) -- (1.5,4) -- (0,1);
\draw [fill, opacity=0.2, blue,path fading=south] (0,1) -- (1,-1) -- (0,-1.5) -- (-1,-1) -- (0,1);
\draw [-open triangle 45, black, thin] (-1,-1) -- (1.5,4); \node [above] at (1.5,4) {$y$};
\draw [-open triangle 45, black, thin] (1,-1) -- (-1.5,4); \node [above] at (-1.5,4) {$x$};
\draw [dashed, thick] (1,3) -- (-1,3);
\draw[-latex, thick] (-2,3) -- (-1.3,3); \draw [thick] (-1.,3) -- (-2,3);  
\draw[-latex, thick] (1,3) -- (1.7,3); \draw [thick] (2.,3) -- (1,3); 
\node [left] at (-2,3) {$i$}; \node [right] at (2,3) {$i'$};
\node [above] at (1.6,3) {$z$};
\node [above] at (-1.75,3) {$z$};
\draw [blue!60!black, ultra thick, fill] (0,-0.3) circle [radius=0.09]; \draw [blue!60!black, ultra thick, fill] (0,2.3) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (0,-0.3) -- (0,1); \draw [blue!60!black, ultra thick] (0,1) -- (0,2.3);
\node [below] at (0,-0.35) {$a$}; \node[above] at (0,2.34) {$b$};
\node [above] at (0,3.55) {$c$};
\node [anchor=north east] at (-1,3) {$\Omega$};
\node [anchor=north west] at (1,3) {$\overline{\Phi}$};
\node [right] at (3.5,1.5) {$=$};
\end{tikzpicture}
\hskip 10mm
\begin{tikzpicture}[scale=0.7,baseline=(current  bounding  box.center)]
\draw [fill, opacity=0.2, blue,path fading=north] (0,2) -- (-1,4) -- (0,4.5) -- (1,4) -- (0,2);
\draw [fill, opacity=0.25, blue,path fading=south] (0,2) -- (1.5,-1) -- (0,-1.5) -- (-1.5,-1) -- (0,2);
\draw [-open triangle 45, black, thin] (-1.5,-1) -- (1,4); \node [above] at (1,4) {$y$};
\draw [-open triangle 45, black, thin] (1.5,-1) -- (-1,4); \node [above] at (-1,4) {$x$};
\draw [dashed, thick] (1,0) -- (0,0); \draw [dashed, thick] (0,0) -- (-1,0);
\draw[-latex, thick] (-2,0) -- (-1.3,0); \draw [thick] (-1.,0) -- (-2,0); 
\draw[-latex, thick] (1,0) -- (1.7,0); \draw [thick] (2.,0) -- (1,0); 
\node [left] at (-2,0) {$i$}; \node [right] at (2,0) {$i'$};
\node [below] at (-1.75,0) {$z$};
\node [below] at (1.6,0) {$z$};
\draw [blue!60!black, ultra thick, fill] (0,0.7) circle [radius=0.09]; \draw [blue!60!black, ultra thick, fill] (0,3.3) circle [radius=0.09];
\draw [-latex, green!20!black, ultra thick] (0,0.7) -- (0,2); \draw [blue!60!black, ultra thick] (0,2) -- (0,3.3);
\node [below] at (0,0.66) {$b$}; \node[above] at (0,3.35) {$c$};
\node [below] at (0,-0.5) {$a$};
\node [anchor=south east] at (-1,0) {$\Omega$};
\node [anchor=south west] at (1,0) {$\overline{\Phi}$};
\end{tikzpicture}
\end{equation}
With the vectors \eqref{Phi-def} and \eqref{omega-def} regarded as an
input, each of the equations \eqref{psieq2} and \eqref{psieq3}
can be considered as an equation for the weight function $V_x(n)$.
Similarly to \eqref{psieq1}, these equation also
lead to the function $V_x(n)$, defined by
\eqref{Vdef}-\eqref{vanish} (up to an overall normalization). All
these equations can be viewed as (Yang-Baxter type) exchange
relations, suitable for defining $Z$-invariant lattices models
\cite{Baxter:1978xr}. The models of this type are formulated 
on irregular lattices
formed by an intersection of straight lines (with at most two lines
intersecting at one point). Their partition function remain unchanged
upon a parallel transport of individual lattice lines. In our
case the lattices shown in \eqref{peq1}, \eqref{peq2} and \eqref{peq3}
are the simplest examples of such lattices. For instance, the lattice
on the left side of \eqref{peq1} is formed by three straight lines
associated with the spectral variables $x$, $y$ and $z$. There are four spins
$a$, $a'$, $b$, $b'$, assigned to external faces and one spin $i$
assigned to the internal edge on the $z$-line. The model contains 
the three-spin interactions, with the Boltzmann 
weights given by \eqref{Phi-def},
\eqref{omega-def}, as well as the two-spin interaction with the
weights \eqref{V-shade}.
The LHS of
\eqref{psieq1} is the partition function of this model for the left lattice
of \eqref{peq1} with fixed boundary spins on the external faces. 
Eq. \eqref{psieq1} states that 
this partition function remains unchanged by
moving the $z$-line down through the intersection point of the other
two lines and transforming the lattice to that shown in the right side of
\eqref{peq1}.  A similar interpretation exists for the other two
relations \eqref{psieq2} and \eqref{psieq3}.

We are now ready to give a simple visual proof of the Yang-Baxter
equation \eqref{RLL-irf}. Substituting the factorized expressions 
\eqref{def3} and \eqref{LL2} and using the above graphical notations 
one can represent this equation as shown in Fig.~\ref{RLL-new-v2}. 
\input RLL_new_v2.tex
Consider the diagram on the left side on this figure. Let us move 
the horizontal $z$-line down through the intersection points of the
other lines and then consequently use: (i) the relation \eqref{peq1},
(ii) the relations \eqref{peq2} and \eqref{peq3} (in any order) and,
finally, (iii) the relation \eqref{peq1} again. In this way the left diagram
in Fig.~\ref{RLL-new-v2} is transformed to the right one, thus proving
\eqref{RLL-irf}. The above calculation completes our proof of the
linear relation \eqref{fourier} 
between $R$-matrices arising from the ``star'' and
``box'' diagrams. 
Relations of this
type commonly appear as (self-) duality transformations in the Ising-type 
lattice model \cite{Wu:1982,Savit:1980}. 
\section{Functional relations for the transfer matrices\label{FR-section}}
\subsection{Definition of the
 ${\bf T}$- and ${\bf Q}$-matrices.\label{TQ-section}}
To discuss the functional relations for the transfer matrices it is 
convenient
to also use the equivalent vertex model, introduced in
Sect.~\ref{irf-vertex}. 
\begin{figure}[h]
\begin{center}
\input fig-irf-vert.tex
\end{center}
\caption{The arrangement of spins for the transfer matrices in the IRF
  and vertex formulations of the model, defined by \eqref{tmat} and
  \eqref{t-vert}, respectively. The correspondence between the spin
  variables is given by \eqref{index}. Note also, that the 
boundary conditions \eqref{quasi1} imply 
that $k_{N+1}=b_{N+1}-a_{N+1}=b_{1}-a_{1}=k_1$. \label{t-vertex}}
\end{figure}
Fig.~\ref{t-vertex} shows the correspondence
between the corner spins of the IRF model and the edge 
indices of the vertex model, 
\begin{equation}\label{index}
i_n=a_{n+1}-a_n\,,\qquad  j_n=b_{n+1}-b_n\,,\qquad
    k_n=b_n-a_n\,,\qquad n=1,2,\ldots N\,,
\end{equation}
which take arbitrary non-negative integer values.
For the quasi-periodic
boundary condition (involving the horizontal field factor $\omega_h^{k_1}$ at the first
column of the lattice) the elements of the transfer matrix read
\begin{equation}\label{t-vert}
\begin{array}{l}
\displaystyle
\big({\bf T}^{(+)}(y,y'|x_N,x_N',\ldots,x_1,x'_1)\big)_{i_N,i_{N-1},
\ldots,i_1}^{j_N,j_{N-1},\ldots,j_1}=\rho(y,y')\,
\sum_{k_1,\ldots,k_N=0}^\infty\,
 \omega_h^{k_1} \,{\mathcal R}(x_N,x'_N,y,y')_{i_N,k_1}^{j_N,k_N}\,\times\\[.8cm]
\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\times
\,
{\mathcal R}(x_{N-1},x'_{N-1},y,y')_{i_{N-1},k_N}^{j_{N-1},k_{N-1}}\ 
\cdots\ 
{\mathcal R}(x_1,x'_1,y,y')_{i_1,k_2}^{j_1,k_1}\,,
\end{array}
\end{equation}
where the sum is taken over all nonnegative integer values of
$k_1,\ldots,k_N$. Actually, there is only a one-dimensional summation
with $k_1=0,1,\ldots,\infty$, while the rest of $k$'s is uniquely
determined by the values of $i$'s, $j$'s and $k_1$ due to 
the conservation laws \eqref{ar-cons} at each vertex.  The superscript
$\scriptstyle{(+)}$ in the notation for the transfer matrix indicates that the sum
over $k_1$ goes over an infinite interval. Using the correspondence 
\eqref{index} and the relation \eqref{R-def} it is easy to see that
the expression inside the sum in \eqref{t-vert} exactly coincides with
the IRF transfer matrix \eqref{tmat} (with $\omega_v=1$). 




The transfer matrix \eqref{t-vert} is an operator acting in the
``quantum'' space, formed by the direct
product 
\begin{equation}\label{qspace}
{\mathcal H}^{\rm (vertex)}
={\mathbb C}^{(\infty)}_N\otimes{\mathbb C}^{(\infty)}_{N-1}\otimes\cdots
\otimes   {\mathbb C}^{(\infty)}_1\,,
\end{equation}
of $N$ identical infinite-dimensional vector spaces ${\mathbb C}^{(\infty)}$, 
each of which is spanned by the basis vectors $|j\rangle$ with
$j=0,1,2,\ldots,\infty$.  
The lower indices in \eqref{qspace} numerate the columns of the lattice, 
as in Fig.~\ref{t-vertex}.
Below we will
regard the variables $(x_N,x'_N,\ldots,x_1,x'_1)$ as fixed and 
will write the transfer matrix simply as ${\bf T}^{(+)}(y,y')$ 
assuming the
dependence on fixed variables implicitly. 
The Yang-Baxter equation \eqref{Rybe} implies that the transfer matrices ${\bf
  T}^{(+)}(y,y')$ with different values of $y,y'$ form a commuting
family,
\begin{equation}
\big[{\bf T}^{(+)}(y,y'),{\bf T}^{(+)}(z,z')\big]=0\,,\qquad 
\forall \qquad y,y',z,z'\,,\label{tcomm2}
\end{equation}
Moreover, they preserve the number of particles $M$, 
previously introduced by \eqref{down}.
Indeed, the matrix elements \eqref{t-vert} vanish, unless
\begin{equation}\label{down2}
i_1+i_2+\cdots+i_N=j_1+j_2+\cdots+j_N\equiv M\,,\qquad M=0,1,2\ldots \,.
\end{equation}
So, the transfer matrices act invariantly in each subspace of
\eqref{qspace} with a fixed particle number $M=0,1,2,\ldots,\infty$. 
It is convernient then to define diagonal matrices 
\begin{equation}
w\;=\; q^{2M}\, \prod_{i=1}^N\, \frac{x_i}{x_i'}\,,\qquad 
Z_0=\sum_{k=0}^\infty\,(w\,\omega_h)^k=\big(1-\omega_h\,w\big)^{-1}
\,,\label{zw-def}
\end{equation}
which also belong to the commuting family \eqref{tcomm2}. The sum in
$Z_0$ converges provided the horizontal field parameter $\omega_h$ 
is such that 
\begin{equation} \label{ow1}
|\omega_h w|<1\,.
\end{equation}

At this point it is worth noting that the
definition of the $R$-matrix \eqref{R-def2} only involves a finite
summation, so it is well defined for arbitrary values of its parameters, in
particular, for an arbitrary value of the ``anisotropy'' parameter
$q$. However the convergence of the infinite sum in the definition of
the transfer matrix \eqref{t-vert} requires a special investigation. 
Plausibly, for $|q|<1$ or $|q|=1$ the restriction \eqref{ow1} for the field
factor is sufficient for convergence of the sum \eqref{t-vert}. On the
other hand, the
considered model is extremely general, so that
a more detailed analysis is, of course, needed for considerations of
its particular cases. For simplicity, in what follows, we assume
$|q|<1$.   

The normalization factor in \eqref{t-vert} is chosen as
\begin{equation}\label{t-norm}
\rho(y,y')=(y'/y)^M\,\prod_{n=1}^N\frac{\big(\,(y/x'_n)^2\,;q^2\big)_\infty}
{\big(\,(y'/x'_n)^2\,;q^2\big)_\infty}\,,\
\end{equation}
where the $q$-products are defined in \eqref{qprod0}.
Then, with an account of \eqref{R-def2} one concludes\footnote{%
The $R$-matrix \eqref{R-def2} is a Laurent
polynomial of the spectral variable $y'$ and a meromorphic function in
$y$, having poles at $y^2={x'}^2 q^{-2n}$,\  $n=0,1,2,\ldots$.} that
the transfer matrix \eqref{t-vert} is a meromorphic function in
the variables $y$ and $y'$. The reason for taking the normalization
\eqref{t-norm} will be clarified below. 
Next, define the factor
\begin{equation}\label{phi-def}
\varphi(y,y')\,=\,
\prod_{n=1}^N\frac{\big(\,(y'/x'_n)^2\,;q^2\big)_\infty}
{\big(\,(y/x_n)^2\,;q^2\big)_\infty}\;,
\end{equation}
and introduce two new transfer matrices 
\begin{equation}\label{Q-def}
{\bf Q}(y')=Z^{-1}_0\,\lim_{y\to 0} 
\,\varphi(y,y')\,{\bf T}^{(+)}(y,y')\,,
\qquad
\overline{{\bf Q}}(q^{-1}\,y)=Z^{-1}_0\,
\lim_{y'\to 0} \,\varphi(y,y')\,{\bf T}^{(+)}(y,y')\,,
\end{equation}
which, as we shall see below, play the role of the Baxter ${\bf
  Q}$-matrices \cite{Baxter:1972}. 
Evidently, they belong to the same commuting family \eqref{tcomm2},
\begin{equation}\label{tq-comm}
[{\mathbf T}^{(+)}(y,y'),{\bf Q}(z)]=[{\mathbf T}^{(+)}(y,y'),\overline{{\bf Q}}(z)]=
[{\bf Q}(z),{\bf Q}(z')]=[\overline{{\bf Q}}(z),\overline{{\bf Q}}(z')]=
[{{\bf Q}}(z),\overline{{\bf Q}}(z')]=0\,,
\end{equation}
where $y,y',z,z'$ take arbitrary values. 




Remind, that ${\mathbf T}^{(+)}(y,y')$ as well as ${\bf Q}(y')$ 
and $\overline{{\bf Q}}(y)$ implicitly depend on the set of fixed variables 
$(x_N,x'_N,\ldots,x_1,x'_1)$. Note also, that if a pair these
variables is constrained as  
\begin{equation}\label{x-const}
x_n'=x_n\, q^{2s_n}\,,\qquad 2s_n\in {{\mathbb Z}}_{\ge0}
\end{equation}
for some particular $n$, then the above transfer matrices admit a
$(2s_n+1)$-dimensional reduction in the $n$-th factor of the 
product \eqref{qspace}.

\subsection{Functional relations}
The elements of the ${\bf Q}$-matrices \eqref{Q-def} 
can be written explicitly by taking required limits of
\eqref{t-vert}. The details are presented in Appendix~\ref{appD}. 
It is worth noting that in case of $U_q(\widehat{sl}(2))$ 
considered here the algebraic origin of the ${\bf Q}$-matrices 
(or the ${\bf Q}$-operators) is very well understood \cite{Bazhanov:1996dr,
Bazhanov:1998dq}.  
Their main algebraic properties can be concisely expressed by a single
factorization relation, which in our case takes the form
\begin{equation}\label{factor}
\varphi(y,y')\,{\mathbf T}^{(+)}(y,y')=
Z_0\,{\bf Q}(y')\,\overline{{\bf Q}}(q^{-1}\,y)
\,,
\end{equation}
derived in Appendix~\ref{appD}. Mention also the normalization conditions 
\begin{equation} \label{Q-norm}
{\bf Q}(0)=\overline{{\bf Q}}(0)=1\,,
\end{equation}
where the ${\bf Q}$-operators reduce to the identity operator. 

Consider now the case when the spectral variables $y$ and $y'$ are
related as
\begin{equation}\label{findim}
y'=y\, q^{2s}\,,\qquad 2s\in{\mathbb
  Z}_{\ge0}\,. 
\end{equation}
This is precisely the finite-dimensional 
reduction case \eqref{twos2}-\eqref{inf}, considered in 
Sect.~\ref{red} (see also
\eqref{RR-rel}). This reduction implies that the sum in \eqref{t-vert} splits
into two parts: 
\begin{enumerate}[(i)]
\item
the first part, where all the indices
$k_1,\ldots,k_N\in\,{\mathcal I}_s=\{0,1,\ldots,2s\}$
take the finite set of
values. Assuming the relation \eqref{findim} we denote this part as 
${\bf T}_s(y)$ without the superscript
$\scriptstyle(+)$. Explicitly, for $2s\in{\mathbb
  Z}_{\ge0}$ one obtains 
\begin{equation}\label{ts-def}
\begin{array}{l}
\displaystyle
\big({\bf T}_s(y)\big)_{i_N,i_{N-1},
\ldots,i_1}^{j_N,j_{N-1},\ldots,j_1}=\rho(y,q^{2s}\,y)\,
\sum_{k_1,\ldots,k_N\in \,{\mathcal I}_s}\,
 \omega_h^{k_1}
 \,{\mathcal R}(x_N,x'_N,y,q^{2s}\,y)_{i_N,k_1}^{j_N,k_N}\,\times
\\[.8cm]
\displaystyle\ \qquad\qquad\qquad\qquad \times
\,
{\mathcal R}(x_{N-1},x'_{N-1},y,q^{2s}\,y)_{i_{N-1},k_N}^{j_{N-1},k_{N-1}}
\cdots
{\mathcal R}(x_1,x'_1,y,q^{2s}\,y)_{i_1,k_2}^{j_1,k_1}\,,
\end{array}
\end{equation}
\item
the second part where all the indices 
$k_1,\ldots,k_N\in
\overline{\mathcal I}_s=\{2s+1,2s+2,\ldots,\infty\}$
take the infinite set
of values. Using \eqref{twos2inf} one can reduce this part to ${\bf
  T}^{(+)}(y' q,y q^{-1})$. 
\end{enumerate}
In this way one obtains
\begin{equation}
\varphi(y,y')\,{\bf T}^{(+)}(y,y')=
\varphi(y,y')\,{\bf T}_s(y)+(\omega_h\,w)^{2s+1}\, 
\varphi(y'q, \,y q^{-1})\, {\bf
  T}^{(+)}(y'q,\,yq^{-1})\,,\qquad y'=y q^{2s}\,,
\end{equation}
where the relation \eqref{findim} is assumed.
The above reduction for the
transfer matrices is well known, see \cite{Bazhanov:1998dq}. 
Combining the last
equation with \eqref{factor} one gets
\begin{equation}\label{T-QQ}
Z_0^{-1}\,\varphi(y,yq^{2s})\, {\bf T}_s(y)={\bf Q}(y')
\overline{{\bf Q}}(y\,q^{-1})
-\big(\omega_h\,w\big)^{2s+1}\,
{\bf Q}(y\,q^{-1})
\overline{{\bf Q}}(y')\,.
\end{equation}


From \eqref{spin0} it follows that ${\bf T}_0(y)\equiv 1$. Thus,
setting $s=0$ in \eqref{T-QQ} one obtains the quantum Wronskian
relation for the ${\bf Q}$-operators 
\begin{equation}\label{qwrons}
Z_0^{-1}\,\varphi(y,y)={\bf Q}(y)
\overline{{\bf Q}}(y\,q^{-1})
-\big(\omega_h\,w\big)\,
{\bf Q}(y\,q^{-1})
\overline{{\bf Q}}(y)\,.
\end{equation}
Next, setting $s=1/2$ in \eqref{T-QQ} one gets
\begin{equation}
Z_0^{-1}\,\varphi(y,yq)\, {\bf T}_{\frac{1}{2}}(y)={\bf Q}(yq)
\overline{{\bf Q}}(y\,q^{-1})
-\big(\omega_h\,w\big)^2\,
{\bf Q}(y\,q^{-1})
\overline{{\bf Q}}(y q)\,.
\end{equation}
Combinig the last relation with \eqref{qwrons} one obtains the famous
$TQ$-equations \cite{Baxter:1972}
\begin{equation}\label{tqeq}
\begin{array}{l}
{\bf T}_{\frac{1}{2}}(y){\bf Q}(y)=f(y)\,{\bf Q}(yq)
+(\omega_h\,w)\,g(y)\,{\bf Q}(yq^{-1})\,, 
\\[.4cm]
{\bf T}_{\frac{1}{2}}(y)\overline{{\bf Q}}(y)=(\omega_h\,w)\,
f(y)\,\overline{{\bf Q}}(yq)
+g(y)\,\overline{{\bf Q}}(yq^{-1})\,, 
\end{array}
\end{equation}
where
\begin{equation}\label{fg-def}
f(y)=\frac{\varphi(y,y)}{\varphi(y,qy)}=\prod_{n=1}^N
\big(1-(y/{x'_n})^{2}
\big)\,,\qquad
g(y)=\frac{\varphi(qy,qy)}{\varphi(y,qy)}=\prod_{n=1}^N
\big(1-(y/{x_n})^{2}
\big)\,.
\end{equation}

To conclude our derivation of the functional relations  
note, that seemingly artificial normalization factors \eqref{t-norm}
and \eqref{phi-def}, entering the definitions of the transfer matrices
and ${\bf Q}$-operators were specially chosen to keep the functional
relations in the standard form, with a minimal number of the scalar
factors. If necessary (for instance, when $|q|=1$) the above  
normalization factors could be removed and the functional
relations then be easily modified. Finally, note that essentially the
same ${\bf Q}$-operators were defined by Mangazeev in
\cite{Mangazeev:2014gwa,Mangazeev:2014bqa} (see Appendix~\ref{appC}
for more details). 

\subsection{Six-vertex model}

By construction the eigenvalues of
$\varphi(y,y')\,{\mathbf T}^{(+)}(y,y')$, ${\bf Q}(y')$ and $\overline{{\bf Q}}(y)$
are polynomial in the
variable $(y')^2$, and, in general, meromorphic functions of $y^2$ 
with poles in the finite part of the complex plane.
An important exception is the finite-dimensional reduction case when
the constraint \eqref{x-const} holds for all columns of the lattice 
$m=1,2,\ldots,N$. 
For an illustration consider the case of the six-vertex model with
two-dimensional quantum space representations for all sites of the
chain; that is possible if the spectral parameters obey the relations
\begin{equation}
x_n'=q x_n\,,\qquad s_n=\textstyle{\frac{1}{2}}\,,\qquad m=1,\ldots,N\,.
\end{equation}
Using the definitions \eqref{t-vert}, \eqref{zw-def} and \eqref{t-norm}
one obtains
\begin{equation}
{\bf T}_{\frac{1}{2}}(y)=z^{-1}\,{\bf
  T}^{(6v)}(q^{-\frac{1}{2}}\,y)\,,\qquad
z=(\omega_h w)^{-\frac{1}{2}}\,,
\end{equation}
where
\begin{eqnarray}
\big({\bf T}^{(6v)}(q^{-\frac{1}{2}}\,y)\big)_{i_N,i_{N-1},
\ldots,i_1}^{j_N,j_{N-1},\ldots,j_1}&=&\\[.4cm]
&&\kern-11em=(\omega_h)^{-\frac{1}{2}}
\sum_{k_1,\ldots,k_N=0,1}
 \omega_h^{k_1} \,{\mathcal R}(x_N/y)_{i_N,k_1}^{j_N,k_N}
\,
{\mathcal R}(x_{N-1}/y)_{i_{N-1},k_N}^{j_{N-1},k_{N-1}}\ 
\cdots\ 
{\mathcal R}(x_1/y)_{i_1,k_2}^{j_1,k_1}\,,
\end{eqnarray}
where ${\mathcal R}(\lambda)$ is given by \eqref{R6v} with
$\rho(\lambda)=q^{-\frac{1}{2}}\,\lambda^{-1}$. 

To present the expressions for the ${\bf Q}$-operators \eqref{Q-def}
for this case consider the $q$-oscillator algebra 
\begin{equation}\label{qosc}
[{\mathcal H},{\mathcal E}_\pm]=\pm2\,{\mathcal E}_\pm\,,\qquad q\, {\mathcal E}_+{\mathcal E}_--
q^{-1}\,{\mathcal E}_-{\mathcal E}_+=q-q^{-1}
\end{equation}
and its two representations $\rho_\pm$, defined by the action on the vectors 
$|j\rangle$, $j=0,1,2,\ldots,\infty$,
\begin{equation}\label{qreps}
\rho_\pm({\mathcal H})|j\rangle=\mp 2j\,|j\rangle\,,\qquad
\rho_\pm({\mathcal E}_\pm)|j\rangle=q^{\mp\frac{1}{2}}
(q^{-j}-q^{j})\,|j-1\rangle\,,\qquad
\rho_\pm({\mathcal E}_\mp)|j\rangle=\mp q^{\mp (j+ 1/2)}|j+1\rangle\,.
\end{equation}
Define operators 
\begin{equation}\label{Apm-def}
{\bf A}_\pm(y)=Z^{-1}_0\,{\rm Tr}_{\rho_\pm}\Big\{\omega_h^{\mp{\mathcal H}/2}\,
L^{(N)}_\pm(y/x_N)\, L^{(N-1)}_\pm(y/x_{N-1})\cdots L^{(1)}_\pm(y/x_1)\Big\}\,,\qquad {\bf A}_\pm(y)=1+O(y^2)\,,
\end{equation}
where 
\begin{equation}\label{Lpm-def}
L_+(y)=\begin{pmatrix}
q^{\frac{{\mathcal H}}{2}}&y\,{\mathcal E}_-\\
y\, {\mathcal E}_+& \ \ q^{-\frac{{\mathcal H}}{2}}-q^{-1}\,y^2\,q^{\frac{{\mathcal H}}{2}}
\end{pmatrix}\,,\qquad
L_-(y)=\begin{pmatrix}
q^{-\frac{{\mathcal H}}{2}}-q\,y^2\,q^{\frac{{\mathcal H}}{2}}\ \ &y\,{\mathcal E}_+\\[.3cm]
y\,{\mathcal E}_-& q^{\frac{{\mathcal H}}{2}}
\end{pmatrix}
\end{equation}
are two by two matrices acting the quantum space, whose entries are
elements of the $q$-oscillator algebra \eqref{qosc}. 

Using \eqref{Qelem}, \eqref{Qelem4} and \eqref{Qelem5} from
Appendix~\ref{appD} it is not difficult to show that 
\begin{equation}
{\bf Q}(y)={\bf A}_+(q^{-\frac{1}{2}}\,y)\,,
\qquad
\overline{{\bf Q}}(y)={\bf A}_-(q^{-\frac{1}{2}}\,y)\,.
\end{equation}
The equations \eqref{tqeq} can now be rewriten a
\begin{equation}
{\bf T}^{(6v)}(y) {\bf A}_\pm(y)=z^{\pm1} \,g(q^{-\frac{1}{2}}\,y)\, 
{\bf A}_\pm(q\,y)+
z^{\mp1} \,g(q^{+\frac{1}{2}}\,y)\, {\bf A}_\pm(q^{-1}\,y)\,,\qquad
z=(\omega_h w)^{-\frac{1}{2}}\,,
\end{equation}
Note, that apart from a trivial change $y\to y^2=\zeta$ of the arguments of the
transfer matrices this relation exactly coincides with Eq.(3.4) of 
ref.~\cite{Bazhanov:2020fbp}. 
\section{The homogeneous case}
Consider the column homogeneous case 
\begin{equation}\label{homo}
x_1=x_2=\ldots=x_N\equiv x\,,\qquad 
x_1'=x_2'=\ldots=x_N'\equiv x'\,,\qquad x'=q^{2s}\,x\,,
\end{equation}
so that we can abbrevate the arguments of the transfer matrix \eqref{t-vert} 
simply as ${\bf T}^{(+)}(y,y'\,|\,x,x')$.
Obviously, the homogeneous model possesses the translational invariance
\begin{equation}\label{TK-com}
\big[{\bf T}^{(+)}(y,y'\,|\,x,x'),{\bf K}\big]=0\,,\qquad
\end{equation}
where ${\bf K}$ is the one-site  translation operator,
whose matrix elements are given by
\begin{equation}\label{K-def}
\big({{\bf K}}\big)_{i_{N}\, i_{N-1}\,\ldots
\,\,  i_1}^{j_{N} \,j_{N-1}\,\ldots \,j_1}=(\omega_h)^{i_1}\,
\delta_{i_N}^{j_{N-1}}\,\ldots\,
\delta_{i_2}^{j_1}\,\delta_{i_1}^{j_N} \,.
\end{equation}
In writing this formula we have taken 
into account the quasi-periodic boundary conditions, involving
the horizontal field factor at the first column of the lattice. 

\subsection{Hamiltonians}
Apart from the translation operator the transfer matrix 
commuting family now contains spin chain
Hamiltonians that are given by a sum of terms, each of which acts
non-trivially only in two consecutive factors of the product \eqref{qspace}.
First consider a generic case, when $2s$ is not taking non-negative
integer values. 
Using the definition \eqref{t-vert} and \eqref{R-def2} 
it is not difficult to check that 
\begin{equation}\label{T-expan}
\Big\{\rho(y,y')^{-1}\,{\bf T}^{(+)}(y,y'\,|\,x,x')\Big\}{\Big|_{y=x e^\varepsilon,\,
y'=x'e^{\varepsilon'}}}
={{\bf K}}\Big(1+\varepsilon\,{\bf H}^{(1)}+
\varepsilon'\,{\bf H}^{(2)}\Big)+O(\varepsilon^2,{\varepsilon'}^2,
\varepsilon \varepsilon')
\end{equation}
where the Hamiltonians ${\bf H}^{(1)}$ and  ${\bf H}^{(2)}$ are
defined as  
\begin{equation}\label{hamdef1}
\big({\bf H}^{(\ell)})_{i_{N}\, i_{N-1}\,\ldots
\,\,  i_1}^{j_{N} \,j_{N-1}\,\ldots \,j_1}=
\sum_{m=1}^{N} \omega_h^{(j_N-i_N)\delta_{m,N}}\,\,
\delta_{i_N}^{j_N}\cdots\,\delta_{i_{m+2}}^{j_{m+2}}\,\,
\big(H^{(\ell)}\big)_{i_{m+1},i_m}^{j_{m+1},j_m}\,\,
\delta_{i_{m-1}}^{j_{m-1}}\cdots\,\delta_{i_{1}}^{j_{1}}\,,
\end{equation}
with $\ell=1,2$ and 
\begin{equation}\label{hamdef2}
\big(H^{(1)}\big)_{i\,,\,i'}^{j\,,\,j'}\;=\;
\delta_{i}^{j}\,\delta_{i'}^{j'} \,F_{i}(s) \;+\;
2\delta_{i+i',j+j'} \,\frac{\theta(i>j)}{[q^{i-j}]} \,\frac{V_{q^{-2s}}(j)}{V_{q^{-2s}}(i)}\,,\qquad \big(H^{(2)}\big)_{i\,,\,i'}^{j\,,\,j'}\;=\;
\big(H^{(1)}\big)_{i'\,,\,i}^{j'\,,\,j}\,.
\end{equation}
Here we use the notations 
\begin{equation}
\theta(a>b)=\left\{\begin{array}{ll} 1, 
 &\mbox{if\ } a>b\,,\\[.3cm]
0,& \mbox{otherwise}\,,\end{array}\right.
\end{equation} 
and 
\begin{equation}
F_i(s)\;=\;-\sum_{k=0}^{i-1} \frac{\displaystyle q^{k-2s}+q^{-k+2s}}{\displaystyle q^{k-2s}-q^{-k+2s}}\;,\qquad
[q^n]\;=\;q^n-q^{-n}\;,
\end{equation}
where $i=0,1,2,\ldots,\infty$.
The weight function $V_x(n)$, appearing in \eqref{hamdef2}, is defined in
\eqref{Vdef}. 
Moreover, the boundary conditions imply that $i_{N+1}=i_1$ and  
$j_{N+1}=j_1$. 

It is worth reminding that the transfer matrix ${\bf T}^{(+)}(y,y'\,|\,x,x')$ has two spectral parameters $y$ and $y'$. This is why there are two independent commuting Hamiltonians in the homogeneous case. 


\subsection{{\bf Q}-operators}
For the homogeneous case it is convenient to define 
different ${\bf Q}$-operators
\begin{equation}\label{Qnew}
{\bf Q}^{(h)}(y')\;=\;\rho(y,y')^{-1}{\mathbf T}^{(+)}(y,y')|_{y=x}\;,\quad \overline{{\bf Q}}^{(h)}(q^{-1} y)\;=\;\rho(y,y')^{-1}{\mathbf T}^{(+)}(y,y')|_{y'=x}\;,
\end{equation}
with matrix elements,
\begin{eqnarray}
{\bf Q}^{(h)}(y')_{i_N i_{N-1} \dots i_1}^{j_N j_{N-1} \dots j_1}&=&
\sum_{k_1,\dots, k_N} 
\omega_h^{k_1} \prod_{n=1}^N \delta_{i_n+k_{n+1},j_n+k_n}^{}
\frac{V_{x/y'}(k_n) V_{y'/x'}(i_n-k_n)}{V_{x/x'}(i_n)}\;,\\[.4cm]
\overline{{\bf Q}}^{(h)}(q^{-1}y)_{i_N i_{N-1} 
\dots i_1}^{j_N j_{N-1} \dots j_1}&=&
\sum_{k_1,\dots, k_N} 
\omega_h^{k_1} \prod_{n=1}^N \delta_{i_{n}+k_{n+1},j_n+k_{n}}^{}
\frac{V_{x/x'}(j_n) V_{y/x}(k_n)}{V_{y/x'}(j_n+k_{n})}\;.
\end{eqnarray}
They satisfy the factorization relation (similar to \eqref{factor}) 
\begin{equation}
\rho(y,y') \,\overline{{\bf Q}}^{(h)}(q^{-1} y) \,
{\bf Q}^{(h)}(y')\;=\;{\mathbf T}^{(+)}(y,y')\,.
\end{equation}
and normalization conditions
\begin{equation}
{\bf Q}^{(h)}(x' e^{\varepsilon'})={\bf K}\big(1+\varepsilon' {\bf H}^{(2)}+
O({\varepsilon'}^2)\big)
\,,\qquad
{\bf Q}^{(h)}(x\,\mbox{e}^{-\varepsilon}) = 1 + \varepsilon {\bf H}^{(1)}+O(\varepsilon^2)\,,
\end{equation}
\begin{equation}
\overline{{\bf Q}}^{(h)}(q^{-1} x e^{\varepsilon})=
1+\varepsilon {\bf H}^{(1)}+O(\varepsilon^2)\;.
\end{equation}
where the Hamiltonians are defined in \eqref{hamdef1}, \eqref{hamdef2} and 
``$1$'' stands for the identity operator. We have added the
superscript $(h)$ to the notations for the new ${\bf Q}$-operators
\eqref{Qnew} to distinguish them from the previously defined operators
\eqref{Q-def}, normalized as \eqref{Q-norm}. The latter can, of
course, be specialized for the homogeneous case \eqref{homo} and
connected with \eqref{Qnew},
\begin{equation}\label{Qhnorm1}
{\bf Q}^{(h)}(y')\;=\;\Big(\frac{x}{y'}\Big)^M\,
{\bf Q}(y')\,{\bf Q}(x)^{-1}\;=\;
\Big(\frac{x'}{y'}\Big)^M\, {\bf K} \,{\bf Q}(y')\,{\bf Q}(x')^{-1} \;,
\end{equation}
and
\begin{equation}\label{nQhnorm1}
\overline{{\bf Q}}^{(h)}(q^{-1}y) \;=\; \Big(\frac{y}{x}\Big)^M
\,\varphi(y,y)^{-1}\,\overline{{\bf Q}}(q^{-1}y) \,
\Big(\lim_{y\to x}\,\varphi(y,y)^{-1} \,
\overline{{\bf Q}}(q^{-1}y) \Big)^{-1}\;.
\end{equation}

For completeness mention also the IRF version of the new 
operators \eqref{Qnew} 
\begin{equation}\label{2a}
\mathbb{Q}^{(h)}(y')_{a_N,\dots,a_1}^{b_N,\dots,b_1} \;=\; \prod_{i=1}^N
\frac{V_{y'/x'}(a_{i+1}-b_i) V_{x/y'}(b_i-a_i)}{V_{x/x'}(a_{i+1}-a_i)}\,,
\end{equation}
\begin{equation}\label{3}
\overline{\mathbb{Q}}^{(h)}(q^{-1}y)_{a_N,\dots,a_1}^{b_N,\dots,b_1} \;=\; \prod_{i=1}^N
\frac{V_{y/x}(b_{i}-a_i) V_{x/x'}(b_{i+1}-b_i)}{V_{y/x'}(b_{i+1}-b_i)}\,,
\end{equation}
which satisfy the factorization relation
\begin{equation}\label{4}
\mathbb{T}(y,y') \;=\; \mathbb{Q}^{(h)}(y')\, \overline{\mathbb{Q}}^{(h)}(q^{-1} y)\;,
\end{equation}
with the IRF transfer matrix \eqref{tmat} (where the field factors are set
as $\omega_h=\omega_v=1$). Note, that the operator \eqref{2a}
possesses quite remarkable locality properties. It is a rather sparse
matrix (similarly to the Hamiltonians \eqref{hamdef2}), which is very 
convenient for numerical calculations. 

\subsection{Finite-dimensional reduction }
Now consider the case when the parameter $2s$ in \eqref{homo} is 
an integer, $2s\in{\mathbb Z}_{\ge0}$. The corresponding Hamiltonians cannot be
obtained from \eqref{hamdef1}, \eqref{hamdef2} just by setting $2s$ to
such a value, because this operation does not commute with taking the
limit \eqref{T-expan}.  Actually, for $2s\in{\mathbb Z}_{\ge0}$ 
the limit \eqref{T-expan} exists only
when $\varepsilon=\varepsilon'$. So, there is only one Hamiltonian in
this case. The quantum space of states \eqref{qspace} becomes
reducible. Each factor of the product \eqref{qspace} splits into two
components
\begin{equation} 
{\mathbb C}^{(\infty)}={\mathbb C}^{(2s+1)}\oplus{\mathbb
  C}^{(\infty)}
\end{equation}
spanned by the vectors $|j\rangle$, $j\in {\cal
  I}_{s}=\{0,1,2,\ldots,2s\}$ and 
$j\in \overline{\cal I}_{s}=\{2s+1,2s+2,\ldots,\infty\}$, respectively,
which according to Sect.~\ref{red} 
support the $(2s+1)$-dimensional representations $\pi_s$ and
the infinite-finite dimensional representation $\pi^{+}_{-s-1}$ of the
$U_q(sl(2))$ algebra \eqref{uqsl2}. The Hamiltonian takes a block
triangular form, acting invariantly in two blocks where one takes
either the same $(2s+1)$-dimensional or the same infinite-finite dimensional
components in all factors of the product \eqref{qspace}.  

First, consider the $(2s+1)^N$ -dimensional block. The Hamiltonian is
given by the above formula \eqref{hamdef1}, but with the superscipt
$(\ell)$ omitted and 
\begin{equation}\label{H0}
\begin{array}{ll}
\displaystyle H_{i,i'}^{j,j'}\;=\; &
\delta_{i}^{j} \delta_{i'}^{j'}  (F_{i}(s) + F_{i'}(s) ) \\
[5mm]
& \displaystyle
+ 2\delta_{i+i',j+j'} \left(
\frac{\theta(i>j)}{[q^{i-j}]} \frac{V_{q^{-2s}}(j)}{V_{q^{-2s}}(i)}
+
\frac{\theta(i'>j')}{[q^{i'-j'}]} \frac{V_{q^{-2s}}(j')}{V_{q^{-2s}}(i')}
\right)\;.
\end{array}
\end{equation}
for $0\leq i+i'\leq 2s$ and 
\begin{equation}\label{H1}
\begin{array}{ll}
\displaystyle H_{i,i'}^{j,j'}\;=\; &
\delta_{i}^{j} \delta_{i'}^{j'}  (F_{2s-i}(s) + F_{2s-i'}(s) ) \\
[5mm]
& \displaystyle
+ 2\delta_{i+i',j+j'} \left(
\frac{\theta(j>i)}{[q^{j-i}]} \frac{V_{q^{-2s}}(j)}{V_{q^{-2s}}(i)}
+
\frac{\theta(j'>i')}{[q^{j'-i'}]} \frac{V_{q^{-2s}}(j')}{V_{q^{-2s}}(i')}
\right)\;.
\end{array}
\end{equation}
for $2s\leq i,j\leq 4s$. Note, that (\ref{H0}) and (\ref{H1}) display 
the reflection symmetry 
\begin{equation}
H_{i,i'}^{j,j'} \;=\; H_{2s-i,2s-i'}^{2s-j,2s-j'}\;.
\end{equation}

\bigskip

For the infinite dimensional block the local Hamiltonian is given by
the same formula \eqref{H1}, but with $i,i',j,j'> 2s$. 
It is worth noting that for $i>2s$ 
\begin{equation}\label{Fneg}
F_{2s-i}(s)\;=\;-\; 
\sum_{k=2s+1}^{i} \frac{q^k+q^{-k}}{q^k-q^{-k}} \;=\;F_{i-2s-1}(-s-1)\,-\,\frac{q^{2s+1}+q^{-2s-1}}{q^{2s+1}-q^{-2s-1}}\;.
\end{equation}
Moreover, with l'H\^opitale rule one obtains 
\begin{equation}\label{lH}
\frac{V_{q^{-2s}}(j)}{V_{q^{-2s}}(i)}\;=\;\frac{V_{q^{2s+2}}(i-2s-1)}{V_{q^{2s+2}}(j-2s-1)}\;.
\end{equation}
Taking into this into account it is easy to see that upon the shift of
all indices as $i\to i-2s-1$ the resulting Hamiltonian up to a
constant coincides with ${\bf H}^{(1)}+{\bf H}^{(2)}$ defined by
\eqref{hamdef1}, \eqref{hamdef2} with $s$ replaced by $-s-1$.  




\section{The Bethe Ansatz}
\subsection{Coordinate Bethe ansatz}
The commuting transfer matrices and ${\bf Q}$-operators 
\eqref{tcomm2}, \eqref{tq-comm} act in the infinite-dimensional quantum space 
${\mathcal H}^{\rm (vertex)}$ defined by \eqref{qspace}. As noted
before (see \eqref{down2}) 
all these matrices act invariantly in all finite dimensional
subspaces of ${\mathcal H}^{\rm (vertex)}$ with a fixed number of
particles $M=0,1,2,\ldots,\infty$. 
Their diagonalization problem can be solved via the coordinate Bethe ansatz 
by generalizing the results of
\cite{Bethe:1931,Lieb:1967zz,Sutherland:1967a,CPYang:1967,Sutherland:1967zz,
Baxter:1971cs}. Let $\widehat\sigma^+$ and $\widehat h$ 
be the ``particle creation'' and ``particle number''  
operators acting on the basis vectors
$|j\rangle \in{\mathbb C}^{\infty}$ as 
\begin{equation}\label{creation}
\widehat\sigma^+\,|j\rangle=|j+1\rangle\,,\qquad
\widehat h\,|j\rangle=j\,|j\rangle\,,
 \qquad j=0,1,2,\ldots,\infty\,.
\end{equation}
The eigenvectors describing $M$-particle states have the form 
\begin{equation}\label{cba-state}
{\bf \Psi}=\sum_{1\le r_1\le r_2 \le\cdots \le r_M \le N}
\Psi(r_1,r_2,\ldots ,r_M)\,\, \widehat\sigma_{r_M}^+\,\cdots\,
\,\widehat\sigma_{r_1}^+\,{\bf \Psi}_0\,,\qquad M\ge0\,,
\end{equation}
where the (unnormalized) Bethe ansatz wave function reads as
\begin{eqnarray}\label{wave}
\Psi(r_1,\ldots,r_M)=\sum_{\hat P}A_{\hat P}\ \prod_{m=1}^M \phi_{{\hat P}m}(r_m)
\end{eqnarray}
with suitably chosen coefficients $A_{\hat  P}$ and functions $\phi_m(r)$. 
The vector ${\bf \Psi}_0$ is the bare vacuum state 
\begin{equation} 
{\bf \Psi}_0=
\underbrace{|0\rangle\otimes|0\rangle\otimes\cdots\otimes|0\rangle}_N
\,,
\end{equation}
which
is the $M=0$ eigenstate of all operators in the commuting family, 
and $\widehat\sigma^+_r$ 
denotes the particle creation operator acting 
in the $r$-th factor of the tensor product \eqref{qspace}.
The summation is taken over all $M!$ permutations ${\hat  P}$ of the
integers $(1,2,\ldots,M)$. Remind that the 
number of particles $M$ is a conserved quantity, preserved by all
operators in the commuting family.
\medskip
The coefficients $A_{\hat P}$ read 
\begin{equation}\label{ap-def}
A_{\hat  P}=
\prod_{1\le j < m\le M}
\frac{q\,y^2_{\hat Pj\vphantom{{|}^{|}}}-q^{-1}\,y^2_{\hat Pm}}
{y^2_{\hat Pj}-y^2_{\hat Pm}}\ \,,
\end{equation}
where the set of complex numbers $\{y_m\}_{m=1}^M$ 
satisfies the system of algebraic equations (called the Bethe ansatz equations)
\begin{equation}\label{BA-eq}
\prod_{i=1}^N \frac{1-y_m^2/x_i^2{}_{\phantom{|}}\,}
{1-y_m^2/{x'}_i^{2}{}^{\phantom{|}}}
=-(\omega_h w)^{-1}\, \prod_{j=1}^M 
\frac{y_j^2-q^{+2}\,y_m^2}{y_j^2-q^{-2}\,y_m^2}\,,\qquad m=1,2,\ldots,M\,,
\end{equation}
and the functions $\phi_m(r)$ in \eqref{wave} have the form 
\begin{equation}\label{one-p}
\phi_m(r)=-[q]\,\left(\prod_{i=1}^{r-1} [y_m/x_i]\right)\ 
\left(\prod_{i=r+1}^{N} [y_m/x'_i]\right)\,.
\end{equation}
Note, that for $M=1$ the coefficient $A_{\hat P}\equiv1$.

The corresponding eigenvalues of 
the ${\bf Q}$-operators 
defined by \eqref{Q-def} 
\begin{equation} 
{\bf Q}(y)\,{\bf \Psi} =Q(y)\, {\bf
  \Psi}\,,\qquad \overline{{\bf Q}}(y)\,{\bf \Psi} =\overline{Q}(y)\,
{\bf \Psi}\,, 
\end{equation} 
(see also \eqref{Qelem},
\eqref{Qelem2} and \eqref{Qelem3} for explicit expressions for their
matrix elements) 
can be expressed via the Bethe roots
$\{y_m\}_{m=1}^M$, satisfying \eqref{BA-eq}. 
The eigenvalues of ${\bf Q}(y)$ are polynomials 
in the variable $y^2$,
\begin{equation} Q(y)=\prod_{m=1}^M
\Big(1-{y^2}/{y_m^2}\Big)\,,\label{Q-eig} 
\end{equation} 
while the eigenvalues of $\overline{\bf Q}(y)$ are, in general,
meromorphic functions of this variable.  They can be expressed through
$Q(y)$ by solving the quantum Wronskian relation
\eqref{qwrons},
\begin{equation}\label{Qbar-eig}
\overline{Q}(y)\;=\;Z_0^{-1}\,Q(y)\,
 \sum_{j=1}^\infty \big(\omega_h\,w\big)^{j-1} 
\frac{\varphi(q^jy,q^jy)}{{Q}(q^{j-1} y)\, {Q}(q^{j} y)}\,,\qquad |q|<1\,.
\end{equation}
Note, that separate terms in the above sum contain 
poles arising from the zeroes of $Q(q^j y)$,
$j=1,2,\ldots,\infty$. However, it is easy to see that all such 
poles cancel out between 
two successive terms of the sum by virtue of the Bethe ansatz
equations \eqref{BA-eq}. Thus, the poles in \eqref{Qbar-eig} could only  
originate from
the function $\varphi(y,y)$, defined in \eqref{phi-def}.
For arbitrary values of 
$\{s_n\}_{n=1}^N$ in \eqref{x-const} 
is has the poles located
at\footnote{%
Of course, these fixed ``kinematic'' 
poles could easily be removed by multiplying $\overline{Q}(y)$ by a
suitable factor. The functional relations should then
be modified accordingly.}  
\begin{equation}
y^2=q^{-2j} x_n^2\,,\qquad n=1,\ldots,N\,,\qquad j=1,2,\ldots,\infty\,.   
\end{equation}
However, if $2s_n\in{\mathbb
  Z}_{\ge0}$ for all $n=1,2,\ldots,N$ the function $\varphi(y,y)$
simplifies to 
\begin{equation}
\varphi(y,y)=\prod_{n=1}^N\, \prod_{j=1}^{2s_n}\,\Big(1-(y/q^{j}x_n)^2\Big)\,.
\end{equation}
and the eigenvalue $\overline{Q}(y)$ also becomes a polynomial in $y^2$.

The above construction for the eigenvectors \eqref{cba-state}-\eqref{one-p} 
is very similar to the Baxter's results \cite{Baxter:1971cs} 
for the eigenvectors of the inhomogeneous six-vertex model. The
difference is that the single particle wave function \eqref{one-p} and
the Bethe ansatz equations \eqref{BA-eq} in our case are slightly more
general than those for the 
spin $s=\tfrac{1}{2}$ six-vertex model where $x_n'=q\,x_n$  ($n=1,2,\ldots,N$).  
Moreover, the sum in \eqref{cba-state} allows any number of the 
particle positions $\{r_1,r_2,\ldots,r_N\}$ to coincide, whereas for
  the six-vertex model the inequalities in \eqref{cba-state} are
  strict, due to the ``exclusion principle'', that no more than one ``down
  arrow'' can be at the same site.   

\medskip\noindent
\subsection{Algebraic Bethe ansatz}
The Bethe state \eqref{cba-state}
may also be constructed within the framework of the QISM
\cite{Faddeev:1979gh,Takhtajan:1979iv}. Introduce the so-called
monodromy matrix 
\begin{equation}\label{M-in}
\boldsymbol{{\cal M}}(y)\equiv
\boldsymbol{{\cal M}}(y\,|\,x_N,x_N',\ldots,x_1,x'_1)=
\boldsymbol{{\cal L}}_{N}(y\,|\,x_N,x_N')\,
\boldsymbol{{\cal L}}_{N-1}(y\,|\,x_{N-1},x_{N-1}')
\cdots\boldsymbol{{\cal L}}_{1}(y\,|\,x_{1},x_{1}')\,,
\end{equation}
where $\boldsymbol{{\cal L}}_{n}(y\,|\,x_n,x'_n)$ 
is the two-by-two matrix 
\eqref{Lop}, whose entries are elements of the $U_q(sl(2))$
algebra \eqref{uqsl2}, acting in the $n$-th
component of the tensor product ${\cal H}^{\rm(vertex)}$,  
defined by \eqref{qspace}.
Each these components realizes an infinite-dimensional highest weight
representation of this algebra given by \eqref{piplus_s1} 
(where $s_1$ for the $n$-th component is 
replaced by $s_n$, such that $q^{2s_n}=x'_n/x_n$).
It is convenient to denote the entries of the monodromy matrix as
\begin{equation}\label{Mmat}
\bm{{\cal M}}(y)=
\left(\begin{array}{cc} 
{ \hat{ {\mathsf A}}}(y) & { \hat{ {\mathsf B}}} (y)\\[0.2cm] 
{ \hat{ {\mathsf C}}} (y) & { \hat{ {\mathsf D}}}(y)
\end{array}\right)\,,
\end{equation}
where $ {\hat{ {\mathsf A}}},\, {\hat{ {\mathsf B}}},\, 
{ \hat{ {\mathsf C}}},\,{ \hat{ {\mathsf D}}}$
are operators acting in the quantum space \eqref{qspace}.
With the above notations, the transfer matrix ${\bf T}_{\tfrac{1}{2}}(y)$,
defined by \eqref{ts-def} with $s=\tfrac{1}{2}$, 
is given by
\begin{equation}\label{Tmat1}
{\bf T}_{\tfrac{1}{2}}(y)=\Big(q^M \prod_{n=1}^N (y/x'_n)\Big)\,
\Big({ \hat{ {\mathsf A}}}(y) +\omega_h\, { \hat{ {\mathsf D}}}(y)\Big)\  ,
\end{equation}
while
the Bethe state \eqref{cba-state}-\eqref{one-p},
with the set $\{y_m\}_{m=1}^M$ solving the Bethe ansatz equations, is expressed as 
\begin{equation}\label{aba-state}
\boldsymbol{\Psi}={ \hat{ {\mathsf B}}}(y_M)\,\cdots\,
{ \hat{ {\mathsf B}}}(y_2)\, { \hat{ {\mathsf B}}}(y_1)\,{\bf {\Psi}}_0\,.
\end{equation}
The corresponding eigenvalue of the transfer matrix \eqref{Tmat1},
\begin{equation}
{\bf T}_{\tfrac{1}{2}}(y)\,{\bf \Psi}={T}_{\tfrac{1}{2}}(y)\,{\bf \Psi}\,,
\end{equation}
reads \cite{Takhtajan:1979iv,Korepin:1982gg} 
\begin{equation}
T_{\tfrac{1}{2}}(y)\;=\;f(y) \prod_{j=1}^M \frac{1-q^2 y^2/y_j^2}{1-y^2/y_j^2} \;+\;
(\omega_h w)\, g(y) \prod_{j=1}^M \frac{1-q^{-2} y^2/y_j^2}{1-y^2/y_j^2} \;.
\end{equation}
where $f(y)$ and $g(y)$ are defined in \eqref{fg-def}.
Combining the last formula with
\eqref{Q-norm} and \eqref{tqeq} one immediately deduces
the expression \eqref{Q-eig} for the eigenvalue of the operator ${\bf Q}(y)$. 

Finally, from the definitions \eqref{M-in}, \eqref{Mmat} it follows that 
\begin{equation}
\hat{{\mathsf B}}(y)=-\sum_{r=1}^N
\left(\prod_{i=1}^{r-1} [q^{-\widehat{h}_i}\,y_m/x_i]\right)
\,[q^{\widehat{h}_r}]\,\widehat
\sigma_r^+\,
\left(\prod_{i=r+1}^{N} [q^{\widehat{h}_i}\,y_m/x'_i]\right)\,,
\end{equation}
where $\widehat h_i$ is the particle number operator \eqref{creation}
acting at the
$i$-th component of the product \eqref{qspace}.
Substituting this into 
\eqref{aba-state} and rearranging similar terms in a special way 
one obtains the coordinate
Bethe ansatz expression for the eigenvectors given by
\eqref{cba-state}-\eqref{ap-def}. 
In the context of the XXX related models 
this statement first appeared in \cite{Izergin:1987} 
(see eq.(4.8) therein). For generalizations to the XXZ case 
and derivations see \cite{Slavnov:2007,Fuksa:2017}.
Note that this equivalence 
between the two forms of the ``ansatz'' for the eigenvectors holds 
for arbitrary values of $\{y_m\}_{m=1}^M$, in particular, it does not   
require them to satisfy the Bethe ansatz equations \eqref{BA-eq}.






\section{Generalizations to other models}


Apparently, all vertex models 
associated with quantized affine Lie algebras \cite{Bazhanov:1984gu,Jimbo:1986}
and superalgebras \cite{Bazhanov:1986av} can
be reformulated as Ising-type models. 
Indeed, the most important algebraic structure required for such a
reformulation --- the 3-dimensional interpretation based on the
tetrahedron equation --- is known for a majority of these models, 
\cite{Bazhanov:2005as,Sergeev:2009,
Bosnjak:2016oze,KunS1,KunS2,KunS3}.  
In support of our statement we present here a new class of  
Ising-type models related to the $U_q(\widehat{sl}(n))$ algebra, 
extending the results of this 
paper which was devoted to the $n=2$ case. 
The formulation of these new
models follows
exactly the same description as in Sect.~\ref{isingtype}, except
that each site of the lattice now carries an $(n-1)$-component integer
spin 
variable $\boldsymbol{a}=(a_1,a_2,\ldots , a_{n-1})\in{\mathbb
  Z}^{(n-1)}$ and the Boltzmann weight function $V$, entering the definition 
of the edge weights \eqref{twotypes}, is now replaced with 
\begin{equation}\label{Veight}
V_x(\boldsymbol{a},\boldsymbol{b}) \;=\; q^{2Q(\boldsymbol{a},\boldsymbol{b})-Q(\boldsymbol{a},\boldsymbol{a})-Q(\boldsymbol{b},\boldsymbol{b})}\;
\left(\frac{q}{x}\right)^{a_0-b_0}\frac{(x^2;q^2)_{a_0-b_0}}{\displaystyle\prod_{i=1}^{n-1}(q^2;q^2)_{a_i-b_i}}\;.
\end{equation}
Note, that this function 
depends on both spin variables $\boldsymbol{a}$ and 
$\boldsymbol{b}$ at the ends of the edge, rather than their
difference, and 
\begin{equation}
a_0=\sum_{i=1}^{n-1}a_i\,,\qquad
\qquad Q(\boldsymbol{a},\boldsymbol{b}) \;=\; \sum_{1\le i<j\le n-1}
a_ib_j\,.
\qquad 
\end{equation}
Similarly, to \eqref{W-def1}, \eqref{W-def2} one can define 
the corresponding IRF weight 
\begin{subequations}\label{SS2}
\begin{eqnarray}
\label{W-def1-2}
\displaystyle \mathcal{W}^{sl(n)}(\boldsymbol{a},\boldsymbol{b},\boldsymbol{c},
\boldsymbol{d}\,|\,x,x',y,x')&=&
\frac{\displaystyle V_{y/y'}(\boldsymbol{b},\boldsymbol{d})}
{\displaystyle V_{y/y'}(\boldsymbol{a},\boldsymbol{c})}\;
\sum_{\boldsymbol{n}=\max(\boldsymbol{b},\boldsymbol{c})}^{\boldsymbol{a}}
\frac{\displaystyle 
V_{x/y'}(\boldsymbol{a},\boldsymbol{n})
V_{y'/x'}(\boldsymbol{n},\boldsymbol{b})
V_{y/x}(\boldsymbol{n},\boldsymbol{c})}
{\displaystyle V_{y/x'}(\boldsymbol{n},\boldsymbol{d})}\,,\\[.3cm]
\label{W-def2-2}
&=&
\frac{\displaystyle V_{x/x'}(\boldsymbol{a},\boldsymbol{b})}{\displaystyle V_{x/x'}(\boldsymbol{c},\boldsymbol{d})}\;
\sum_{\boldsymbol{n}=\boldsymbol{d}}^{\min(\boldsymbol{b},\boldsymbol{c})}
\frac{\displaystyle 
V_{x/y'}(\boldsymbol{n},\boldsymbol{d})
V_{y'/x'}(\boldsymbol{c},\boldsymbol{n})
V_{y/x}(\boldsymbol{b},\boldsymbol{n})}
{\displaystyle V_{y/x'}(\boldsymbol{a},\boldsymbol{n})}\;.
\end{eqnarray}
\end{subequations}
where the summation is taken over $(n-1)$ independent integer variables $n_i$, 
such that $\max(b_i,c_i)\leq n_i \leq a_i$ in (\ref{W-def1-2}) and,
similarly, in \eqref{W-def2-2}. The corner spins 
should obey the relations 
$a_i\geq b_i,c_i\geq d_i$  $(i=1,2,\ldots,n-1)$, otherwise the weights  
$\mathcal{W}^{sl(n)}(\boldsymbol{a},\boldsymbol{b},\boldsymbol{c},
\boldsymbol{d})$ 
are assumed to vanish identically.
The equality of two different
expressions in \eqref{SS2} is the star-star relation, which ensure the
integrability of the model. It and has exactly the same graphical representation as in Fig.~\ref{sspic} and reduces to \eqref{ssrel} for $n=2$. For $n>3$ 
the star-star relation \eqref{SS2} is a new identity which we
currently claim as a conjecture (though we have thoroughly verified 
it for many particular cases).  

It is not difficult to check that the IRF weights
$\mathcal{W}^{sl(n)}(\boldsymbol{a},\boldsymbol{b},
\boldsymbol{c},\boldsymbol{d}\,|\,x,x',y,x')$
only depend on the differences between the corner spins, so one can
trivially convert these weights into a vertex $R$-matrix (similarly to
\eqref{R-def}). The resulting solutions of
the Yang-Baxter equation are associated with infinite-dimensional
highest weight evaluations representations of $U_q(\widehat{sl}(n))$
algebra, namely, the symmetric tensor representations
\cite{Bosnjak:2016oze}.      

The vertex $R$-matrix, corresponding to \eqref{SS2} reads
\begin{equation}
{\mathcal
  R}^{sl(n)}(x,x',y,y')^{\bm{j}_1,\bm{j}_2}_{\bm{i}_1,\bm{i}_2}=\delta_{\bm{{i}}_1+\bm{i}_2,
\bm{j}_1+\bm{j}_2}\, 
{{\mathcal W}}^{sl(n)}(\bm{i}_1+\bm{i}_2,\bm{j}_2,\bm{i}_1,0\,|x,x',y,y')\,.\label{Rsln-def}
\end{equation}
where $\bm{j}_1,\bm{j}_2,\bm{i}_1,\bm{i}_2\in {\mathbb
  Z}_{\ge0}^{(n-1)}$.
Note, if $x'=q^{\mu_1}x$ and $y'=q^{\mu_2}y$, with $\mu_1,\mu_2\in{\mathbb
  Z}_{\ge0}$, then the above $R$-matrix admits a finite-dimensional
reduction. For instance, 
the Cherednik-Perk-Shultz $R$-matrix \cite{Cherednik:1980ey,Perk:1981}
(vector representations of $U_q(sl(n))$) 
appears in the case $x'=qx$, $y'=qy$.
Let
\begin{equation}
\boldsymbol{e}_0\;=\;(0,\ldots,0)\;,\quad \boldsymbol{e}_i\;=\;(0,\ldots,\underbrace{1}_{i^{th} \;\textrm{place}},\ldots,0)\;,
\end{equation}
and $\varepsilon_{i,j}=-\varepsilon_{j,i}=1$ for $i<j$, except the
case $i=0$ or $j=0$ where $\varepsilon_{i,j}\;=\;0$. 
Then, 
\begin{equation}
\mathcal{R}_{\boldsymbol{e}_i,\boldsymbol{e}_i}^{\boldsymbol{e}_i,\boldsymbol{e}_i}=1\;,\qquad
\mathcal{R}_{\boldsymbol{e}_i,\boldsymbol{e}_j}^{\boldsymbol{e}_i,\boldsymbol{e}_j}\;=\;
q^{\varepsilon_{i,j}} \frac{[x/y]}{[qx/y]}\;,\qquad
\mathcal{R}_{\boldsymbol{e}_i,\boldsymbol{e}_j}^{\boldsymbol{e}_j,\boldsymbol{e}_i}\;=\;
\left(\frac{y}{x}\right)^{\varepsilon_{i,j}}  \frac{[q]}{[qx/y]}\;,\qquad
i\not=j\,,
\end{equation}
where we have omitted the superscript ``$sl(n)$''.
      
\section{Conclusion}
In this paper we have presented a new approach to the six-vertex model and
its higher spin generalizations. We have reformulated these models as
Ising-type models with only a two-spin interaction across edges of the
square lattice. This reformulation leads to a significant simplification 
of the algebraic theory of the higher spin six-vertex model, since 
the edge Boltzmann weights are given by a very simple formula
\eqref{Vdef}. As a result explicit expressions for associated 
$R$-matrices, as well as their properties, become much more transparent. 


An Ising-type model can be equivalently converted into a
vertex model in (at least) two different ways. In the first way, 
the $R$-matrix defining the vertex weights is deduced from the
Bolzmann weights of the four-edge star (like the one shown in
Fig.~\ref{fig-star}) by using the star-square transformation followed
by the IRF-vertex correpondence. In our case
this procedure leads to Eq.\eqref{R-def2}. 
In the second way the vertex-type $R$-matrix is identified with ``box
diagram'', shown in Fig.~\ref{fig-box}, that leads to Eq.\eqref{def3}.
In Sect.~\ref{descen} we have shown these two $R$-matrices are related
by a linear relation \eqref{fourier}, which is not unnatural to expect,
since they both describe exactly the same model. 

From the quantum group point of view these $R$-matrices 
intertwine equivalent 
representations of the $U_q(\widehat{sl}(2))$ algebra, connected 
by a similarity transformation in the representation space. For the first
representation \eqref{R-def2} the Cartan element is
diagonal, while for \eqref{Lop-new} it is realised as a shift operator. 
It should be stressed, that 
it is the shift operator realization of the Cartan elements
that leads to the factorized $R$-matrices. This fact has been previously
observed \cite{Bazhanov:1989nc} for the chiral Potts model and its
generalizations \cite{BKMS}, which are related to the
$U_q(\widehat{sl}(n))$ algebra. We expect a similar phenomenon to take place
for all other quantized affine algebras.  
One disadvantage of the factorized $R$-matrices 
is the necessity to work with
infinite-dimensional representations, whereas the models with finite
number of discrete spin states only appear via a reduction. In
practice, however, this
does not create real problems. Actually, it
would be fair to say that certain infinite-dimensional representations (in
particular, the $q$-oscillator representations) 
cannot be avoided, since they 
are required for the transfer matrix
construction of the Baxter ${\bf Q}$-operators 
\cite{Bazhanov:1996dr,Bazhanov:1998dq}. 
Moreover, they play an important role of ``prefundamental'' 
representaions in the theory of quantized affine algebras
\cite{Hernandez:2012,Frenkel:2012}. 


As is well known, the six-vertex model has a vast number of important
applications to various problems of physics and mathematics (for reviews
see \cite{Lieb:1980ix,Reshetikhin:2010}). 
Mention, in particular, the enumeration of
the alternating sign matrices \cite{RazStr,BGN}, asymmetric stochastic
exclusion processes 
and stochastic lattice models \cite{Borodin:2014,Garbali:2016}.    
Most recently the (inhomogeneous) six-vertex model was involved in
fascinating connections to black hole sigma models in quantum field
theory %
\cite{Jacobsen:2005xz,Ikhlef:2008zz,Ikhlef:2011ay,Frahm:2012eb,Candu:2013fva,Frahm:2013cma,Bazhanov:2020,Bazhanov:2020uju}.
It would be interesting to reconsider these problems and
connections in view of the Ising-type structures studied in our work.  

\section{Acknowledgements}
The authors thank S.L.Lukyanov, V.V.Mangazeev and N.Y.Reshetikhin for 
stimulating discussions. One us (VVB) thanks the Sydney Mathematical
Research Institute and especially Prof. A.Molev 
for hospitality at Sydney when some parts of this work were completed.

\medskip
\noindent
SMS acknowledges the support of the Australian Research Council grant 
DP190103144. 

\app{The star-star and Yang-Baxter relations}\label{appA}
\addcontentsline{toc}{section}{A. The star-star and Yang-Baxter relations}
\sapp{Proof of the star-star relation.}
In this Appendix we will prove the star-star \eqref{ssrel} and 
the Yang-Baxter \eqref{ybe}
relations. 
For readers' convenience, let us reproduce here the definitions of the
IRF-type weights given by \eqref{W-def} and \eqref{W-def2}, 
\begin{equation}
{\mathcal W}(a,b,c,d\,|\,x,x',y,y')=
\frac{V_{y/y'}(b-d)}{V_{y/y'}(a-c)}\,
\ 
\sum_{n=\max(b,c)}^{a}\frac{V_{x/y'}(a-n)V_{y'/x'}(n-b)
V_{y/x}(n-c)}{V_{y/x'}(n-d)}\,.\label{star1}
\end{equation}
\begin{equation}\label{star2}
{\overline {\mathcal W}}(a,b,c,d\,|\,x,x',y,y')
=\frac{V_{x/x'}(a-b)}{V_{x/x'}(c-d)}
\,\ \sum_{{n}=a}^{\min(b,c)} \ \ \frac{V_{x/y'}({n}-d)
V_{y'/x'}(c-{n}) 
V_{y/x}(b-{n})}
{V_{y/x'}(a-{n})}\,.
\end{equation}
We need to prove that the two definitions lead to the same result
\begin{equation}\label{ssrel2}
{\mathcal W}(a,b,c,d\,|\,x,x',y,y')=\overline{{\mathcal W}}(a,b,c,d\,|\,x,x',y,y')\,.
\end{equation}
This is the star-star relation, which states the equality of the Boltzmann
weights of the two ``stars'' shown in Fig.~\ref{sspic}. 
Below it will be
more convenient to modify the graphical representations
\eqref{twotypes} for the edge weights. We will abolish the
rapidity lines, but instead indicate the rapidity ratio arguments of the
weights for all edges of the lattice. These arguments are shown
near the edges, as in the pictures below,
\begin{equation}\label{gnote2}
\input fig2nc.tex
\end{equation}
Otherwise, the rules remain unchanged: the edge arrow pointing from
$b$ to $a$ means that the spin difference argument of the weight
function is 
$(a-b)$ rather than $(b-a)$, 
moreover, the colour distiguishes the type $(i)$
(blue) and the type $(ii)$ (red) edges with the Boltzmann weights
$V_x(a-b)$ and $(V_x(a-b))^{-1}$, respectively. 
Remind, that the
function $V_x(n)$ is defined in \eqref{Vdef}. 
With the new 
convensions the star-star relation \eqref{ssrel2} is presented in
Fig.~\ref{fig-star2}. 
\begin{figure}[ht]
\begin{center}
\input stst2.tex
\caption{The star-star relation \eqref{ssrel2} 
with alternative graphical notations
  for the edge weight, as in \eqref{gnote2}.
\label{fig-star2}}
\end{center}
\end{figure}




Below we will use the second Sears's transformation formula for the basic
hypergeometric series (see \cite{Sears:1951}), 
\begin{equation}\label{S2}
\phantom{|}_4\phi_3\left(\begin{array}{cccc|}
q^{-2n} & {\mathsf a} & {\mathsf b} & {\mathsf c}\\ & {\mathsf d} & {\mathsf e} & {\mathsf f}\end{array}\,\,q^2,q^2\right) \;=\;
\frac{\displaystyle \left({\mathsf a},\,\frac{{\mathsf e}\, {\mathsf f}}{{\mathsf a}\, {\mathsf b}},\,\frac{{\mathsf e}\, {\mathsf f}}{{\mathsf a}\,
    {\mathsf c}};\,q^2\right)_n}{\displaystyle \left({\mathsf e},\,{\mathsf f},\,\frac{{\mathsf e}\, {\mathsf f}}{{\mathsf a}\, {\mathsf b}\, {\mathsf c}};\,q^2\right)_n}\ \phantom{|}_4\phi_3\left(\begin{array}{cccc|}
\displaystyle q^{-2n} & \displaystyle \frac{{\mathsf e}}{{\mathsf a}} & \displaystyle\frac{{\mathsf f}}{{\mathsf a}} & \displaystyle \frac{{\mathsf e}\,
  {\mathsf f}}{{\mathsf a}\, {\mathsf b}\, {\mathsf c}} 
\\&&&  \\
\displaystyle & \displaystyle \frac{{\mathsf e}\, {\mathsf f}}{{\mathsf a}\, {\mathsf b}} & \displaystyle \frac{{\mathsf e}\, {\mathsf f}}{{\mathsf a}\, {\mathsf c}}  & \displaystyle
\frac{q^{2-2n}}{{\mathsf a}}\end{array} \,\,q^2,q^2\right)\;,
\end{equation}
where the integer $n\ge0$ and the other
parameters satisfiy the constraint
\begin{equation}
{\mathsf a}\, {\mathsf b} \,{\mathsf c} \;=\;q^{2n-2}\, {\mathsf d}\, {\mathsf e}\, {\mathsf f}\;.
\end{equation}
Here we use the standard notations 
\begin{equation}\label{phi43}
\phantom{|}_{r+1}\phi_r\left(\begin{array}{ccccc|}
q^{-2n}, & {\mathsf a}_1, &{\mathsf a}_2, & \ldots\ , & {\mathsf a}_r\\ & {\mathsf b}_1, &{\mathsf b}_2 & \ldots\ , &
{\mathsf b}_r\end{array}\,\,q^2,q^2\right) 
\;=\;
\sum_{k=0}^n \frac{(q^{-2n},{\mathsf a}_1, {\mathsf a}_2, \ldots,
  {\mathsf a}_r;q^2)_k}{(q^2,{\mathsf b}_1,{\mathsf b}_2, 
  \ldots, {\mathsf b}_r;q^2)_k}\, q^{2k}\;,
\end{equation}
and 
\begin{equation}\label{qprod}
({\mathsf a};q^2)_k=\prod_{\ell=0}^{k-1}(1-{\mathsf a}\, q^{2\ell})\,;\qquad 
({\mathsf a}_1,{\mathsf a}_2,\ldots,{\mathsf a}_m;q^2)_k=({\mathsf a}_1;q^2)_k\,({\mathsf a}_2;q^2)_k\,\cdots\,
({\mathsf a}_m;q^2)_k\,.
\end{equation}

We are now ready to prove the star-star relation
\eqref{ssrel2}. First, consider the case 
$b\ge c$.
From the definition \eqref{Vdef} it is easy to verify, that
\begin{equation}\label{cor}
V_x(a-n) \;=\; V_x(a) \left(\frac{q}{x}\right)^n
\frac{(q^{-2a};q^2)_n}
{(q^{2-2a}\,x^{-2};q^2)_n}\;.
\end{equation}
Using this formula and the definition \eqref{phi43} 
one can rewrite \eqref{star1} and \eqref{star2} in the form
\begin{equation}\label{SL}
\begin{array}{l}
\displaystyle {\mathcal W}\;=\;\frac{V_{y/y'}(b-d)}{V_{y/y'}(a-c)} 
\frac{V_{y'/x'}(a-b) V_{y/x}(a-c)}{V_{y/x'}(a-d)}\\
[7mm]
\displaystyle \times
\phantom{|}_4\phi_3\left(\begin{array}{rrrr|}
\displaystyle q^{-2(a-b)}; & \displaystyle q^{-2(a-c)}, & \displaystyle (x/y')^2, & \displaystyle q^{2(1-a+d)} 
(x'/y)^{2}\\
[5mm]
\displaystyle  & \displaystyle q^{2(1-a+b)} (x'/y')^{2}, & \displaystyle q^{-2(a-d)}, & 
\displaystyle q^{2(1-a+c)} (x/y)^{2} \end{array}\,\,q^2,q^2\right) \;,
\end{array}
\end{equation}
and
\begin{equation}\label{SR}
\begin{array}{l}
\displaystyle \overline{{\mathcal W}}\;=\;\frac{V_{x/x'}(a-b)}{V_{x/x'}(c-d)} \frac{V_{y'/x'}(c-d) 
V_{y/x}(b-d)}{V_{y/x'}(a-d)}\\
[7mm]
\displaystyle \times
\phantom{|}_4\phi_3\left(\begin{array}{rrrr|}
\displaystyle q^{-2(c-d)}; & \displaystyle q^{-2(b-d)}, & \displaystyle (x/y')^2, & \displaystyle q^{2(1-a+d)} 
(x'/y)^{2}\\
[5mm]
\displaystyle  & \displaystyle q^{2(1-c+d)} (x'/y')^{2}, & \displaystyle q^{-2(a-d)}, & \displaystyle
q^{2(1-b+d)} 
(x/y)^{2} \end{array}\,\,q^2,q^2\right) \;.
\end{array}
\end{equation}
where we have omitted the arguments of ${\mathcal W}$ and $\overline{{\mathcal W}}$ since
they do not change during these calculations.
Applying the second Sears transformation (\ref{S2}) to both (\ref{SL}) and (\ref{SR}), one obtains 
\begin{equation}
\begin{array}{l}
\displaystyle {\mathcal W}\;=
\;\frac{V_{x/x'}(a-b)\, V_{y/y'}(a-d)\, V_{y/x}(b-c) }
{ V_{y/y'}(a-c)\, V_{y/x'}(a-d)} \\
[7mm]
\displaystyle \times
\phantom{|}_4\phi_3\left(\begin{array}{rrrr|}
\displaystyle q^{-2(a-b)}; & \displaystyle q^{-2(c-d)}, & \displaystyle q^2\,(x/y)^{2}, & \displaystyle (y'/x')^2\\
[5mm]
\displaystyle  & \displaystyle q^{2(1-a+d)}\, (y'/y)^{2}, & \displaystyle (x/x')^2, & \displaystyle
q^{2(1-c+b)} 
\end{array}\,\,q^2,q^2\right) \;,
\end{array}\label{SL2}
\end{equation}
and
\begin{equation}
\begin{array}{l}
\displaystyle \overline{{\mathcal W}}\;=\;\frac{V_{x/x'}(a-b)\, V_{y/y'}(a-d) \,
V_{y/x}(b-c) }{ V_{y/y'}(a-c)\, V_{y/x'}(a-d)} \\
[7mm]
\displaystyle \times
\phantom{|}_4\phi_3\left(\begin{array}{rrrr|}
\displaystyle q^{-2(c-d)}; & \displaystyle q^{-2(a-b)}, & \displaystyle q^2\,(x/y)^{2}, & \displaystyle (y'/x')^2\\
[5mm]
\displaystyle  & \displaystyle q^{2(1-a+d)}\, (y'/y)^{2}, & \displaystyle (x/x')^2, & \displaystyle q^{2(1-c+b)} \end{array}\,\,q^2,q^2\right) \;,
\end{array}\label{SR2}
\end{equation}



The resulting expressions coincide due to an obvious symmetry of
${}_4\varphi_3$ upon permutations of its upper line
arguments. Evidently, this proves the star-star relation \eqref{ssrel2}.
The above computation is valid for $b\geq c$. 
The case $b<c$ can be obtained by the transformation  $b\leftrightarrow c$,
$x'\leftrightarrow y^{-1}$ and $y'\leftrightarrow x^{-1}$, which maps
the star-star relation into itself. 


\sapp{Proof of the IRF-type Yang-Baxter equation\label{app-irf-ybe}}
Let us now prove the Yang-Baxter relation \eqref{ybe}. 
The star-star
relation \eqref{ssrel2} could be used in two ways. Firstly, one could
express the sum in \eqref{star1} via the sum in \eqref{star2} (of
course, taking into account the pre-factors in front of the sums). 
We will call this
operation as the ``{\bf direct star-star map centered at $n$}'', where $n$
denotes the summation index in the first sum \eqref{star1}, which is
the central spin in the left star depicted in Fig.~\ref{fig-star2}. 
Similarly, one could
express the second sum \eqref{star2} via the sum in \eqref{star1}.  
We will call this
operation as the ``{\bf reverse star-star map centered at ${n}$}'',
where ${n}$ now refers to the central spin in the right star 
in Fig.~\ref{fig-star2}.

Return now to the Yang-Baxter equation \eqref{ybe}. For convenience we
reproduce it here
\begin{equation}\label{ybe2}
\begin{array}{l}
\displaystyle \sum_{a\in{\mathbb Z}} 
{\mathcal W}(g,a,b,f\,|\,{\boldsymbol x},{\boldsymbol y})\ 
{\mathcal W}(c,e,g,a\,|\,{\boldsymbol x},{\boldsymbol z})\ 
{\mathcal W}(e,d,a,f\,|\,{\boldsymbol y},{\boldsymbol z})\;=\;\\
[3mm]
\displaystyle \qquad\qquad\qquad\qquad\qquad=\;
\sum_{h\in{\mathbb Z}} 
{\mathcal W}(c,h,g,b\,|\,{\boldsymbol y},{\boldsymbol z})\ 
{\mathcal W}(h , d ,b , f\,|\,{\boldsymbol x},{\boldsymbol z})\ 
{\mathcal W}(c , e ,h , d\,|\,{\boldsymbol x},{\boldsymbol y})
\end{array}
\end{equation}
where the bold symbols denote the pairs of spectral variables 
${\boldsymbol x}=(x,x')$, ${\boldsymbol y}=(y,y')$ and ${\boldsymbol z}=(z,z')$. Consider the LHS of this
equation and use the expression \eqref{star1} for each of the weights 
${\mathcal W}(g,a,b,f\,|\,{\boldsymbol x},{\boldsymbol y})$,
${\mathcal W}(c,e,g,a\,|\,{\boldsymbol x},{\boldsymbol z})$ and 
${\mathcal W}(e,d,a,f\,|\,{\boldsymbol y},{\boldsymbol z})$
denoting the central spins by $n_1$, $n_2$ and $n_3$,
respectively. Then, with the graphical representation of the weight
\eqref{star1} 
shown on the left side of Fig.~\ref{fig-star2}, the complete LHS of
\eqref{ybe2} will be represented as shown below
\input ybepic1.tex
Note, that the external sites, corresponding to fixed spins in
\eqref{ybe2}, are represented by open circles, while the internal
sites, corresponding to the summation spin indices, are shown by
filled circles. 

Now, make the reverse star map centered at $a$. This transformation
involves six edges, the four edges surrounding the site $a$ and two boundary
edges $(f, n_1)$ and $(f, n_3)$ (see the shaded
area on the figure \eqref{ybepic1} above). The result is shown in the
left side of \eqref{ybepic2}. 
\input ybepic2.tex
\noindent
The shading on this figure now 
indicate the area of the next transformation, which
is the reverse star map centered at $n_3$. The result is 
shown on the right side of \eqref{ybepic2}. Next, make the 
reverse star map centered at $n_1$ followed by the direct star map 
centred at $n_2$. The results are shown on the left and right sides of 
\eqref{ybepic3}, respectively. 
\input ybepic3.tex
Finally, redenoting the summation indices
\begin{equation}
n_3\to n_1,\qquad a\to n_2,\qquad n_1\to n_3,\qquad  n_2\to h\,,
\end{equation}
one converts the LHS of \eqref{ybe2} to its RHS. This proves the
Yang-Baxter relation \eqref{ybe} and commutativity of transfer
matrices \eqref{tcomm} in the main text.
\sapp{Proof of the vertex-type Yang-Baxter equation\label{box-YBE}}
Now consider the Yang-Baxter equation \eqref{ybe4} 
for the $R$-matrix \eqref{def3}. For convenience let us reproduce it here,
\begin{equation}\label{ybe4a}
\sum_{a',\,b',\,c'}{\mathbb R}_{a\,,\,b}^{a',\,b'}({\boldsymbol x},{\boldsymbol y}) \,
{\mathbb R}_{a'\;,\,c}^{a'',\,c'}({\boldsymbol x},{\boldsymbol z})\,
{\mathbb R}_{b'\;,\,c'}^{b'',\,c''}({\boldsymbol y},{\boldsymbol z})
\;=\;
\sum_{a',\,b',\,c'}
{\mathbb R}_{b\,,\,c}^{b',\,c'}({\boldsymbol y},{\boldsymbol z})\,
{\mathbb R}_{a\,,\,c'}^{a',\,c''}({\boldsymbol x},{\boldsymbol z})\,
{\mathbb R}_{a'\;,\,b'}^{a'',\,b''}({\boldsymbol x},{\boldsymbol y})\,.
\end{equation}
With the graphical notations \eqref{gnote2} the LHS of this equation can be represented as
\input ybepic4.tex  
\noindent
The proof of this equation consists of the repeated application of
the star-star relation \eqref{ssrel2}, represented graphically in
\eqref{gnote2}. The working is very similar to the proof of
\eqref{ybe2} in previous subsection. In particular we will use the
notions of the ``direct'' and ``reverse'' star-star maps introduced in
the first paragraph of Sect.~\ref{app-irf-ybe}.

Start from the left hand side of (\ref{ybe4a}) and then 
\begin{enumerate}
\item make the reverse star map centred at $a'$,
\item make the reverse star map centred at $b'$,
\item make the reverse star map centred at $c'$,
\item make the direct star map centred at $a'$ again,
\end{enumerate}
and then making the change
\begin{equation}
a'\to a'\;,\quad b'\to c'\;,\quad c'\to b'\;,
\end{equation}
one converts the left hand side of (\ref{ybe4a}) to its right hand side.



\app{Star-triangle relation as a primary integrability condition\label{appB}}
\addcontentsline{toc}{section}{B. Star-triangle relation as a primary
  integrability condition} 
\sapp{Star-triangle relation.}
Here we prove the star-triangle relation \eqref{str} by reducing it  
to the Pfaff-Saalshch\"utz-Jackson summation formula
(see eq.(3.5.1) in \cite{Gasper:1995}),
\begin{equation}\label{pfaff}
\phantom{|}_3\phi_2\left(\begin{array}{ccc|} \displaystyle q^{-2n}, & \displaystyle {\mathsf a}, &
  \displaystyle {\mathsf b}\\ \displaystyle & \displaystyle {\mathsf c}, & \displaystyle \ q^{2-2n}\,{\mathsf a}\,{\mathsf b} /{{\mathsf c}}\end{array} \,\,
q^2, q^2\right) 
\;=\; 
\frac{\displaystyle \left({{\mathsf c}}/{{\mathsf a}},\,{{\mathsf c}}/{{\mathsf b}}\,;\,q^2\right)_n}
{\displaystyle \left(\,{\mathsf c},\ {{\mathsf c}}/{{\mathsf a}\,{\mathsf b}}\,;\,q^2\right)_n}
\end{equation} 
where the notations are defined in \eqref{phi43} and \eqref{qprod}.
Here $n$ is a non-negative integer $n\ge 0$, 
and $q$, ${\mathsf a}$, ${\mathsf b}$ and ${\mathsf c}$ stand for
arbitrary complex parameters. 

For readers' convenience we reproduce the relations \eqref{str} here
\begin{subequations}
\label{stra}
\begin{equation}
\sum_d \,V_{y/z}(d-a)\, W_{z,x}(b-d)\, V_{x/y}(c-d)  \;=\;
W_{z,y}(b-c) \,V_{x/z}(c-a)  \,W_{y,x}(b-a)\,,\label{str1a}
\end{equation}
and
\begin{equation}
\sum_d V_{y/z}(a-d)\,W_{z,x}(d-b)\,V_{x/y}(d-c)   \;=\;
W_{z,y}(c-b)\, V_{x/z}(a-c)\, W_{y,x}(a-b)\,.\label{str2a}
\end{equation}
\end{subequations}
The above two relations are equivalent (see 
\eqref{str-tr}), so it is enough to
consider one of them. For definiteness 
we choose the second one \eqref{str2a}. Due to \eqref{vanish} both
sides of this relation vanish identically unless $a\ge c$. So in the
following we assume $a\ge c$. Next, we change the summation
variable there from $d$ to $n=a-d$, which takes the values
in the interval $0\le n\le a-c$ (since otherwise the summand
in \eqref{str2a} vanish). Using the above definitions one can 
transform the LHS of \eqref{str2a} as follows
\begin{subequations}
\begin{eqnarray}
\mbox{LHS\ of\ \eqref{str2a}}=&&\nonumber\\[.4cm]
&&\hskip -3cm
=V_{x/y}(a-c)\,W_{z,x}(a-b)\phantom{|}_3\varphi_2
\left(\begin{array}{ccc|}q^{-2(a-c)},&(y/z)^2,&q^{2(1-a+b)} x^{-2}\\[.2cm]
&q^{2(1-a+c)}(y/x)^2,&q^{2(1-a+b)} z^{-2}
\end{array}\,q^2,q^2\right)\label{pf1}\\[.3cm]
&&\hskip -3cm
=\displaystyle V_{x/y}(a-c)\,W_{z,x}(a-b)\,
 \frac{\left(q^{2(1-a+c)}\,(z/x)^2,\,q^{2(c-b)}\,y^2\,;q^2\right)_{a-c}}
{\left(q^{2(1-a+c)}\,(y/x)^2,\,q^{2(c-b)}\,x^2\,;q^2\right)_{a-c}}\label{pf2}
\\[.4cm]
&&\hskip -3cm
=W_{z,y}(c-b)\, V_{x/z}(a-c)\, W_{y,x}(a-b)\,,\label{pf3}
\end{eqnarray}
\end{subequations}
where in the first line we used \eqref{Vdef}, \eqref{Wxy} and
\eqref{cor} as well as the formula 
\begin{equation}\label{Wcor}
W_{x,y}(a-n) \;=\; W_{x,y}(a) \left(\frac{y}{x}\right)^n
\frac{(q^{2-2a}\,y^{-2};q^2)_n}
{(q^{2-2a}\,x^{-2};q^2)_n}\;.
\end{equation}
In the next line \eqref{pf2} we simply 
substituted the summation formula \eqref{pfaff}. Finally 
Eqs.\eqref{cor} and \eqref{Wcor} were  used again to rewrite 
the result in the form 
\eqref{pf3}, which is identical to the RHS of \eqref{str2a}, thus
completing the proof of this relation.  

\sapp{Star-star relation.} 
Next, we will show how to derive the star-star relation \eqref{ssrel} from
the star-triangle relation. Setting $y=q$ in  \eqref{stra}
one obtains
\begin{subequations}
\label{tr}
\begin{eqnarray}
\frac{V_{z}(b-c)\, V_{x/z}(c-a)}{V_x(b-a)}&=&
\sum_k V_{q/z}(k-a)\, W_{z,x}(b-k)\, V_{x/q} (c-k)\,,\label{tr1}\\[.3cm]
&=& 
\sum_k W_{x/z,x}(k-a)\, V_{qz/x}(b-k)\,  V_{x/q} (k-c)\,,\label{tr2}
\end{eqnarray} 
\end{subequations}
where the first equality follows from \eqref{str1a} while the second
one follows from \eqref{str2a} after an additional substitution
$a\leftrightarrow b$ and $z\to x/z$.

The LHS of \eqref{ssrel} can be transformed as 
\begin{equation}
\begin{array}{l}
\displaystyle \mbox{LHS\ of\ \eqref{ssrel}}=\displaystyle \sum_n 
\left(\frac{\displaystyle V_{y'/x'}({n-b})\, V_{y/y'}(b-d)}
{\displaystyle V_{y/x'}(n-d)}\right)\;
\frac{\displaystyle V_{x/y'}(a-n)\, V_{y/x}(n-c)}{\displaystyle V_{y/y'}(a-c)}\\
[6mm]
\displaystyle =\;\sum_{k,n} V_{qx'/y'}(k-d)\,V_{y/qx'}(b-k)\,
W_{y'/x',\,y/x'}(n-k)\,\,
\frac{\displaystyle V_{x/y'}(a-n)\, V_{y/x}(n-c)}{\displaystyle V_{y/y'}(a-c)}\\
[8mm]
\displaystyle =\;\sum_k V_{qx'/y'}(k-d)\,V_{y/qx'}(b-k)\, 
W_{x/x',y/x'}(a-k) \, W_{y'/x',x/x'}(c-k)\;.
\end{array}\label{SS1}
\end{equation}
where at first \eqref{tr1} is applied to the product of three $V$'s
enclosed by the parenthesis and
then the summation over $n$ is performed with help of \eqref{str2a}. 
Consider now the RHS of \eqref{ssrel}, given by \eqref{W-def2}, 
\begin{equation}
\begin{array}{l}
\displaystyle \mbox{RHS\ of\ \eqref{ssrel}}
=\displaystyle \sum_n \left(\frac{\displaystyle V_{x/x'}(a-b) \,V_{y/x}(b-n)}{\displaystyle
  V_{y/x'}(a-n)}\right)\;
\left(\frac{\displaystyle V_{y'/x'}(c-n)\, V_{x/y'}(n-d)}{\displaystyle V_{x/x'}(c-d)}\right)\\
[6mm]
\displaystyle =\;\sum_{n,k,k'} V_{q x'/x}(k-n)\, V_{y/q x'}(b-k)\,
W_{x/x',y/x'}(a-k) 
\,V_{qx'/y'}(k'-d) \, V_{x/q x'}(n-k') \, W_{y'/x',x/x'}(c-k')\\
[6mm]
\displaystyle =\; \sum_k \,
V_{y/q x'}(b-k)\,
W_{x/x',y/x'}(a-k) 
\,V_{qx'/y'}(k-d) \,  W_{y'/x',x/x'}(c-k)\,,
\end{array}\label{SS4}
\end{equation}
where at first \eqref{tr1} is applied to both products of three $V$'s
grouped in the parentheses and then the summation over $n$ is
removed by using the inversion relation \eqref{inv}. The resulting
expression coincides with \eqref{SS1}, thereby proving \eqref{ssrel}.
Thus the star-triangle relation \eqref{str} 
implies the star-star relation \eqref{ssrel} and, consequently, the
Yang-Baxter equation \eqref{ybe} and the commutativity of the transfer
matrices \eqref{tcomm}.
Therefore, the star-triangle relation
\eqref{str} does, indeed, plays the role of a primary integrability condition
in the model.


\sapp{Alternative graphical notations.}
To present the star-triangle relation graphically we need to introduce
a new variant of graphical notations for the edge Bolzmann
weights. They involve the edges of the type $(i)$ and type $(iii)$
\begin{equation}\label{type3}
\input{fig2dc.tex}
\end{equation}
which are distiguished by relative orientation of the edge and associated
spectral parameter (dashed) lines. 
The corresponding Boltzmann weights $V_{x/y}(a-b)$ 
and $W_{x,y}(a-b)$ are defined by \eqref{Vdef} and \eqref{Wxy}, respectively. 
Note that for the type $(iii)$ edges the directed rapidity lines
cross the edge from the same side, unlike the type $(i)$ edges where they
come from different sides. For better visibility the type $(iii)$
edges will be shown by double lines. 

As in \eqref{twotypes} the weights in \eqref{type3} 
depend on the difference of spins at the ends of the edge, but, in
general, they do not have the spectral paramerer ``difference
property''. The function $W$ 
depends non-trivially on both rapidities $x$ and $y$,
but just on their ratio,
With the graphical rules \eqref{type3} the first
relation in \eqref{str} is represented graphically as in Fig.~\ref{figstr}. 
\begin{figure}[ht]
\begin{center}
\input{strpic.tex}
\end{center}
\caption{Graphical representation of the star-triangle relation
  \eqref{str1}.
\label{figstr}}
\end{figure}
For 
second relation one needs to reverse all arrows (both edge arrows and
the dashed line arrows).

\sapp{New solutions of the Yang-Baxter equations\label{sappb}}
The existence of the star-triangle relation \eqref{str} 
allows one to define new integrable Ising-type models in the standard
way. First, let us present them in the IRF form. To make the notations
uniform it is convenient to redenote weights \eqref{W-def1} as
\begin{equation}
\mathcal{W}^{(1,1)}(a,b,c,d\,|\,\boldsymbol{x},\boldsymbol{y})\;=\;
\mathcal{W}(a,b,c,d\,|\,\boldsymbol{x},\boldsymbol{y})
\end{equation}
where $\mathcal{W}$ is given by \eqref{W-def1} and 
${\boldsymbol x}=(x,x')$, ${\boldsymbol y}=(y,y')$ denote pairs of spectral
variables. 

Define new weight functions,
\begin{equation}
\mathcal{W}^{(2,1)}(a,b,c,d\,|\,\boldsymbol{x},\boldsymbol{y})\;=\;
\theta(a\ge c)\,\sum_{n}
W_{x,y}(a-n) V_{y/x'}(b-n) W_{y',x}(c-n)  V_{x'/y'}(n-d)\,,
\end{equation}
\begin{equation}
\mathcal{W}^{(1,2)}(a,b,c,d\,|\,\boldsymbol{x},\boldsymbol{y})\;=\;
\theta(a\ge b)\,\sum_{n}
W_{x,y}(a-n) W_{y,x'}(b-n) V_{y'/x}(n-c)  V_{x'/y'}(d-n)\,,
\end{equation}
and
\begin{equation}
\mathcal{W}^{(2,2)}(a,b,c,d\,|\,\boldsymbol{x},\boldsymbol{y})\;=\;\sum_{n}
W_{x,y}(a-n) V_{y/x'}(b-n) V_{y'/x}(n-c) W_{x',y'}(n-d)\;.
\end{equation}
Here $\theta(a\ge b)=1$ if $a\ge b$ 
and vanishes otherwise. 
Using the notations \eqref{type3} we can represent these weights graphically
\newcommand{\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);}{\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);}
\begin{equation}\label{W1}
\begin{tikzpicture}[scale=0.8,baseline=(current  bounding  box.center)]
\draw [-open triangle 45, black, thin, dashed] (1,-0.5) -- (1,4.5);
\draw [-open triangle 45, black, thin, dashed] (3,-0.5) -- (3,4.5);
\draw [-open triangle 45, black, thin, dashed] (-0.5,1) -- (4.5,1);
\draw [-open triangle 45, black, thin, dashed] (-0.5,3) -- (4.5,3);
\begin{scope}[very thick, every node/.style={sloped,allow upside down}]
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (4,0)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (2,2);
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (2,2)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (0,4);
\end{scope}
\draw [-latex, blue!60!black, ultra thick] (0,0) -- (1,1); 
\draw [-latex, blue!60!black, ultra thick] (1,1) -- (3,3);
\draw [blue!60!black, ultra thick] (3,3) -- (4,4);  
\draw [white, fill, ultra thick] (0,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (2,2) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (2,2) circle [radius=0.09];
\draw [white, fill, ultra thick] (0,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,4) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,4) circle [radius=0.09];
\node [above] at (1,4.5) {$x$};
\node [above] at (3,4.5) {$x'$};
\node [right] at (4.5,3) {$y$};
\node [right] at (4.5,1) {$y'$};
\node [above] at (0,4.1) {$a$};
\node [above] at (4,4.1) {$b$};
\node [below] at (0,-0.1) {$c$};
\node [below] at (4,-0.1) {$d$};
\node [left] at (-1,2) {$\displaystyle \mathcal{W}^{(2,2)}(a,b,c,d| x,x';y,y')\;=$};
\end{tikzpicture}
\end{equation}

\begin{equation}\label{W2}
\begin{tikzpicture}[scale=0.8,baseline=(current  bounding  box.center)]
\draw [-open triangle 45, black, thin, dashed] (1,-0.5) -- (1,4.5);
\draw [-open triangle 45, black, thin, dashed] (3,-0.5) -- (3,4.5);
\draw [-open triangle 45, black, thin, dashed] (-0.5,3) -- (4.5,3);
\draw [-open triangle 45, black, thin, dashed] (4.5,1) -- (-0.5,1);
\begin{scope}[very thick, every node/.style={sloped,allow upside down}]
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (2,2)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (0,0);
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (2,2)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (0,4);
\end{scope}
\draw [-latex, blue!60!black, ultra thick] (4,0) -- (3,1);
\draw [blue!60!black, ultra thick] (3,1) -- (2,2); 
\draw [-latex, blue!60!black, ultra thick] (2,2) -- (3,3);
\draw [blue!60!black, ultra thick] (3,3) -- (4,4);  
\draw [white, fill, ultra thick] (0,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (2,2) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (2,2) circle [radius=0.09];
\draw [white, fill, ultra thick] (0,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,4) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,4) circle [radius=0.09];
\node [above] at (1,4.5) {$x$};
\node [above] at (3,4.5) {$x'$};
\node [right] at (4.5,3) {$y$};
\node [right] at (4.5,1) {$y'$};
\node [above] at (0,4.1) {$a$};
\node [above] at (4,4.1) {$b$};
\node [below] at (0,-0.1) {$c$};
\node [below] at (4,-0.1) {$d$};
\node [left] at (-1,2) {$\displaystyle \mathcal{W}^{(2,1)}(a,b,c,d| x,x';y,y')\;=$};
\end{tikzpicture}
\end{equation}

\begin{equation}\label{W3}
\begin{tikzpicture}[scale=0.8,baseline=(current  bounding  box.center)]
\draw [-open triangle 45, black, thin, dashed] (1,-0.5) -- (1,4.5);
\draw [-open triangle 45, black, thin, dashed] (3,4.5) -- (3,-0.5);
\draw [-open triangle 45, black, thin, dashed] (-0.5,3) -- (4.5,3);
\draw [-open triangle 45, black, thin, dashed] (-0.5,1) -- (4.5,1);
\begin{scope}[very thick, every node/.style={sloped,allow upside down}]
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (2,2)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (4,4);
\draw[double=none, double distance=0.5pt, line join=round, very thick, blue!60!black] (2,2)--node [pos=.4] {\tikz \draw[-{Latex[length=3.5mm]}](0,0) -- +(.1,0);} (0,4);
\end{scope}
\draw [-latex, blue!60!black, ultra thick] (0,0) -- (1,1);
\draw [blue!60!black, ultra thick] (1,1) -- (2,2); 
\draw [-latex, blue!60!black, ultra thick] (2,2) -- (3,1);
\draw [blue!60!black, ultra thick] (3,1) -- (4,0);  
\draw [white, fill, ultra thick] (0,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (2,2) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (2,2) circle [radius=0.09];
\draw [white, fill, ultra thick] (0,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (0,4) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,0) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,0) circle [radius=0.09];
\draw [white, fill, ultra thick] (4,4) circle [radius=0.09];
\draw [blue!60!black, ultra thick] (4,4) circle [radius=0.09];
\node [above] at (1,4.5) {$x$};
\node [above] at (3,4.5) {$x'$};
\node [right] at (4.5,3) {$y$};
\node [right] at (4.5,1) {$y'$};
\node [above] at (0,4.1) {$a$};
\node [above] at (4,4.1) {$b$};
\node [below] at (0,-0.1) {$c$};
\node [below] at (4,-0.1) {$d$};
\node [left] at (-1,2) {$\displaystyle \mathcal{W}^{(1,2)}(a,b,c,d| x,x';y,y')\;=$};
\end{tikzpicture}
\end{equation}
There are eight different Yang-Baxter equations for $\mathcal{W}^{(i,j)}$:
\begin{equation}
\begin{array}{l}
\displaystyle \sum_{a} 
\mathcal{W}^{(i,j)}(g,a,b,f \,|\, \boldsymbol{x},\boldsymbol{y}) \,
\mathcal{W}^{(i,k)}(c,e,g,a \,|\, \boldsymbol{x},\boldsymbol{z}) \,
\mathcal{W}^{(j,k)}(e,d,a,f \,|\, \boldsymbol{y},\boldsymbol{z}) \;=\\
[5mm]
\displaystyle \qquad \qquad
=\; \sum_{h} 
\mathcal{W}^{(j,k)}(c,h,g,b \,|\, \boldsymbol{y},\boldsymbol{z})\,
\mathcal{W}^{(i,k)}(h,d,b,f \,|\, \boldsymbol{x},\boldsymbol{z})\,
\mathcal{W}^{(i,j)}(c,e,h,d \,|\, \boldsymbol{x},\boldsymbol{y})
\end{array}
\end{equation}
for all possible choices of the indices $i,j,k=1,2$. 

All four IRF weights $\mathcal{W}^{(i,j)}$ 
allow a corresponding vertex form $\mathcal{R}^{(i,j)}$ 
by means of \eqref{R-def} \eqref{ij-def}. From the quantum group point
of view these $R$-matrices are related with infinite-dimensional
evaluation representations of the $U_q(\widehat{sl}(2))$ algebra. The
indices $i,j=1$ in the notation $\mathcal{R}^{(i,j)}$
indicate representations with highest weights, while 
$i,j=2$ indicate representations without highest or lowest weights.
As in Sect.~\ref{red} the $R$-matrices, associated with the 
highest weight representations admit
finite-dimensional reductions. For instance, 
when $y'=q\,y$, the $R$-matrix 
$\mathcal{L}\;=\;\mathcal{R}^{(2,1)}(x,x';y,qy)$ has a
block-triangular structure in the second space with a $2$-dimensional
block given by 
\begin{equation}\label{Lop-inf}
\begin{array}{ll}
\displaystyle \mathcal{L}_{i_1,0}^{j_1,0}\;=\;\delta_{i_1,j_2} \frac{[q^{i_1}y]}{[y]}\;,
& \displaystyle \mathcal{L}_{i_1,0}^{j_1,1}\;=\;-\delta_{i_1,j_1+1} \frac{[q^{j_1} x']}{[y]}\;,\\
[7mm]
\displaystyle \mathcal{L}_{i_1,1}^{j_1,0}\;=\;\delta_{i_1+1,j_1}\frac{[q^{i_1}x]}{[y]}\;, &
\displaystyle \mathcal{L}_{i_1,1}^{j_1,1}\;=\;-\delta_{i_1,j_1} \frac{[q^{i_1-1}xx'/y]}{[y]}\;.
\end{array}
\end{equation}
It can be written in the canonical form (the same as \eqref{Lop}) 
\begin{equation}
\boldsymbol{{\cal L}}=\frac{1}{[\mu\,q^{s+s'}]}\,
\begin{pmatrix}
\mu\, q^{\frac{H}{2}}-\mu^{-1} \,q^{-\frac{H}{2}}&
F\\[3mm]
E& \mu \,q^{-\frac{H}{2}}-\mu^{-1} \,q^{\frac{H}{2}}
\end{pmatrix}\,,\label{Lop1}
\end{equation}
where $E,F,H$ stands for the generators of 
the quantum universal enveloping algebra $U_q(sl_2)$, defined in
\eqref{uqsl2},
and the matrix elements \eqref{Lop-inf} correspond to the
infinite-dimensional representation $\pi^{+}_{s,s'}$ of this algebra
\begin{equation}
\pi^+_{s,s'}[H]\,|j\rangle=(2s+2s'-2j)\,|j\rangle\,,
\quad \pi^+_{s,s'}[E]\,|j\rangle=[q^{2s+1-j}]\,|j-1\rangle\,,\quad 
\pi^+_{s,s'}[F]\,|j\rangle=[q^{j+1-2s'}]\,|j+1\rangle\,,
\end{equation}
where $j\in {\mathbb Z}$ and the parameters are defined as 
\begin{equation}
y=\big(\mu q^{s+s'}\big)^{-1}\,,\qquad x=q^{-2s}\,,\qquad x'=q^{1-2s'}\,.
\end{equation}
Note, that for the generic case this representation is irreducible and 
does not have highest or lowest weights.  





\app{Mangazeev's $R$-matrix for higher spin 6-vertex model\label{appC}}
\addcontentsline{toc}{section}{C. Mangazeev's $R$-matrix for higher
  spin 6-vertex model} 
In ref.\cite{Mangazeev:2014gwa} Mangazeev obtained a remarkably simple
formula for the most general higher spin $R$-matrix related to
6-vertex model, expressing it via ${}_4\varphi_3$ $q$-hypergeometric series.
The result of \cite{Mangazeev:2014gwa} 
can be represented in various equivalent forms. Here we
will use its variant from \cite{Bosnjak:2016oze}, eq.(3.22) therein,
which we quote below preserving their notations 
\begin{equation}\label{R-Man}
\begin{array}{lll}
\displaystyle \bigl[R_{I,J}(\lambda)\bigr]_{i,j}^{i',j'} &=& \displaystyle \delta_{i+j,i'+j'} q^{i'j'-ij-iJ-Ij'}\, \frac{(q^2,q^2)_{i+j}}{(q^2;q^2)_i\,(q^2,q^2)_j}\\
[6mm]
&\displaystyle \times &\displaystyle  \frac{(\lambda^{-2}q^{I-J};q^2)_{j'}\, 
(\lambda^{-2}q^{J-I};q^2)_i\, (q^{-2J};q^2)_j}
{(\lambda^{-2}q^{-I-J};q^2)_{i+j}\,(q^{-2J};\,q^2)_{j'}}\\
[6mm]
&\times & \displaystyle \sum_{n\geq 0} \frac{\displaystyle (q^{-2i},\, q^{-2j'},\, 
\lambda^2q^{-I-J},\, \lambda^2q^{2+I+J-2i-2j};q^2)_n}{\displaystyle (q^2,\, 
q^{-2i-2j},\, \lambda^2q^{2+I-J-2i},\, \lambda^2q^{2+J-I-2j'};q^2)_n}\, q^{2n}\,.
\end{array}
\end{equation}
Apart from the parameter $q$, this expressions contains three
independent variables $q^I$, $q^J$ and $\lambda^2$.
They only appear in certain combinations, so it is convenient 
to introduce the variables $x,x',y,y'$, such that 
\begin{equation}
\left(\frac{x}{y'}\right)^2\;=\;\lambda^2 q^{-I-J}\;,\quad 
\left(\frac{x'}{y}\right)^2\;=\;\lambda^2 q^{I+J}\;,\quad 
\left(\frac{x}{y}\right)^2\;=\;\lambda^2 q^{J-I}\;,\quad 
\left(\frac{x'}{y'}\right)^2\;=\;\lambda^2 q^{I-J}\;.
\end{equation}
Clearly, there are  
only three independent ratios of these new variables,
but the fourth equation above is a corollary of the other three.
Next, let us make a substitution of the integer indices in \eqref{R-Man},
\begin{equation}
i=a-b\;,\quad j= b-d\;,\quad i'=c-d\;,\quad j'= a-c\,,\qquad
a,b,c,d\in {\mathbb Z}\,,
\end{equation}
which automatically solves the $\delta$-function constraint. 
Further, let us combine all $q$-product factors in front of the sum in
\eqref{R-Man} into a product of the weight functions $V_x(n)$, using the
definition \eqref{Vdef}, and then write the sum in
\eqref{R-Man} as ${}_4\varphi_3$ series using \eqref{phi43}. Then,
after all these substitutions, the resulting expression for \eqref{R-Man}
almost literally coincides with the RHS of \eqref{SL}. More precisely,
one obtains,
\begin{equation}
\mbox{\Big[RHS of \eqref{R-Man}\Big]}=q^{i'j'-ij}\, 
\big(\lambda q^{(I-J)/2}\big)^{i'-i}\,\times\,
\mbox{\Big[RHS of \eqref{SL}\Big]}\,.
\end{equation}
Taking now into account the definitions \eqref{R-def}, \eqref{pardef},
\eqref{RR-rel} and the fact that \eqref{SL} coincides with \eqref{star1} 
one concludes that the $R$-matrix \eqref{R-Man} is simply related to that
defined by \eqref{RR-rel}, \eqref{R-def2},  
\begin{equation}
\bigl[R_{I,J}(\lambda)\bigr]^{i_1,i_2}_{j_1,j_2}=
q^{(j_1j_2-i_1i_2)}\,\big(q^{s_1-s_2}\,\lambda\big)^{i_1-j_1}\,
{\cal R}(\lambda\,\,|\,\,s_1,s_2)_{i_1,i_2}^{j_1,j_2}
\end{equation}
where
\begin{equation}
I=2s_1\,,\qquad J=2s_2\,,\qquad
\end{equation}
and $\lambda$ is the same in both sides. 
So the two $R$-matrices are related by the matrix
transpostion in both spaces, followed by a multiplication by 
some gauge transformations factors. Note, that the latter 
cancel out the Yang-Baxter equation. The reason we have dropped these
factors is that it is impossible to 
consistently attribute them to the edge Boltzmann weights
in the Ising type formulation of the model.


Let now present the connection with Mangazeev's definitions of the
${\bf Q}$-operators. First compare two definitions, our
Eq.\eqref{Qelem3}, defining
$\overline{\mathfrak{R}}_{i_1,i_2}^{j_1,j_2}$, and Eq.(6.21) of
\cite{Mangazeev:2014gwa}, which defines 
$[A_{-}^{(I)}(\lambda)]_{\,n,\,i}^{n',i'}$. 
These quantities are related as follows,
\begin{equation}\label{1}
[A_{-}^{(I)}(\lambda)]_{\,n,\,i}^{n',i'}\;=\;(-\lambda)^{I-i'} \; ((y/x')^2;q^2)_{I} \; \overline{\mathfrak{R}}(y|x,x')_{i',n'}^{\,i,\,n} \; q^{n'i'-ni}\;\frac{c_{n'}}{c_n}\;,
\end{equation}
where 
\begin{equation}\label{22}
\left(\frac{y}{x'}\right)^2\;=\;\lambda^{-2}q^{1-I}\;,\quad
\left(\frac{y}{x}\right)^2\;=\;\lambda^{-2}q^{1+I}\;,\quad
c_n\;=\;q^{-n(1+I)/2} (q^2;q^2)_n\;.
\end{equation}
Extra power $(-\lambda)^{I-i'}$ in (\ref{1}) appear because in Mangazeev' formulas $Q_{-}$ is the Laurent polynomial. 

Thus, Mangazeev's $Q_{-}(\lambda)$ is related to $\overline{\mathbf{Q}}(q^{-1}y)$ by
\begin{equation}
Q_{-}(\lambda) \;=\; \lambda^{NI-M} \overline{\mathbf{Q}}(q^{-1}y)\;.
\end{equation}
for integer $I=2s_1$.



\app{{\bf Q}-matrices and functional relations\label{appD}}
\addcontentsline{toc}{section}{D. {\bf Q}-matrices and functional relations}
The elements of the ${\bf Q}$-matrices \eqref{Q-def} 
\begin{equation}
\label{Qelem}
\begin{array}{rcl}
\displaystyle \big({\bf Q}(y')\big)_{i_N,\ldots,i_1}^{j_N,\ldots,j_1}&=&
\displaystyle Z_0^{-1}\,\sum_{k_1=0}^\infty \omega_h^{k_1} \prod_{m=1}^N
{\mathfrak R}(y'|x_m^{},x'_m)_{i_m,k_{m+1}}^{j_m,k_m}\;,\\[.5cm]
\displaystyle
\varphi(y,y)^{-1}\big(\overline{{\bf Q}}(q^{-1}\,y)\big)_{i_N,\ldots,i_1}^{j_N,\ldots,j_1}&=&
\displaystyle Z_0^{-1}\,\sum_{k_1=0}^\infty \omega_h^{k_1} \prod_{m=1}^N 
 \overline{{\mathfrak R}}(y|x_m^{},x'_m)_{i_m,k_{m+1}}^{j_m,k_m}\;,
\end{array}
\end{equation}
where $k_{N+1}=k_1$, are expressed through the 
$R$-matrices, 
\begin{equation}\label{Qelem2}
\begin{array}{l}
\displaystyle{\mathfrak R}(y'\,|\,x,x')_{i_1,i_2}^{j_1,j_2}=
\displaystyle q^{i_2-j_2}\,\frac{(q^2;q^2)_{j_2}}{(q^2;q^2)_{i_2}}\,
\lim_{y\to0}\Big(\frac{y'}{y}\Big)^{j_1}\,{\cal R}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}
\\[7mm]
\qquad=\displaystyle\delta_{i_1+i_2,j_1+j_2} \left(q\frac{x'}{y'}\right)^{-i_1}
\sum_{n=0}^{\min(i_2,j_1)}
\left(\frac{x'}{x}\right)^{n-i_2} \frac{(q^2;q^2)_{i_1+i_2-n}}{(q^2;q^2)_{i_2-n}}
V_{x/y'}(n) V_{y'/x'}(j_1-n)
\end{array}
\end{equation}
and 
\begin{equation}\label{Qelem3}
\begin{array}{l}
\displaystyle\overline{{\mathfrak R}}(y\,|\,x,x')_{i_1,i_2}^{j_1,j_2}=
\displaystyle (-)^{i_2-j_2}\,q^{i_2^2-j_2^2}\,\frac{(q^2;q^2)_{j_2}}{(q^2;q^2)_{i_2}}\,
\lim_{y'\to0}\Big(\frac{y'}{y}\Big)^{i_1}\,{\cal R}(x,x',y,y')_{i_1,i_2}^{j_1,j_2}
\\[7mm]
\qquad=\displaystyle\delta_{i_1+i_2,j_1+j_2}  \left(q\frac{x'}{y}\right)^{j_1}\sum_{n=0}^{\min(i_2,j_1)} 
\left(-q\frac{x'}{x}\right)^{-n} \frac{q^{n^2}}{(q^2;q^2)_{n} (q^2;q^2)_{j_1-n}} \frac{V_{y/x}(i_2-n)}{V_{y/x'}(i_1+i_2-n)}\;.
\end{array}
\end{equation}
obtained by suitable limits of \eqref{R-def2}. 

At $x=q^{-1}x'$, what is $s=1/2$, one has
\begin{equation}\label{Qelem4}
\begin{array}{ll}
\displaystyle \mathfrak{R}_{0,i_2}^{0,j_2}\;=\;\delta_{i_2,j_2} \, q^{-j_2}\;,& 
\displaystyle \mathfrak{R}_{0,i_2}^{1,j_2}\;=\; - \,\delta_{i_2,j_2+1}\, \frac{y'}{x'} \,q^{-j_2}\;,\\
[7mm]
\displaystyle \mathfrak{R}_{1,i_2}^{0,j_2}\;=\; -\, \delta_{i_2+1,j_2} \, \frac{y'}{x'}\, [q^{j_2}]\;,& 
\displaystyle \mathfrak{R}_{1,i_2}^{1,j_2}\;=\;\delta_{i_2,j_2} \, \frac{y'}{x'} \, [ \frac{x'}{y'} q^{j_2} ] \;,
\end{array}
\end{equation} 
and
\begin{equation}\label{Qelem5}
\begin{array}{ll}
\displaystyle \overline{\mathfrak{R}}_{0,i_2}^{0,j_2}\;=\;\delta_{i_2,j_2} \, \frac{\displaystyle  [ \frac{y}{x'} q^{j_2} ]}{\displaystyle [\frac{y}{x'}]}\;,& 
\displaystyle \overline{\mathfrak{R}}_{0,i_2}^{1,j_2}\;=\; - \,\delta_{i_2,j_2+1}\, \frac{\displaystyle q^{j_2}}{\displaystyle  [\frac{y}{x'}]}\;,\\
[7mm]
\displaystyle \overline{\mathfrak{R}}_{1,i_2}^{0,j_2}\;=\; \delta_{i_2+1,j_2} \, \frac{\displaystyle [q^{j_2}]}{\displaystyle [\frac{y}{x'}]} \;,& 
\displaystyle \overline{\mathfrak{R}}_{1,i_2}^{1,j_2}\;=\;-\, \delta_{i_2,j_2} \, \frac{\displaystyle \frac{x'}{y} q^{j_2}}{\displaystyle [\frac{y}{x'} ]} \;,
\end{array}
\end{equation} 
In the short notations
\begin{equation}
\mathfrak{R}\;=\;\left(\begin{array}{cc}
\displaystyle q^{H/2} & \displaystyle \mu \mathcal{E}_{-} \\
[5mm]
\displaystyle \mu \mathcal{E}_{+} & \displaystyle q^{-H/2} - q^{-1}\mu^2 q^{H/2}
\end{array}\right)
\end{equation}
where
\begin{equation}
\mu\;=\;q^{1/2}\frac{y'}{x'}\;,\quad \langle i |q^{H/2}| j\rangle \;=\; \delta_{i,j} q^{-j}\;,
\end{equation}
and
\begin{equation}
\langle i | \mathcal{E}_{-} | j \rangle \;=\; -\delta_{i,j+1} q^{-j-1/2}\;,\quad 
\langle i | \mathcal{E}_{+} | j \rangle \;=\; \delta_{i+1,j} (q^{-j} - q^{j}) q^{-1/2}\;.
\end{equation}
Similarly,
\begin{equation}
\overline{\mathfrak{R}}\;=\;\frac{1}{1-q\mu^2} \,
\left(\begin{array}{cc}
\displaystyle q^{-H/2} - q\mu^2 q^{H/2} & \displaystyle \mu\mathcal{E}_{+}\\
[5mm]
\displaystyle \mu\mathcal{E}_{-} & \displaystyle q^{H/2}
\end{array}
\right)\;,
\end{equation}
where
\begin{equation}
\mu\;=\;q^{-1/2}\frac{y}{x'}\;,\quad \langle i | q^{H/2} | j \rangle \;=\; \delta_{i,j} q^j\;,
\end{equation}
and
\begin{equation}
\langle i | \mathcal{E}_{-} | j \rangle \;=\; \delta_{i+1,j} (q^{-j}-q^j) q^{1/2}\;,\quad 
\langle i | \mathcal{E}_{+} | j \rangle \;=\; -\delta_{i,j+1} q^{j+1/2}\;.
\end{equation}
In both cases
\begin{equation}
[H,\mathcal{E}_{\pm}]\;=\;\pm 2 \mathcal{E}_{\pm}\;,\quad
q\mathcal{E}_{+}\mathcal{E}_{-} - q^{-1} \mathcal{E}_{-}\mathcal{E}_{+}\;=\;q-q^{-1}\;.
\end{equation}


Using the correspondence \eqref{index} one can write the IRF version
of the above operators \eqref{Qelem} as
\begin{equation}\label{Qcal-def}
\begin{array}{l}
\displaystyle
\big({\mathbb
  Q}(
y')\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}=\\[.5cm]
\qquad=
\displaystyle\prod_{i=1}^N \left(\frac{q\,x'_i}{y'}\right)^{a_i-a_{i+1}}\,
\left(\sum_{m_i} \left(\frac{x'_i}{x_i}\right)^{a_{i+1}-m_i}\,
\frac{(q^2;q^2)_{m_i-a_i}\,
V_{x_i/y'}(b_{i+1}-m_i)\, V_{y'/x'_i}(m_i-b_i)}
{(q^2;q^2)_{m_i-a_{i+1}}
}
\right)\,,\\[.7cm]
\displaystyle\big(\overline{{\mathbb
  Q}}(q^{-1}\,
y)\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}=\\[.5cm]
\qquad=\displaystyle
\prod_{i=1}^N \left(\frac{q\,x'_i}{y}\right)^{b_{i+1}-b_i}
\left(
 \sum_{n_i}
\left(-q\frac{x'_i}{x_i}\right)^{n_i-b_{i+1}}\,
\frac{q^{(b_{i+1}-n_i)^2}\,V_{y/x_i}(n_i-a_{i+1})}
{(q^2;q^2)_{b_{i+1}-n_i}\, (q^2;q^2)_{n_i-b_i}\, V_{y/x'_i}(n_i-a_i)}
\right)\,,
\end{array}
\end{equation}
where the summations run over the intervals 
\begin{equation}\label{intervals}
\max(a_{i+1},b_i)\le
n_i\le b_{i+1}\,,\qquad \max(a_{i+1},b_i)\le
m_i\le b_{i+1}\,,
\end{equation}
and 
\begin{equation}
a_{N+1}=M+a_1\,,\qquad b_{N+1}=M+b_1\,,
\end{equation}
with $M$ being
the conserved number of up spins \eqref{down2}. 
With the same correspondence between indices \eqref{index} the matrix
elements \eqref{Qelem} can now be written as
\begin{equation}
\begin{array}{rcl}
\displaystyle \big({\bf Q}( y')\big)_{i_N,\ldots,i_1}^{j_N,\ldots,j_1}&=&
\displaystyle Z_0^{-1}\,\sum_{a_1\le b_1\le \infty}
\omega_h^{b_1-a_1} \,
\big({\mathbb
  Q}(
y')\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}\\[.7cm]
\displaystyle \varphi(y,y)^{-1}\big(\overline{{\bf Q}}(q^{-1}\, y)\big)_{i_N,\ldots,i_1}^{j_N,\ldots,j_1}&=&
\displaystyle Z_0^{-1}\,\sum_{a_1\le b_1\le \infty}
\omega_h^{b_1-a_1} \,
\big(\overline{{\mathbb
  Q}}(q^{-1}\,
y )\big)_{a_{N+1},a_{N},\ldots,a_{1}}^{b_{N+1},\,b_{N},\ldots,b_{1}}
\end{array}
\end{equation}
where the sum is taken over the spin $b_1$. Consider now the product
\begin{equation}\label{QQ-prod}
\begin{array}{l}
\displaystyle \varphi(y,y)^{-1}
\big(\overline{{\bf
  Q}}(q^{-1}\,
y)\,{{\bf
  Q}}(
y')\big)_{i_N,\ldots,i_1}^{j_N,\ldots,j_1}=\\[.5cm]
\qquad=\displaystyle
Z_0^{-2}\sum_{a_1\le b_1\le \infty}\,
\omega_h^{b_1-a_1}\sum_{c_1,c_2,\ldots,c_N}
\overline{{\mathbb
  Q}}(q^{-1}\,
y)_{a_{N+1},\ldots,a_{1}}^{c_{N+1},\ldots,c_{1}}
\,
{\mathbb
  Q}(
y')_{c_{N+1},\ldots,c_{1}}^{b_{N+1},\ldots,b_{1}}\,.
\end{array}
\end{equation}
Note that the elements \eqref{Qcal-def} have the same vanishing
conditions as those of the IRF transfer matrix
\eqref{nonzero}. Therefore, the summation over $c_1,c_2,\ldots,c_N$ is
restricted to a finite domain 
\begin{equation}
c_{i+1}\ge c_i\,,\qquad a_i\le c_i\le b_i\,,\qquad i=1,2,\ldots,N\,. 
\end{equation}
where $c_{N+1}=M+c_1$. Using \eqref{Qcal-def} one obtains 
\begin{equation}
\mbox{\eqref{QQ-prod}}=
\displaystyle
Z_0^{-2}\left(\frac{y'}{y}\right)^M\sum_{b_1=a_1}^\infty\,
\omega_h^{b_1-a_1}
\sum_{\{n_i,m_i\}}\left(\prod_{i=1}^N
\,\frac{V_{x_i/y'}(b_{i+1}-m_i)\, V_{y'/x_i'}(m_i-b_i)\,
  V_{y/x_i}(n_i-a_{i+1})}{V_{y/x'_i}(n_i-a_i)}G_i \right),
\label{QQ-prod2}
\end{equation}
where 
\begin{equation}\label{G-def}
{G}_i\;=\;
\left(\frac{x'_{i-1}}{x_{i-1}}\right)^{(n_{i-1}-m_{i-1})} 
\ \sum_{c_i=n_{i-1}}^{\min(n_i,m_{i-1})}
\frac{q^{(n_{i-1}-c_i)^2}\,(-q)^{n_{i-1}-c_i}\, 
(q^2;q^2)_{m_i-c_i}}{(q^2;q^2)_{c_i-n_{i-1}} (q^2;q^2)_{n_i-c_i}
  (q^2;q^2)_{m_{i-1}-c_i}}\,,
\end{equation}
with the boundary conditions 
\begin{equation}\label{n-bound}
n_0=n_N-M\,,\qquad m_0=m_N-M\,.
\end{equation}
One can show that for all the non-vanishing contributions
to \eqref{QQ-prod2} the spins variables $\{n_i\}$ and $\{m_i\}$ are
ordered as
\begin{equation}\label{m-order}
m_{i+1}\ge n_i\,,\qquad n_{i+1}\ge n_i\,,\qquad m_i\ge n_i\,.
\end{equation}
The considerations are similar to those presented 
in Sect.~\ref{irf-form}, where the ordering of white spins in
Fig.~\ref{fig-lattice} was discussed. Note that spins $\{n_i\}$ and
$\{m_i\}$ correspond to black sites  in the same figure. With 
the substitution 
\begin{equation} \label{abc}
{\mathsf a}=n_i-n_{i-1}\,,\qquad {\mathsf b}=m_{i-1}-n_{i-1}\,,\qquad {\mathsf c}=m_i-n_{i-1}\,,
\end{equation}
the sum in \eqref{G-def} can be written as 
\begin{equation}\label{tozh1}
{\mathsf F}({\mathsf a},{\mathsf b},{\mathsf c})=\sum_{n=0}^{\min({\mathsf a},{\mathsf b})}  \frac{(-)^n
  q^{n(n-1)}\,
(q^2;q^2)_{{\mathsf c}-n}}{(q^2;q^2)_n(q^2;q^2)_{{\mathsf a}-n}(q^2;q^2)_{{\mathsf b}-n}}\;=\;  
q^{2{\mathsf a}{\mathsf b}} 
\frac{(q^2;q^2)_{{\mathsf c}-{\mathsf a}} (q^2;q^2)_{{\mathsf c}-{\mathsf b}}}{(q^2;q^2)_{\mathsf a} (q^2;q^2)_{\mathsf b}
  (q^2;q^2)_{{\mathsf c}-{\mathsf a}-{\mathsf b}}}\;. 
\end{equation}
It is evaluated exactly by using the Heine's $q$-analog
of the Gauss summation formula for the hypergeometric function
\cite{Gasper:2004}.  
It is not difficult to see that
\begin{equation}\label{tozh}
{\mathsf F}({\mathsf a},{\mathsf b},{\mathsf c})\;=\;
\left\{\begin{array}{ll}
0\;,\qquad\qquad\quad &\textrm{if}\quad \max({\mathsf a},{\mathsf b})\leq {\mathsf c} <
{\mathsf a}+{\mathsf b}\,,\\[.1cm] 
\displaystyle q^{2{\mathsf a} {\mathsf b}}\;,\quad &\textrm{if}\quad {\mathsf c}={\mathsf a}+{\mathsf b}\,,\\[.1cm] 
\textrm{non-zero}\;,\quad &\textrm{otherwise}\,.
\end{array}\right.
\end{equation}
Note that with the substitution \eqref{abc} 
the first inequality $\max({\mathsf a},{\mathsf b})\leq {\mathsf c}$ in \eqref{tozh}
is satisfied for all sums in \eqref{G-def} by virtue of \eqref{m-order}.
Let us now examine the difference $\Delta={\mathsf c}-{\mathsf a}-{\mathsf b}$. In terms of the
variables \eqref{abc} it reads 
\begin{equation}
\Delta_i=(m_i-n_i)-(m_{i-1}-n_{i-1})\,.
\end{equation}
Evidently, with the boundary condition \eqref{n-bound} one gets 
\begin{equation}
\sum_{i=1}^N \, \Delta_i=0\,.
\end{equation}
Therefore, either $\Delta_i=0$ for all $i=1,2,\ldots,N$ or at least
one $\Delta_i<0$ (for some particular value of $i$). In the latter
case the corresponding $G_i$ vanishes and so does the whole product in
\eqref{QQ-prod2}. Thus the non-vanishing contributions to
\eqref{QQ-prod2} could only arise when $\Delta_i=0$ for all
$i=1,2,\ldots,N$. This gives $N-1$ independent constraints between the
spins $\{n_i\}$ and $\{m_i\}$ which allow one to easily perform the
summation over all but one spins $\{m_i\}$. Then, instead of $m_1$ we choose the integer 
variable 
\begin{equation}
\delta=m_1-n_1\,, \qquad 0\le \delta \le b_1-a_1,
\end{equation}
whose allowed values are deternimed by the inequalities
\eqref{intervals} and \eqref{m-order}. Eq.\eqref{QQ-prod2} then
takes the form
\begin{equation}
\mbox{\eqref{QQ-prod}}=\displaystyle
Z_0^{-2}\left(\frac{y'}{y}\right)^M\sum_{b_1=a_1}^\infty
\omega_h^{b_1-a_1} 
\displaystyle\sum_{\delta=0}^{b_1-a_1} w^\delta \prod_{i=1}^N 
{\cal W}(b_{i+1},b_i,a_{i+1}+\delta,a_i+\delta\,|\,x_i,x'_i\,|\,y,y')
\label{QQ-prod3}
\end{equation}
where we used \eqref{W-def1} and
\begin{equation}
Z_0=\sum_{k=0}^\infty\,(w\,\omega_h)^k\,,\qquad 
w\;=\; q^{2M}\, \prod_{i=1}^N\, \frac{x_i}{x_i'}\,
\end{equation}
Remembering the definitions \eqref{R-def}, \eqref{t-vert}, \eqref{t-norm}), using a simple identity 
\begin{equation}
(y'/y)^M \varphi(y,y)=\varphi(y,y')\,\rho(y,y')\,,
\end{equation}  
and rearanging the summations
one gets from \eqref{QQ-prod3},
\begin{equation}\label{QQ-prod4}
\displaystyle
\overline{{\bf
  Q}}(q^{-1}\,
y)\,{{\bf
  Q}}(
y')
=\displaystyle
Z_0^{-1}\,\varphi(y,y')\,
{\mathbf T}^{(+)}(y,y')
\end{equation}
as stated in \eqref{factor} in the main text.




























