Boundary matrices for the higher spin six vertex model

In this paper we consider solutions to the reflection equation related to the higher spin stochastic six vertex model. The corresponding higher spin $R$-matrix is associated with the affine quantum algebra $U_q(\widehat{sl(2)})$. The explicit formulas for boundary $K$-matrices for spins $s=1/2,1$ are well known. We derive difference equations for the generating function of matrix elements of the $K$-matrix for any spin $s$ and solve them in terms of hypergeometric functions. As a result we derive the explicit formula for matrix elements of the $K$-matrix for arbitrary spin. In the lower- and upper- triangular cases, the $K$-matrix simplifies and reduces to simple products of $q$-Pochhammer symbols.


Introduction
In the last two decades years there was a significant growth of interest in applications of quantum integrable systems to KPZ universality [1], stochastic processes and non-equilibrium statistical mechanics [2][3][4]. The asymmetric simple exclusion process (ASEP) [5] is one of the most studied examples, both on a line and with open boundary conditions (see, for example, [6][7][8][9]). It is intimately connected to the higher spin stochastic six vertex model which has been studied on a quadrant or a semi-infinite line with simple open boundary conditions [10][11][12]. The R-matrix of the higher spin six vertex model is related to the higher weight representations of the U q ( sl(2)) algebra and the explicit formula was derived in [13].
One of the main approaches to quantum integrable systems with general open boundary conditions is Sklyanin's method [14]. This method relies on solutions of the reflection equation [14,15]. In principle, solutions of the reflection equation for higher spins can be obtained using the fusion procedure [16,17] but such formulas are not explicit and quite complicated.
In this paper we attempt to find a general explicit expression for reflection matrices for the higher spin six vertex model in a stochastic gauge. Starting with the four-parametric solution of the reflection equation for the spin s = 1/2 [18], we get explicit formulas for matrix elements of the K-matrix for any higher spin.
The paper is organized as follows. In Section 2 we review the theory of reflection equations and the construction of commuting transfer-matrices with open boundary conditions. We also slightly generalize it to include models with R-matrices lacking a difference property. In Section 3 we review the construction of the R-matrix for the higher spin six vertex model and its factorization properties. In Section 4 we derive the recurrence relations for matrix elements of the higher spin K-matrices. In Section 5 we solve these recurrence relations in special low-and upper-triangular cases. In Section 6 we introduce equations for the generating function of the matrix elements of the K-matrix in a non-degenerate case. In Section 7 we solve these equations and find a solution for the generating function in terms of the terminating balanced 4 φ 3 series. We also obtain the explicit formula for matrix elements in the form of a double sum. Finally, in Conclusion we discuss the obtained results and outline directions for further research.

Reflection equation and commuting transfer-matrices
Reflection equation [14,15] plays a fundamental role in constructing quantum integrable systems with open boundary conditions. For a given solution R 12 (x, y) of the Yang-Baxter equation R 12 (x, y)R 13 (x, z)R 23 (y, z) = R 23 (y, z)R 13 (x, z)R 12 (x, y), (2.1) the reflection equation has the following form R 12 (x, y)K 1 (x)R 21 (y,x)K 2 (y) = K 2 (y)R 12 (x,ȳ)K 1 (x)R 21 (ȳ,x). (2.2) Here we assume that the R-matrix R 12 (x, y) is a linear operator acting nontrivially in the tensor product of vector spaces V 1 ⊗ V 2 and K i (x) acts nontrivially in V i , i = 1, 2. In general, the R-matrix R 12 (x, y) does not have a difference property and the variablesx,ȳ are the "reflected" spectral parameters. For trigonometric R-matrices we have If the R-matrix is regular, i.e. R 12 (x, x) = P 12 (2.4) with P 12 being a permutation operator, then R 12 (x, y) also satisfies the unitarity condition Using the R-matrix R 12 (x, y) and the boundary matrix K(x) we can construct a double row monodromy matrix acting in the tensor product with V a being the auxiliary space. It is easy to show that as a consequence of (2.1-2.2) the double row monodromy matrix satisfies the relation Let us assume that R t 2 12 (x, y) is non-degenerate and define a linear operator [19] We define the dual reflection equation by For a given solutionK(x) of (2.10) and the monodromy matrix (2.6) we define a double row We notice that we do not require the crossing unitarity of the R-matrix, only the existence of R 12 (x, y) andK(x). The proof goes as follows (2.13) where we used (2.5), (2.9) and the fact that Tr(AB) = Tr(A t B t ) for any matrix operators A and B. Now using (2.7) and (2.10) we transform (2.13) to i.e. we showed a commutativity of two transfer-matrices (2.12).
If the R-matrix satisfies the difference property (2.3) and the crossing unitarity condition for some ρ ∈ C and a constant matrix M ∈ End(V ), then R 12 (x, y) = R 12 (x/y) exists and is given by In general, (2.15) is a stronger condition than (2.8) even for R-matrices with a difference property [19]. Notice that (2.18) can be used to construct solutions of (2.10) from any solution K(x) of (2.2) provided that R 12 (λ) is given by (2.16). There are other automorphisms between solutions of the reflection equation and its dual [14] but we will not consider them here.

The higher spin six vertex model
In this section we start with explicit formulas for the higher-spin R-matrix R I,J (λ) related to the U q ( sl(2)) algebra following [13,20]. For arbitrary complex weights I, J ∈ C we define a linear operator R I,J (λ) ∈ End(V ⊗ V ) by its action on the basis |i , i ∈ Z + of V where matrix elements R I,J (λ) i ′ ,j ′ i,j are given by the following expression [20] Here we used standard notations for q-Pochhammer symbol, q-binomial coefficients and the basic hypergeometric function 4 φ 3 (see Appendix A). The R-matrix (3.2) satisfies the Yang-Baxter equation (2.1) with three arbitrary weights I, J, K ∈ C associated with V 1,2,3 . Let us notice that an apparent singularity coming from q −2i−2j for i, j ∈ Z + in (3.2) never happens, since the sum terminates earlier either at i ≤ i + j or j ′ ≤ i + j due to a conservation law i + j = i ′ + j ′ . Therefore, the hypergeometric function in (3.2) does not require a regularization. The representation (3.2) is equivalent to (5.8) from [13] after a Sears' transformation (A.11).
From now on we will assume that weights I and J are positive integers and the R-matrix acts in the tensor product V I ⊗ V J , where V I is a finite-dimensional module with the basis |i , i = 0, . . . , I. Therefore, we will be looking at finite-dimensional solutions of the Yang-Baxter and reflection equation unless explicitly stated otherwise.
The reason for this is that the Sklyanin approach to integrable systems with open boundaries [14] relies on the crossing relation. As shown in the previous section this can be relaxed to the existence of the operator R 12 in (2.8). A sufficient condition for the operator R 12 to exist is the crossing symmetry of the R-matrix (or a weaker condition of the crossing unitarity).
To our knowledge the crossing symmetry for the R-matrix (3.2) is only known when I, J ∈ Z + . To write it down it is convenient to define a symmetric version of (3.2) In particular,R I,J (λ) with I = J = 1 is proportional to the R-matrix of the symmetric 6-vertex model.
Let us use the standard notationR 12 (λ) for (3.3) and assume that the first and the second spaces correspond to representations with weights I and J, respectively.
Then we have two relations P 12R12 (λ)P 12 =R 12 (λ), I, J ∈ Z + , (3.4) and The relation (3.4) was proved in [13]. The second relation (3.5) with I ∈ Z + can be proved by using Sears' transformation (A.11). It is done in three steps. First we apply to the RHS of (3.5) the transformation (A.11) with q → q 2 and the following choice of parameters Second, we use relation (3.4) to interchange I and J and all indices between the first and the second spaces. The result coincides with (5.8) from [13] up to a certain factor. Applying Sears' transformation again we come to the LHS of (3.5).
We can rewrite the relation (3.5) as a crossing relation where V 1 is a (I + 1) × (I + 1) matrix with matrix elements As an immediate consequence of (3.7) and the inversion relation we have a crossing unitarity relation The easiest way to see that the inversion factor in (3.9) is equal to 1 is to rewrite it in terms of the stochastic R-matrix (3.12) and use the relation (3.17) below. Following [20,21] we introduce a stochastic version of the higher-spin six-vertex model with the Using the conservation laws i + j = i ′ + j ′ for all R-matrices in (2.1) one can easily show that the twist in (3.12) does not affect the Yang-Baxter equation Let us define the following function This function was introduced in [21] for an arbitrary rank n of the U q (A n ) algebra. Here we only consider the case n = 1.
The stochastic R-matrix (3.12) admits the following factorization [20] in terms of Φ functions: (3.14) The R-matrix (3.14) satisfies the stochasticity condition [20,21] i,j The proof immediately follows from the identity which we apply twice to (3.14). Let us notice that the inversion relation for S follows from the Yang-Baxter equation and (3.15).
It is easy to see that the crossing unitarity relation (3.10) for the stochastic R-matrix S 12 (λ) takes the following form One can ask whether the relation (3.18) can be generalized to arbitrary I, J ∈ C, since the Rmatrix (3.14) is well defined in this case [21]. The answer is apparently negative. If we substitute (3.14) directly into (3.18), we get a triple sum, with two summations coming from (3.14) and a single sum coming from the summation over matrix indices in (3.18). After straightforward calculations one can see that this last sum is given again by a balanced 4 φ 3 series which terminates when either I or J is a positive integer. Then we can use Sears' transformations to prove (3.18) directly. When both I, J ∈ C, no transformation between two non-terminating 4 φ 3 series exists. A simple numerical check shows that (3.18) does not hold in this case. However, the operator (2.8) may still exist and can be used to define the dual reflection equation.
Finally we notice that there are several choices of the spectral parameter λ, when the R-matrix S I,J (λ) simplifies to a factorized form. First, it is easy to check two properties of the function Φ q (3.20) Substituting λ = q ±(J−I)/2 and λ = q (I+J)/2 we obtain The reduction (3.21) was first noticed in [22] and then generalized to the higher rank case in [21]. The weights I and J can take complex values and play the role of spectral parameters. This case corresponds to the Povolotsky model [23].
Note that S t 1 12 is no longer invertible in (3.21-3.22) and we can not define the dual reflection equation. This is similar to the TASEP model where we can still define integrable boundary conditions for TASEP as a limit from the more general ASEP model [24].

Recurrence relations for K-matrices
We are interested in finding a general solution of the reflection equation (2.2) with the R-matrix (3.12) for arbitrary higher weights.
The reflection equation (2.2) takes the form We are only interested in non-diagonal solutions of (4.1), since any diagonal K-matrix satisfying (4.4) below will be proportional to the identity matrix. It is well known that for the case of I = J = 1 the equation (4.1) admits a 4-parametric solution of 2 × 2 matrix K 1 (x) [18]. In addition, compatibility conditions of (4.1) lead to the following restrictionx A general non-diagonal solution for K 1 (x) has the following form where t + , t − , µ, ν ∈ C are arbitrary complex parameters. This parametrization of K 1 (x) naturally appears from solving equations for has two solutions µ = 1 and µ = −t − /t + . It is easy to see that these two solutions are equivalent up to a reparametrization of the remaining parameters t ± , ν. It is convenient to choose a solution Later on we will see that the K-matrix depends on the parameter µ in a simple way and one can set µ = 1 without loss of generality. Let us notice that the stochastic K-matrix (2.25) with parameters α, γ from [24] is obtained from (4.3) by a specialization Now we substitute I = 1 into the reflection equation (4.1) and obtain This is a linear system of recurrence relations for the matrix K J (y) with arbitrary J. Moreover, we can keep J as a complex parameter, since L-operators S 1,J (x) and S J,1 (x) are well defined even for J ∈ C.
In principle, its solution for integer J is known and given by the fusion procedure [16,17]. However, we are interested in finding explicit formulas for matrix elements of K J (x) or their generating function.
Whether the solution K J (x) of (4.7) will satisfy (4.1) with both I, J ∈ C is not clear. Most likely the answer is negative because the equation (4.1) contains a double sum which is terminated either by I or J. If this double sum is infinite we can not use hypergeometric identities similar to Sears' transformations. We have already seen this phenomenon with the crossing unitarity relation.
We note that a situation with the Yang-Baxter equation is different. Due to the conservation law in (3.2) internal sums in the Yang-Baxter equation are terminated by external indices. This is the reason why the solution (3.2) can be analytically continued to complex I and J [21].
To find equations for K J (y) in (4.7) we need to derive formulas for the L-operators S 1,J (x) and S J,1 (x). Specifying I = 1 and J = 1 in (3.2) and using (3.12) one can obtain after straightforward calculations where all indices i, j, i ′ , j ′ ∈ Z + (or ≤ J for integer J) and we used a notation Substituting (4.3, 4.8-4.9) into (4.7) we obtain a set of equations polynomial in x. Decoupling with respect to x we get 12 recurrence relations for matrix elements K J (y) l j . After some algebra one can see that only two of them are linearly independent (4.12) A detailed analysis of these relations for arbitrary complex J shows that any solution contains two arbitrary parameters K J (y) 0 0 and K J (y) 0 1 and we can consistently choose For any J = 1, 2, 3, . . . we impose a terminating condition (4.14) The condition (4.14) determines K J (y) 0 1 in terms of the normalization factor K J (y) 0 0 . Once (4.14) is satisfied, a simple analysis of (4.11-4.12) shows that K J (y) l J+1+j = 0, for j, l ≥ 0. In particular, for J = 1 we reproduce a solution (4.3) and for J = 2 explicit formulas for the K-matrix are given in Appendix B. For the J = 2 untwisted R-matrix (3.2) the corresponding K-matrix was first obtained in [25]. Now we will show that the condition (4.4) with µ = 1 is compatible with (4.11-4.12) for any J = 1, 2, 3, . . ..
First, we introduce two quantities Summing up (4.12) over j and taking (4.13) and (4.15) into account we can express T l in terms of S l and S l±1 . Summing up (4.11) over j and substituting T l we observe that a constant solution S l = S exists provided that µ = 1 or µ = −t − /t + independent of J. Indeed, we solved (4.11-4.12) for J = 1, 2, 3 and checked that up to an overall normalization the K-matrix is stochastic at µ = 1.

Special solutions of the reflection equation
In this section we first analyze lower-and upper-triangular solutions of the reflection equation. Let us notice that the defining relations (4.11-4.12) become trivial for diagonal K-matrices. So we assume that either t + or t − is not equal to 0.
First we set t + = 0. Then it is easy to see that a solution to (4.11-4.12) has an upper-triangular form and a simple analysis shows that The parameter t − becomes an overall factor in (4.11-4.12) and can be set to 1. The solution (5.1) is well defined even for J ∈ C. When J ∈ Z + , both indices run the values 0 ≤ j, l ≤ J. Due to the property (3.16), the matrix (5.1) at µ = 1 is stochastic where the sum terminates at j = l. Similarly, if we set t − = 0, then the solution of (4.11-4.12) has a lower-triangular form where c J is the normalization factor. All matrix elements in (5.3) become zero for j > J, J ∈ Z + due to the factor (q −2J ; q 2 ) j . If we choose i.e. the matrix (5.3) is also stochastic for any J ∈ Z + . From (5.1-5.3) we can construct upper-and lower-triangular solutions of the dual reflection equation using the mapping (2.18).
One can specify spectral parameters x and y in the reflection equation (4.1) such that all R-matrices degenerate to a single Φ function as in (3.21-3.23). This is achieved by setting Under this specialization the R-matrices in (4.1) degenerate to different limits, i.e. (3.22), (3.23) in the LHS and (3.23), (3.21) in the RHS. In particular, the R-matrix S J,I (y/x) degenerates into (3.23) which is no longer invertible. One can ask whether it is possible to start with the degenerate R-matrix (3.21) without difference property and construct solutions to the reflection equation We found that the equation (5.8) admits the following upper-triangular solution where z ∈ C and parametersx andȳ in (5.8) are not constrained by (4.2) and remain free. The reflection equation (5.8) reduces to the 4th degree relation for Φ functions where we dropped a subscript q 2 of the function Φ, α, β, γ, δ ∈ Z + , x, y, z, u, v ∈ C and all summations are finite and restricted by external indices. Surprisingly (5.10) is very hard to prove. It reduces to some transformation of double generalized hypergeometric series which we failed to identify. Moreover, this identity can be directly generalized to a higher rank n > 1 by replacing the function Φ with its U q (A (1) n ) version from [21] with all indices replaced by their n-component analogs. We checked a generalization of (5.10) for n = 1, 2, 3 and external indices ≤ 3 and leave it as a conjecture.

A non-degenerate case
In this section we study the general off-diagonal K-matrices when both parameters t ± = 0. First, we notice that any off-diagonal solution of (4.11-4.12) possesses a symmetry This can be established by substituting K J (y) l j from the LHS of (6.1) to (4.11-4.12) and showing that the resulting recurrence relations are equivalent to original ones.
Using this fact let us define a new variable t as and introduce matrices N j,l by It is easy to see from (6.1) that N j,l is symmetric Recursion relations (4.11-4.12) can be rewritten for N j,l If we impose boundary conditions N j,l = 0 for any j, l < 0, then a solution to (6.5-6.6) will be symmetric in j, l and depend on two initial conditions, say N 0,0 and N 1,0 . Moreover, N j,l will depend only on three parameters, y, t and ν and we will omit this dependence from now on.
Recently very similar equations for higher rank K-matrices were derived using a special coideal algebra of U q (A (1) n ) [26]. The authors of [26] solved the analog of (6.5-6.6) using a matrix product of local operators acting in the auxiliary q-oscillator algebra. This approach is inspired by a 3D structure of the R-matrix (3.14) and was developed by several authors [13,[27][28][29][30].
However, the equations for K-matrices in [26] depend only on the spectral parameter with no free parameters similar to t and ν above. It would be very interesting to understand whether their approach can be extended to find a matrix product solution of (6.5-6.6). Unfortunately, we failed to do this and developed an alternative approach using techniques coming from basic hypergeometric functions.
Let us introduce a generating function for matrix elements N j,l From (6.5-6.6) we can derive a system of coupled q-difference equations for F (u, v) When we derive equations for generating functions from recurrence relations, one can expect extra boundary terms in difference equations corresponding to initial conditions in (6.7), see for example equation (4.4) in [31]. However, since recurrence relations for N j,l are consistent with terminating conditions N j,l = 0 for j < 0 or l < 0, no boundary terms appear in (6.9-6.10). Expanding (6.9-6.10) in series in u and v one can check that coefficients for any solution of the form (6.7) will solve recurrence relations for N j,l . We also note that a special choice of parameters in the K-matrix (4.3) ensures that all coefficients in (6.9-6.10) factorize. This was the main reason for using such a parametrization. 7 Construction of the generating function F (u, v) Instead of solving the system (6.9-6.10) in two variables u, v we can exclude shifts in v and derive a 2nd order difference equation in u only. The result reads where [x] for x = 0 is defined in (4.10).
Our goal is to construct a solution to (7.1) which is a polynomial in u, v of the degree J for J ∈ Z + . Difference equations similar to (6.9-6.10) and (7.1) have been studied by several authors [31][32][33][34]. Their general solution is given in terms of very well-poised non-terminating 8 W 7 series. If 8 W 7 series terminates, then due to Watson's transformation formula (III.18) in [35] it can be transformed into a terminating balanced 4 φ 3 series. So for J ∈ Z + one can expect the answer in terms of 4 φ 3 series.
A realization of this program has several difficulties. First, the 2nd order difference equation for 8 W 7 (see (2.1) in [33]) has the same structure as (7.1) but with all coefficients factorized. This is not the case for (7.1). However, this can be repaired in the following way. Let us assume that a solution to (7.1) has the form where 8 W 7 (u) solves the 2nd order equation in one variable with other parameters fixed (see (2.1) in [33]). We will not give this equation here because its explicit form is not important for further discussion. We also assume that Ψ(u) satisfies the recurrence relation with ρ(u) being a rational function. We aim at finding ρ(u) which has a structure of a simple product of a ratio of linear factors. Then the function Ψ(u) can be expressed in terms of q-Pochhammer symbols. Now we substitute (7.2) into (7.1) and use the equation for 8 W 7 to exclude the term with 8 W 7 (q 2 u). It results in the relation (A(u)ρ(u) + B(u)) 8 W 7 (u) + (C(u)ρ(u)ρ(u/q 2 ) + D(u)) 8 W 7 (u/q 2 ) = 0, (7.4) where A(u), B(u), C(u), D(u) are known factors. Since 8 W 7 can not satisfy the first order difference equation, both terms in (7.4) should be identically zero. Solving these two relations with respect to ρ(u) and ρ(u/q 2 ) we get two compatibility conditions for the function ρ(u). Further analysis shows that one can choose parameters of 8 W 7 series in such a way that ρ(u) is completely factorized and Ψ(u) is given by a product of q-Pochhammer symbols. In this way one can find two linearly independent solutions F ± (u, v) which are symmetric in u, v and solve (7.1).
The main difficulty of this general approach is that both solutions F ± (u, v) are nonterminating even for integer J. We can form a linear combination and demand that F (u, v) is a polynomial in u of the degree J for J ∈ Z + . This terminating condition can be written in terms of 3 φ 2 series and is very complicated. Substituting this back to (7.5) we should obtain the desired polynomial solution but the level of technical difficulties is so extreme that we did not succeed in finding it in any reasonable form. At least we can learn from the above calculations that a polynomial solution symmetric in u and v maybe expressible in terms of terminating very well-poised 8 W 7 series, i.e. terminating balanced 4 φ 3 series. This is indeed the case as we will see below.
Let us start with the simpler case v = 0 and construct F (u, 0). If we substitute v = 0 into (7.1), we get a difference equation for F 0 (u) = F (u, 0) In fact, it is easier to solve difference equations for coefficients N j,0 themselves. Choosing l = 0 in (6.5-6.6) and using N j,−1 = 0 we get where [x] defined in (4.10). Excluding N j,1 from (7.7) we get a three-term recurrence relation for This relation is similar to a recurrence relation for Al-Salam-Chihara polynomials [36] (see (14.8.4) in [37]) and admits a terminating solution with N j,0 = 0 for j > J in terms of 2 φ 1 series. Using a contiguous relation (A.12) it is not difficult to check that 9) where N J is the normalization factor which we will fix later from the stochasticity condition (5.5). The generating function F 0 (u) is given by To calculate F 0 (u) we first apply the transformation (A.8) to (7.9). The result reads Expanding 2 φ 1 into series in k and substituting the result into (7.10) we can calculate the sum over j using (A.6) As a result we get the following expression for F 0 (u) Applying the transformation (A.9) with q → q 2 and a = q −J y 2 νt 2 , b = u/t, c = −νy 2 q 2−J , d = q −2J u/t, e = −q 2 /t 2 (7.14) we bring F 0 (u) back to the polynomial in u Finally applying (A.10) we obtain The purpose of these calculations is to show how we arrived at (7.16). Using contiguous relations for 3 φ 2 [31] one can show that (7.16) indeed satisfies (7.6). However, it is almost impossible to guess this formula from contiguous relations for 3 φ 2 .
Having the result (7.16) one can try to generalize it to the full generating function F (u, v). We expect it to be a terminating balanced 4 φ 3 series symmetric in u and v. The only possible candidate which reduces to F 0 (u) at v = 0 is where we used (A.6). Therefore, a correct normalization of the generating function F (u, v) is given by Now we note that a pre-factor in (7.17) has a zero at u = q 2 t and only the last term with k = J in the series expansion of 4 φ 3 has a pole. Therefore, only this last term survives in the limit Comparing the result with (7.20) we obtain This is the correct stochastic normalization of the generating function (7.17). To calculate matrix elements N j,l we need to expand the generating function (7.17) back into series in u and v. We can do this using the identity ; q 2 k (q −2(J−k) ; q 2 ) s (q 2 ; q 2 ) s m={j,l} (q −2k ; q 2 ) m−s (q 2 ; q 2 ) m−s .

(7.24)
In Appendix B we give explicit formulas for N j,l with J = 1, 2. The K-matrix is now obtained from (6.3). It automatically satisfies the stochasticity condition (5.5).

Conclusion
In this paper we considered the problem of finding a general explicit solution for the reflection equation related to the higher spin representations of the stochastic six vertex model. We found a set of recurrence relations for matrix elements of the K-matrix and solved them explicitly in lower-and up-triangular cases. In the general case we expressed the generating function for matrix elements in terms of the terminating balanced 4 φ 3 series. By expanding it we obtained the expression for matrix elements of the K-matrix in the form of a double sum (7.24). It would be interesting to understand whether this formula can be rewritten in the form of a single sum, i.e. some basic hypergeometric function.
Here we did not address the problem of positivity of matrix elements of the K-matrix. However, for the case J = 1 it is well known that all elements of the R-matrix and K-matrix can be chosen in a positive regime. Since the higher spin R-and K-matrices can be built with fusion from elementary ones, we expect that positivity holds for any J.
Another interesting problem is a connection of our results with a 3D approach [13,[27][28][29][30]. In [26] a matrix product solution to the reflection equation associated with a certain coideal subalgebra of U q (A (1) n ) was constructed. Their defining equations for n = 1 are very similar to our (6.5-6.6) but do not have any free parameters except the spectral one. It is an important question whether it is possible to find a matrix product solution to (6.5-6.6) with arbitrary ν and t equivalent to (7.24). If the answer is positive, then a generalization to higher ranks should be possible 1 . We plan to address these questions in the next publication.