On the solutions to the multi-parametric Yang-Baxter equations

A unified approach is applied in the consideration of the multi-parametric (colored) Yang-Baxter equations (YBE) and the usual YBE with two-parametric R-matrices, relying on the existence of the arbitrary functions in the general solutions. The colored YBE are considered with the R-matrices defined on two and three dimensional states. We present an exhaustive study and the overall solutions for the YBE with $4 \times 4$ colored R-matrices. The established classification includes new multi-parametric free fermionic solutions. In the context of the given approach there are obtained the colored solutions to the YBE with $9 \times 9$ R-matrices having 15 non-zero elements.


Introduction
The Yang-Baxter equations (YBE), being formulated in the early works [1,2] (originated in [3]), have crucial role in the theory of the integrable models in low dimensional statistical physics, quantum field theory and are the non-dividable parts of the theory of the quantum groups [4,5,6,7,8,9]. Investigations of the solutions to YBE are still actual and are involved in newer fields of the theoretical and mathematical physics.
The usual quantum YBE are formulated as the system of the equations Here R ij (u, v) is n 2 × n 2 matrix, u and v are called spectral parameters. The matrix R ij (u, v) is considered as an operator acting on the tensor product of two n-dimensional states V i ⊗ V j . The simplest and well studied solutions to YBE are the 4 × 4-matrices with eight non-zero entries. The symmetric solution is presented in the paper [5], which is the solution corresponding to the 2d classical statistical eight-vertex model (or the XY Z Heisenberg model). This solution has difference property -R ij (u, v) = R ij (u − v), so the R-matrix is actually one-parametric, and the corresponding YBE can be presented as R ij (u)R ik (u + w)R jk (w) = R jk (w)R ik (u + w)R ij (u). (1.2) In the context of the 1 + 1 quantum field theory [7,8], where the R-matrix plays the role of the scattering matrix of two particles, the YBE just ensures the factorization of the many-particle scattering matrices into the products of two-particle scattering matrices. Here the spectral parameters are simply the rapidities of the relativistic particles, and the difference property reflects the relativistic invariance of the system. Besides of the rapidity the scattering matrix can be depended also from the extra characteristics of the particles -"colors". In the papers [13,14,19] the multi-parametric YBE are studied with colored R-matrices R ij (u; p, q) (the color parameters p and q are attached to the the i-th and j-th spaces respectively) R ij (u; p, q)R ik (u + w; p, r)R jk (w; q, r) = R jk (w; q, r)R ik (u + w; p, r)R ij (u; p, q). (1. 3) The corresponding solutions have the so called free-fermionic property [11], R 00 00 R 11 11 + R 01 01 R 10 10 − R 10 01 R 01 10 − R 11 00 R 00 11 = 0, (1.4) which means that the corresponding 1d quantum theories can be presented by means of the scalar free fermionic chains. In the paper [17] there is discussed a two-parametric solution, which is just the solution to the YBE with R-matrices having two pairs of rapidities, so that for each pair the R-matrix has difference property -R(u, v; u ′ , v ′ ) ≡ R(u − v, u ′ − v ′ ). Then YBE has the following form R ij (u, u ′ )R ik (u + w, u ′ + w ′ )R jk (w, w ′ ) = R jk (w, w ′ )R ik (u + w, u ′ + w ′ )R ij (u, u ′ ). (1.5) In the recent paper [18] there are presented multi-parametric solutions, which being freefermionic, do not correspond to the mentioned cases. The new 4 × 4 solutions [18] (see the subsection 4.4 therein) are formulated as two-parametric solutions, which in general case have no difference property and contain arbitrary functions.
The one of the purposes of this paper is to fill the gap, which exists in the study of the multiparametric systems of YBE and in the hierarchy of the free-fermionic solutions for the case n = 2. Then we show that the YBE put the restriction on the number of the possible "colors", which can be "visible" in the solutions for the given n. For analyzing the set of YBE we propose and apply a straightforward way, formulated precisely in section 3. We classify here the all types of the multi-parametric solutions for the matrices of the general eight-vertex kind, analyzing YBE with the 4 × 4-matrices. Meanwhile the whole approach could be valid also for the cases with higher dimensional matrices, and particularly, we consider the YBE with n = 3 too and derive a family of multi-parametric R-matrices. This set of the solutions entirely describes the situation when the R-matrices have 15 non-zero elements. The extensions of the solutions to the cases with the matrices having more non-zero elements will be presented further. In particular cases the obtained solutions correspond to the known three-parametric colored solutions, the eight-vertex model's solution and the other already observed YBE solutions both for n = 2, 3 [5,13,14,19,18,22,23,24,25,26,27,28].
The paper is organized in the following way. In the Section 2 the YBE is presented in general multi-parametric formulation. In the subsections 2.1 and 2.2 the particular cases are discussed, corresponding to the situations b i = 0 and d i = 0. In the subsection 2.3 the general case is discussed. There the generalized multi-parametric version of the eight-vertex model's matrix is obtained, and we see that beside of these solutions the all remaining solutions have free-fermionic property. We are presenting the general free-fermionic solutions as matrices with two arbitrary functions (it corresponds to the four-parametric solutions) and two arbitrary constants. When one of the constants vanishes, the solution coincides with the known colored solutions [13,14]. The sl q (2)-invariant solution in [18] at q = i for the case c 1,2 = 0 (the eigenvalues of the quadratic Casimir operator of the representations spaces on which R-matrix acts), belongs to the family of the four-parametric solutions when there are two arbitrary constants and one arbitrary function, while the other function is parameterized by exponential function in a specific way imposed by the condition of the algebra invariance. In the Section 3 we present all the results and conclusions followed from the performed analysis. Therein a table is presented with the classification of the all principal types of the 4 × 4 solutions. In the context of the given approach in the next section we have solved multi-parametric YBE for general colored 9 × 9 R-matrices with 15 non-zero elements. The first part of the Appendix is devoted to the detailed discussion of the analysis of YBE for a particular case typical for the general solutions with arbitrary functions, and in the second part the main different types of the 4 × 4-solutions are presented in apparent matrix formulations.

General multiparametric solutions to YBE with 4 × 4 Rmatrices
The general form of the quantum Yang-Baxter equations with multi-parameters can be presented as Here under the spectral parameters u, v, w, written in the "bold" shrift, we mean the all possible set of the parameters, which the R matrix can acquire - As the matrix acts on the tensor product of two vector spaces, it means that the sets of the spectral parameters are attached to the corresponding vector spaces. Let us explore the matrices having the property of "particle number" conservation by mod(2) (Z 2 grading symmetry). In matrix representation it means R kr ij = 0 if i + j + k + r = 0 (mod 2). If the indexes i, j... take only two values 0, 1, then R has the following 4 × 4 matrix form which is similar to the R-matrix of the eight-vertex model. We use the commonly adopted notations (2.2) for the matrix elements. We are considering the complete and motivated cases, with a i = 0, c i = 0 anď 3) whereŘ = P R, P is a permutation operator changing the positions of the states, and I is the unit operator. The simplest equations followed from YBE (2.1) are The general solutions to them we can parameterize by means of an arbitrary function g(u) and an arbitrary constant d 0 As it is known there is a possibility to redefine the matrix elements of the matrix R ij by the following transformations [9] R p i p j induced from the following change of the vector basis e n i,j (n i,j = 0, 1, p i,j = 0, 1) of the space V i,j -e n i,j → f n i,j e n i,j . The transformations affect only the elements R 10 01 , R 01 10 , R 00 11 , R 11 00 . Taking the functions in this way f 0 (u) , we can make the elements c(d) i , i = 1, 2 equal. So, we can consider c 1 = c 2 and d 1 = d 2 afterwards.
We shall explore all the possible cases in detail. We can set c 1 = c 2 = 1 taking into account the normalization freedom of the R-matrix. The set of the independent equations of YBE is brought in the Appendix, (A.2-A.7) (the conditions (A.8) ensue from the YBE for each discussed case, as it will be shown).
In this subsection we are demonstrating a plain case with the conditions a 1 = a 2 , b i = 0, i = 1, 2, for which the extra colored parameters are actually absent, and can be introduced only by taking into account the transformation freedom (2.8). From the analysis of the whole set of the equations the following general solution follows The first two equations are just the relations expressing the full functions f i (u, w), f = a, d via the elementary functions f i (u). We see, that the functions a i (u) and d(u) are defined so, that a i (u) = a i (u, u 0 ) and d(u) = d(u, u 0 ), if at the point u 0 we have a(u 0 ) = 1, d(u 0 ) = 0. Note, that u = w corresponds to the normalization point a(u, u) = 1, d(u, u) = 0. Further we can fix u 0 = 0 = {0, 0, ...} symbolically. The last equation (2.11) is the consistency condition of the YBE. This relation is true also for the functions a i (u, w) and d i (u, w). Any R-matrix having the elements (2.9, 2.10) with arbitrary function d(u) and with a(u) = ± d(u) 2 + 1 − d(u)d 0 is YBE solution. As here we have only one arbitrary function, then the dependence from the multi-spectral parameters u one can encode in the argument of that function. We see, that the consistency condition implies trigonometric parameterization. If to impose the constraints a i (u, w) = a i (u − w) and d i (u, w) = d i (u − w), then the relations (2.9, 2.10) become functional equations on the corresponding functions and the solution is the following (2.12) which is the solution [17] describing the XYZ-model with the coupling parameters J x = −J y , Note, that the restriction that the dependence of the matrix elements from the spectral parameters must be in difference form R ij (u − w) fixes the functions a(u), d(u), meanwhile in general these two functions are connected by one relation (2.11). Hence one of the functions, say a(u), is arbitrary and we can choose any parameterization for it. However by means of definite parameterization (2.12) any two-parametric solution R ij (a(u), a(w)) with the conditions (2.9-2.11) can be brought to the actually one-parametric form R ij (u − w).

d
We have a i (u, u) = 1, b i (u, u) = 0 from the relation (2.3). Let us again take a symbolic fixed point u = 0, and set a i (u) = a i (u, 0), b i (u) = b i (u, 0). The YB equations give the following expressions for the complete functions f i (u, w) by means of the elementary functions The all relations between the functions a i (u), b i (u), coming from the YBE, can be expressed by the following constraints of familiar type The relations (2.14, 2.15) in general case are not valid for the full matrix elements f i (u, w). Here the parameters ∆, k are constants, and the second equation (2.15) must take place when ∆ = 0 (otherwise it is not necessary condition). Indeed, one could also immediately this analyzing the YBE in the common way [5], considering the equations (2.1) as homogeneous linear equations in respect of the functions f i (u, v) (or either f i (u, w) or f i (v, w)). Then the relations on the matrix elements are arisen from the consistency conditions of the homogeneous equations (formulated as vanishing of the determinants of the matrices composed by the corresponding coefficient functions). Taking now the YBE as homogeneous linear equations in respect of the functions depending, say on the parameters (v, w), and writing down the consistency conditions of the equations as 16) we see that when the free fermionic condition takes place (the second row in (2.16) vanishes), then it is not necessity for the vanishing of the first row.
In general the solutions contain two arbitrary functions, say a 2 (u) and b 1 (u) (the other two functions, a 1 (u) and b 2 (u), can be obtained from the constraints (2.14, 2.15)) and two arbitrary parameters ∆ and k (2.17) As it was noted, we can encode dependence from two arbitrary functions into the arguments -two sets of independent spectral parameters (u = {u, p}, w = {w, q}), f i (u, w) = f i (u, w; p, q). The choice of the appropriate parameterization, in which the dependence from some variables would have difference property, say R(u, w; p, q) = R(u − w; p, q), brings to differential equations. Taking then 0 = {0, 0}, we shall have Combining the relations (2.13, 2.14, 2.15), and expanding the functions f i (u, w; p, 0) near the point w = 0, we shall come to f ′′ i (u; p) ≈ f i (u, ; p) (the second differential is taken over the variable u). It means that the parameterization must be taken by trigonometric functions over the variables u, w. So, the functions will have the form f i (u; p) = g i (p) sinh [u + u i ], with appropriate chosen functions g i (p) and constants u i , depending on k, ∆. From another hand the desired parameterization one could obtain in this way. We can see that the functions g(u) and g(w), with g(u) 2 , the solution can be written as The arguments in the previous discussion correspond to u = {u, p}, w = {w, q}, and φ(u) = u, g(u) = p, φ(w) = w, g(w) = q. If to take into account also the parameters followed from the transformations (2.8), then after the following notations Of course, also another choices are possible for parameterization, say φ(u) = u, g(u) = e u (p+u), f 0 (u) = e αu t, and so on. The dependence from the arbitrary constant c 0 can be eliminated by redefinition g(u) → √ c 0 g(u) (p → √ c 0 p, q → √ c 0 q). From the ordinary XXZ the parameterization (2.20) differs by transformations of the basis states similar to the discussed one (2.8), if to redefine the matrix elements, as it was in the paper [17] for the solution (5.25) it corresponds to the transformations (2.8) after interchanging the indexes 0 and 1).
The generalized relations similar to the eq.s (2.14) and (2.15) for the complete functions f i (u, w) are the following ones Here the first relation in (2.14) is enough for the functions in (2.13) with arbitrary a i (u), b i (u), connected by one constraint, to be solutions to YBE. The matrix elements of R(u, w), taken into account the consistency conditions, can be written as As in the previous discussion, we can interpret the solutions in two equivalent ways. We can say, that we have three independent and arbitrary functions (e.g. a 1 , a 2 , b 1 ) in the solution, or we can say, that we have a solution with six parameters {a 1 (u), a 2 (u), b 1 (u); a 1 (w), a 2 (w), b 1 (w)} (three pairs of independent parameters). In general the two-parametric solution with n arbitrary functions is equivalent to the 2n-parametric solution. Writing the functions in terms of the composite parameters u, w we unify and extend two interpretations. The demand that the dependence from the spectral parameters to be in difference form (with one spectral parameter) brings the solutions to the known cases, which are the XX-model's When u 0 = π/2 the second case corresponds to the first one. In general we can take as example with arbitrary multiparametric functions f (u), g(u), φ(u) and arbitrary constant φ 0 . Note, that in the unique constraint (2.25) there is no any constant, and hence φ 0 is not a relevant constant and can be eliminated by appropriate reparameterizations, and the same is valid also for the constant u 0 introduced below in (2.28).
With the parameterization (2.27) the independent functions are expressed by the arbitrary functions φ(u), f (u), g(u). When f (u) = 1/g(u) theŘ-matrix in terms of the argumentfunction φ(u) acquires difference property. If to set the composite arguments consisting of three parameters, u = {u, p,p} and w = {w, q,q}, and fix the functions in this way we can write the matrix elements of R(u, w; p, q;p,q) as Taking into account also the parameters followed from the automorphism (2.8), and the notations (2.19), we can introduce new parameters t, s and write as well (below now u = {u, p,p, t} and w = {w, q,q, s}) the expressions for the c i elements of the matrix R(u, w; p, q;p,q; t/s).

d
The case with the choice d(u, w) = 0 brings to the following relations for the full functions via the elementary ones (belowī = i + 1 mod (2)) As previously, we set Analyzing the YBE with the functions (2.30) we obtain several expressions for d(u, w) by means of the elementary functions f i (u), f i (w). Consistency conditions for these expressions to be equal one to another constitute the relations between the elementary functions.
Before proceeding farther, let us reveal the nature of the constants included in the solutions by means of the derivatives of the elementary functions at the normalization point 0. For the multiparametric functions f i (w) we can suggest that derivation d dw is taken along a path We can see, that the constant d 0 in (2.9-2.11) can be expressed as . The constants ∆ and k in (2.14, 2.15) are equal to .
It means that the classification of the solutions by means of the constant parameters of the solutions can be formulated in the language of the derivatives taken at the normalization point. Particularly, we have seen for the case d i (u, w) = 0, that the solutions with a ′ 1 (0) + a ′ 2 (0) = 0 describe free-fermionic models. We shall ascertain that this is a general property.
So, below, we shall use derivatives at the point 0 for describing the constants. The consistency conditions, mentioned at the beginning of this subsection, contain the following general relation valid for each case Let us discuss different cases separately. The next part of this section is organized as follows.
In the first subsection we discuss the case when the constant in the constraint (2.31) (a ′ 1 (0) − a ′ 2 (0))/d ′ (0) vanishes. Here we discuss in detail the two possible situations b 1 (u) = ±b 2 (u) (b ′ 1 (0) = ±b ′ 2 (0)), which after appropriate parametrizations are corresponding respectively to the solution of the XYZ model and to the two-parametric free-fermionic R-matrix. In the next subsections the solutions with non vanishing constant (a ′ 1 (0) − a ′ 2 (0))/d ′ (0) = 0 is presented. It appears that all the solutions in this case have the free fermionic property and for them (a ′ 1 (0) = −a ′ 2 (0)) in general. In the subsection 2.3.2 we separate two cases with b ′ 1 (0) = ±b ′ 2 (0). One solution here can be considered as the generalization of the free-fermionic XY model's Rmatrix (for which a ′ i (0) = 0), the second one corresponds to the known colored thee-parametric solution [13]. In the subsection 2.3.3 the general free-fermionic solutions are presented with two arbitrary functions and two non-vanishing constants.
Here we consider the situation with a ′ 1 (0) = a ′ 2 (0). The analysis of the next equations brings to the relation a 1 (u) 2 = a 2 (u) 2 . This means that we must consider the case a 1 (u) = a 2 (u), as a 1 (0) = a 2 (0) = 1 (2.3). One can verify, that the situation a 1 (u) = −a 2 (u) brings to the rather trivial solutions.
Let us omit by now the indexes i = 1, 2. So, in fact, this case is equivalent to the Baxter's discussion of YBE for the XYZ-model. The functions f (u) obey to the following constraints: Similar relations take place for the functions f (u, w) too, with the same as it follows from the relations (2.30). Thus we have the following expressions for the functions f i (u, w) with one arbitrary function a(u) and two constants, as the other functions d(u) and b(u) can be expressed by means of them using the relations (2.32, 2.33).
We can set here the usual XYZ elliptic parameterization, placing the arbitrariness of the solution on the argument function φ(u), and the constants φ 0 and k (elliptic module), and writing So, here actually we have only one-parameteric solution, as the functions f i (u, w) acquire difference property by this parameterization f (φ(u), φ(w)) = f (φ(u) − φ(w)).
Here we found that a ′ (0) = 0 and This is a free-fermionic condition, as it immediately implies The relation (2.36) appears to be enough for theŘ(u, w)-matrix to satisfy the Yang-Baxter equations. So, there are two arbitrary functions in this solution. One can choose a parameterization, which will bring to the familiar solution (see, e.g. [17]).
, then the full functions become Setting we shall obtain the two-parametric solution [17] For the case a ′ 1 (0) = a ′ 2 (0), and hence a 1 (u) = a 2 (u), besides of the the relation (2.31) also the free-fermionic property for the functions ( 2.40) This relation takes place for the functions f i (u, w) too. Analyzing the next part of the YBE we obtain the following equation nevertheless. One could obtain this constraint also directly from the relation (2.40) expanding it around the point 0, i.e. it is the peculiarity of the free-fermionic property providing that the condition (2.3) takes place. Thus the relation (2.41) in some sense is the analog of the relation (2.16), although all the solutions here are free-fermionic. Anyway, let us at first consider the situation with the constraint a 1 (u)b 1 (u) = a 2 (u)b 2 (u), which is the feature of the usual eightvertex model and gives the solution with one arbitrary function and two constants.
. The functions f i (u) satisfy the following relations The first relation is just the free-fermionic condition.
The parameters x f and f 0 are arbitrary constants and are chosen so, that x f = a ′ 1 (0)/d ′ 1 (0) and f 0 = d ′ 1 (0)/b ′ 1 (0). Note that for this case b ′ 1 (0) = b ′ 2 (0). The solutions can be parameterized in the following form (below (2.46) We have introduced the function f x (u) as an independent function which has the property f x (0) = 0. Of course it is possible to choose different parameterizations, taking for independent function, as example the function d(u). It could seem that here (taking into account the existence of the factors with square roots in (2.46)) it is reasonable to impose the elliptic functions, but note, that there are two arbitrary constants in the square root factors and the elliptic parameterization is significant only when one of them vanishes. Indeed, when we try to parameterize so, that the difference property to be for a pair of the parameters (u, w) from the set {u, w}, f i (u, w; ...) = f i (u − w; ...), then we come to the simple elliptic parameterization with the particular case a ′ 1 (0) = −a ′ 2 (0) = 0, i.e. x 0 = 0, and a 1 (u) = a 2 (u), b 1 (u) = b 2 (u), which is the solution corresponding to the XY-model, containing one arbitrary constant.
One can prove (after some tangled calculations), that the functions above satisfy to the whole set of YBE, although they can be obtained only by considering the particular cases of them, when one of the composite arguments has been taken to be 0. This is a natural result, as the point 0 is chosen arbitrarily, we could take any point w 0 for defining the elementary functions f i (u) = f i (u, w 0 ), and then the primary values f i (w 0 ) for them will be fixed from the normalization condition (2.3). Also the free-fermionic condition can be proved for the entire functions f i (u, w), f = a, b, c, d. Another parameterization would be given in the Appendix, with a short description of a receipt how to solve YBE. So, for this particular case we have the solutions with one arbitrary function and two arbitrary constants, and in general this is a two-parametric solution. • The solution discussed above is the special case of the general solutions with the property a ′ 1 (0) = −a ′ 2 (0). And in general from the overall analysis of the YBE we come to the next important constraint This relation together with (2.31) and (2.40) is enough for f i (u) to satisfy the whole set of YBE. The corresponding general solutions will be discussed in the subsection 2.3.3. Here we should like to investigate separately the case with the vanishing constant (b ′ 1 (0) + b ′ 2 (0))/d ′ (0) = 0 (i.e. a 1 (u)b 2 (u) + a 2 (u)b 1 (u) = 0) (2.47) which as we shall see, corresponds to the elliptic three-parametric solution [24].
There are three relations on the five functions a 1 (u), a 2 (u), b 1 (u), b 2 (u), d(u), and hence the solutions contain two arbitrary functions (four-parametric solution) and one arbitrary constant ( , which comes from (2.31)). From the study of the relations it appears that this is just the "colored" solution presented in the works [13,14], provided that we require difference property for one of the pairs of the parameters. Let us describe this solution.
Fixing the arbitrary functions as d(u) and f z (u) = b 2 (u)/a 2 (u) we find for the remaining functions (below two auxiliary functions g x (u) and d f (u) are introduced and instead of d(u) the function d f (u) could be considered as elementary function) The matrix elements now look like as , a 2 (u, w) = a 1 (w, u). (2.53) The function d(u, w) depends only on the elementary functions d f (u) and d f (w). If to demand that d(u, w) = d(u − w) then we shall have the following differential equation on the function the solution of which is the Jacobi's elliptic function sn[u, k] with k 2 = −x 2 f . This corresponds to the colored parameterization presented in [14], if to take the composite parameters in the matrix R(u, w) to be u = {u, p}, w = {v, q}, and and to fix the next arbitrary functions as f z (u) = p, f z (w) = q. Thus one can check that the whole matrix has difference property in respect of the variables u and w -R(u, w) = R(u, w, p, q) = R(u − w, p, q). The full matrix elements are given in the Appendix. The case with x f = 0 has been discussed in the subsection 2.3.1.
The function d(u, w) looks similar to the expression of the previous case, only the constant parameter

(2.58)
This means that the constraint d(u, w) = d(u − w) leads to the Jacobi's elliptic function sn[u, k] with k 2 =x 2 0 − x 2 f . If to request also that the remained functions have the difference property, we shall come to a particular case, with the following differential relation for the next arbitrary , then we shall have the non-homogeneous free fermionic solution of the paper [17] (corresponding to the 2d Ising model, the matrix form of which is presented in the Appendix, at the particular limit, [10]), with the requirement b 1 = b 2 and the parameters k 2 =x 2 0 − x 2 f and dn[u 0 , k] = x f /x 0 . At the general case the solution can be represented by the following matrix elements (the expression of the element d(u, w) is written above) Here we have two arbitrary and independent functions d(u) and f g (u) and two constantx 0 and x f , the other functions are expressed by them through the relations In general we can write the R-matrix as The sl q (2)-invariant solution in the paper [18] containing one arbitrary function h(u, ε), two parameters ε 1,2 (the characteristics of the cyclic irreps of the sl q (2) algebra at q = i) and two arbitrary constants (see the subsection 4.4 in the mentioned work) belongs to the discussed here case (2.66) -the functions f g (u), d(u) are dependent from the functions h(u, ε) and e ε in a specific nonlinear way (it can be obtained just by comparing the matrix elements R kr ij ), which is given in the Appendix.

Main Conclusions and Outlook
We can summarize the results and conclusions following from the analysis performed in the Section 2 regarding the 4 × 4 solutions of the YBE in the following interrelated points (statements).
• Yang-Baxter equations define the number of the arbitrary colored parameters (independent functions), on which the R-matrices can be dependent. For the XY Z-type 4×4 multi-parametric R-matrices, the maximal number of the relevant parameters is six (or eight, if we take into account also the parameters connected with the automorphisms (2.8)).
free-fermionic solutions with one arbitrary function with two arbitrary functions and two constants: a 1 (u) = a 2 (u) free-fermionic solutions with two arbitrary functions and two constants R(p, q, u, w) the complete set of the independent equations on the functions f 1 (u) the complete set of the independent x 0 the special cases of the solutions with difference property particular cases R(u, w; p, q,p,q) (2.28) Table 1: Classification of the R 22 -solutions • As it was said in the Introduction the scattering matrices of the relativistic particles in 1 + 1 theory depends on the difference of the rapidities and the YB equations, which in this case constitute the factorization behavior of the multi-particle scattering matrices [7,8,9], actually depend only on two spectral parameters. We have shown, that the YBE has the mentioned "relativity" property, even if/when the matrices R(u, w) = R(u−w), in the following sense. When we obtain the solutions to YBE (2.1) for the fixed values of one of the composite parameters, e.g. w = 0, with defined initial conditions followed from (2.3), the solution is valid also for arbitrary w.
• Starting from the usual YBE with two parametric R(u, w) we can obtain all the multiparametric solutions, as when in the solution we get the arbitrary, non-fixed function, we can regard it as another arbitrary parameter.
• The difference property of the dependence of the R-matrices from the spectral parameters is natural property in the scattering theory, when the matrix plays the role of the scattering matrix of two particles, and the spectral parameters are simply the rapidities (u, w). When the particles have additional extra symmetries (are "colored"), the R-matrix in respect of the parameters, which describe their extra characteristics ("colored" parameters p, q) may have no difference property, as it was in the well known solutions [13,14,19]. The performed analysis shows, that the number of the independent extra characteristics, which can be shown in the Rmatrix, is restricted (maximally two kind of "colored" parameters, if we do not take into account the automorphism considered in the first part of the paper), conditioned by the number of the arbitrary functions in the YBE solutions. In the cases, when we have no difference property for any pair of the parameters, but there are some special values of the constants, at which for a pair of the parameters the difference property can be recovered, we may consider again such parameters as "rapidities".
• Now let us turn to the constants existing in the solutions. We see, that there can be utmost two arbitrary constants, which are arisen naturally in the solutions. In the discussions they are noted by k, ∆ or x 0 , f 0 (we are neglecting the constants, which can be introduced via the arbitrary functions, e.g. as in (2.28), and which do not play any role in the classification of the solutions). When all the eight matrix elements are non-zero, then the classification of the solutions leads to definite relations on the elementary functions and constants formulated by means of the derivatives given at the normalization point. A principle one is a ′ 1 (0) = ±a ′ 2 (0). For the "plus" sign, this relation implies a 1 (u) = a 2 (u) and b 1 (u) = ±b 2 (u). For the "minus" sign we have a 1 (u) = a 2 (u) and free-fermionic conditions. The algebraic equations on the functions f i (u), f = a, b, d, i = 1, 2 for both cases are three and are dependent on the constants (b ′ 1 (0) + b ′ 2 (0))/d ′ (0) and a ′ 1 (0)/d ′ (0). In the table Tab.1 we collect together the obtained results on the principal types of the 4 × 4 solutions. The cases with eight and six non-zero entries are considered separately, as just by taking d(u) → 0 we could not recover all the multi-parametric solutions.
Note, that one can choose parameterizations so, that f ′ i (0) = 0 for all the functions f i , but the constants do not vanish, as they must not depend on the parameterization, and then one can use instead the higher derivatives at that point. From this observation it could seem, that the classification in terms of the derivatives of the functions at the given point is dependent from the parameterization, but as the relations, which we use for classification for the solutions are correct relations for any and each parameterization, this classification is justified and must be understood in that sense.

General statements
In usual, solving the system of YBE we take into account, that in respect of the functions f i with any of the arguments (u, w), (u, v) or (v, w), the YB equations form system of homogeneous linear equations. And it immediately gives consistency conditions for the solutions following from the requirement of vanishing of the determinants of the matrices formed by means of the appropriate coefficients, as it was for the symmetric eight-vertex model [5]. For large matrices, with considerable amount of different matrix elements, such determinants in general can be not so easily factorizable and hence not so much informative about the nature of the solutions. We shall adhere to another way for solving YBE -a rather simple algebraic way, described below, and implemented in the analysis for the case n = 2 performed in the previous sections. Let us do the notations YB equations now read as K ijk [u, v, w] = 0. We shall consider the matrices with the property (2.3). Let us recall once again, that the parameters u, w are synthesized, joint parameters with meaning of the sets of the parameters {u}, {w}, connected with the corresponding states.
The matrix elements at the point (u, 0), where 0 symbolically denotes a normalization point, we regard as elementary functions f i (u, 0) = f i (u). Solving YBE, we consider the following scheme.
• We express the two-composite parametric functions in terms of the elementary functions by solving the equations ( * ) • We find the all set of the equations put on the elementary functions by the YBE ( * ) and also (K ijk [u, 0, w] = 0, K ijk [0, u, w] = 0), after inserting there the complete functions The number of the independent equations G k = 0 defines the number of the arbitrary elementary functions which can be in solutions, or, in another words, the number of the independent parameters -possible colors. As the point 0 is chosen arbitrarily, the solutions must satisfy also to the complete equations K ijk [u, v, w]. For the all obtained in this paper solutions it is verified by direct calculations. The sets G k for the 4 × 4-matrices are presented apparently in the Table  1 for each case.

Colored 9 × 9 solutions to YBE
The simplest known 9×9 solution is the sl 2 -invariant solution of the form R(u, w) = P +(u−w)I, P is the permutation operator and I is the unit matrix [22,23,24]. We shall consider here the simplest colored generalization of this solution, with the non-vanishing elements R ij ij and R ji ij , where the indexes i, j... take the values 1, 2, 3.
We study the cases, with a i = 0, c i = 0 and with the condition (2.3). The simplest equations followed from YBE (2.1) are for i = 1, 2, 3. Then we have the following equations with similar equations when instead of c i the functionsc i stand. The next simple equations are The general solutions to the above equations we can parameterize by means of the elementary functions f i (u) ≡ f i (u, 0), according to the previous discussions, and the constants in the solutions we shall parameterize by the derivatives of the elementary functions taken at the point u = 0. The relations which follow from the above equations (4.4-4.9), supposing that R ji ij = 0, when i = j, can be presented as Via the the transformations (2.8), now with n i,j = 1, 2, 3, p i,j = 1, 2, 3, which affect only the elements R ji ij , i = j, we can make two pairs of the symmetric matrix elements identical one to another. Taking the functions in this way f 2 (u) c 3 (u) the following elements become equal -c i (u) = c i (u), i = 1, 3. For i = 2 after this choice that equality does not take place in general, as the fraction f 2 (u) f 3 (u) is fixed by the previous relations, and by such basis transformations we could makec 2 = c 2 , only if c 2 ≈ c 1 c 3 . So, we shall consider in the followingc 1 = c 1 ,c 3 = c 3 ,c 2 = c 2 . Of course there is an arbitrariness in our choice, as we could take to be equal the other two pairs (i = 2, 3 or i = 1, 2) as well and consider the situations with c 1 = c 1 orc 3 = c 3 .
As any solution to YBE is defined up to a multiplicative function, we can fix also a 1 (u, w) = 1. So, initially we have 12 elementary functions f i (u). Then the analysis of the whole set of YBE by means of the scheme defined in the Section 3 shows, that after definition of the functions F i ( * * ), the remaining equations give nine constraints G i ( * * * ) on the functions f i (u), including the three relations (4.11). We shall not present the large expressions of the constraints G i and will present immediately the solutions. So, there are three independent functions in the solutions. We take for distinctness the functions c i (u) as arbitrary ones. The result, after doing the following successive redefinitions non-vanishing elements, constrained with definite symmetry relations on the matrix elements. The N = 3 case of the general U (1) N −1 -symmetric R-matrices, discussed in [28], can coincide with the matrix (4.12) after some symmetry transformations, as one of the authors of [28] has kindly checked 2 . If to require that difference property takes place then the matrix (4.12) would be equivalent to the Perk-Schultz solution for the three-state case [25,26] (we are thankful to prof. J. Perk for drawing our attention to this point), where the function f (u) is fixed by Then setting x f = e η , and with appropriate choice of the remaining functions and constants, the solution in (21) of the encyclopedia article [26] can be recovered. The case x f = 1 corresponds to the rational limit.
Here the general solution has actually one important arbitrary constant - , one rapidity parameter, the role of which plays the function f (u), and two type of colors, described by functions p = c 1 (u), q = c 1 (w) andp = c 3 (u),q = c 3 (w). We can denote the obtained matrices as R αᾱ 33 (u, w; p, q;p,q), α,ᾱ = ±1. Also the irrelevant colors s = f 2 (u) one could take into account, if to recall the transformation freedom connected with the basis renormalization (2.8). Then the matrix elements will be changed in this way w), forming the "eight-parametric" matrix R αᾱ 33 (u, w; p, q;p,q; s/s; t/t). In the scheme of the Quantum Inverse Scattering Problem [4] the structure of the onedimensional quantum spin Hamiltonian corresponding to the R 33 -matrix can be seen presenting theŘ 33 in the operator form, preliminary doing the notations where I is the unite operator and S ± = √ 2J ± , S z = J z . J i are the normalized 3 × 3 spin-1 generators of the sl 2 -algebra, J + = 1 Then the R-matrix has the general structure (below we are omitting the arguments u, w of the functions)Ř 33 = a 1 e + ⊗ e + + a 2 e 0 ⊗ e 0 + a 3 e − ⊗ e − + (4.15) c 1 e 0 ⊗ e + +c 1 e + ⊗ e 0 + c 2 e − ⊗ e 0 +c 2 e 0 ⊗ e − + c 3 e − ⊗ e + +c 3 e + ⊗ e − .
Let us adopt a convention that we have one kind spectral parameter u = u, f (u) = u + 1 and two arbitrary functions c 1,3 (u). Then expanding by the variable w theŘ 33 (u, w)-matrix near the point u, the linear terms in the expansion will correspond to the local cell terms H i,i+1 of the respective spin-Hamiltonian with the interactions between the nearest-neighbors spins (acting on the spaces . For the simplicity we can take c i to be constantc 1 (u) = e ε 1 /2 and c 3 (u) = e ε 3 /2 , ε i are numbers. The extension for the general family of the Hamiltonian operators with arbitrary functions c 1,3 (u) would be obvious.
The operators P α,ᾱ,γ,ε i at the values γ = ±1, ε i = πn i (n i are integers) just correspond to the permutation operators -the signs α,ᾱ = ±1 and γ, e ε i are responsible to the different gradings of the spaces. Such that the ordinary non graded case α =ᾱ = 1, γ = 1, ε i = 0 corresponds to the permutation operator of the sl(2)-invariant spin-1 spaces. The values α = −1,ᾱ = 1, γ = 1, ε i = 0 correspond to the osp(1|2)-invariant case with the fundamental three dimensional spin-1/2 representation spaces with two even and one odd parity vectors,. When α = −1,ᾱ = 1, γ = 1, e ε 1 = −1, e ε 3 = 1, the corresponding operator, multiplied by a minus sign, is the permutation acting on the spaces with two odd and one even parity vectors. For the three dimensional solutions of the mentioned algebras and their quantum extensions see the papers [29,30,31,32,33], and it is interesting to mention that taking into account the correspondence between the representations of the quantum algebras osp q (1|2) and sl i √ q (2), the mentioned graded matrices can be obtained from the sl t (2)-invariant matrices at t = 1, i. The next summand withP (4.18) exists when x f = 1 and spoils the mentioned symmetries even at the points γ = 1, ε i = 0. We must menton once again, that in the course of the solution we have taken c 2 =c 2 , and this choice induces the mentioned Hamiltonian. The solutions with c 1 =c 1 or c 3 =c 3 would bring to slight changes in the structure of the asymmetric partP , conditioned by the appropriate interchanges of the coefficients in (4.15).

Summary
In this paper we have completed the list of the colored R 22 solutions by multi-parametric freefermionic solutions. Colored R 33 -matrices are obtained for the matrices with 15 non-zero entries The given approach for the solving the multi-parametric YBE, being very simple and clearcut, gives an opportunity to find colored YBE solutions for higher dimensional matrices. It is shown, that the number of the possible extra colors, on which the R-matrix can be dependent, is restricted: initially the number of the elementary functions (the matrix elements of R(u, 0) taken at the normalization point 0) is defined from the number of the non-zero elements, and then the set of the independent relations following from YBE imposed on the elementary functions determines the number of the arbitrary functions (or the number of the possible colors) existing in the solutions. The solutions for YBE ( * ) K(u, w, 0) = 0 (K(u, 0, w) = 0, K(0, u, w) = 0) are sufficient to solve the whole set of the YBE, as the taken point 0 is chosen arbitrarily. The discussed multi-parametric R-matrices being the solutions of the Yang-Baxter equations, can have usage and treatment in all areas of the theoretical and mathematical physics, where YBE are involved -integrable models, high energy physics, quantum groups, quantum information theory, statistical physics. We have seen that the four-parametric 4 × 4-solutions presented in the section 2.3.3 have their interpretation as intertwiner matrices in the theory of the quantum algebra [18] (i.e. all free-fermionic solutions can be presented as sl q (2)-invariant matrices at q = i), and the colored parameters are the characteristics of the representations of the quantum algebra. For the cases with higher dimensions also one could expect the existence of respective underlying symmetry. In the framework of the Algebraic Bethe Ansatz the multiparametric solutions bring to the richer families of integrable Hamiltonians. The cases when there are more than one pair of the spectral parameters with "difference" property (which can be interpreted as "rapidities"), as it was in (2.38, A.23), are of particular interest.
The all obtained solutions are unitary in the following sensě This relation can be proved for all the cases, using only the symmetry property of the functions f i (u, w) under the interchange of the variables u and w, and the compatibility conditions of the Yang-Baxter equations. Particularly, the unitarity relation for the R 22 -matrix for the free-fermionic solutions just means the free-fermionic property. If to take the usual unitarity condition for the matrices, then we shall have some restrictions on the solutions, but the number of the arbitrary functions does not changed, as the relation (5.21) means and bring to rather trivial solutions). As we have learnt from [17,18] such cases can include separate solutions, together with the limit cases of the general solutions which one can obtain after taking the appropriate limits (f i → 0, i = 1 or i = 2), and however all the solutions have functional dependence (i.e. existence of the definite number of arbitrary functions, or in the particular cases, elliptic, trigonometric and rational future of the solutions) similar to the case of general solutions.
to be large and rather complicated and for simplifying the evaluation process we could suggest the following reparameterizations. Let us introduce new functions f (u) and g(u), such that We can parameterize them by means of the function d(u) and two constants f 0 , x f The two parametric functions followed from (2.30) can be written in rather compact formulas Here we have introduced the following functions, which can be parameterized only by the func- D(u, w) = (d(u) − d(w)) (1 − g(u)g(w)) (1 + d(u)d(w)) (g(u) − g(w)) (A. 19) We see that the functions factorize into two parts which contain the function f (u) and the functions d(u), g(u). In the equations we can expand the relations on the series in terms of the constant f 0 and the function-factors f x1 = f g (u) 2 − 1 , as they meet only in the functions f (u), f (v), f (w). Also one can notify that the functions a 2 (u) either have been factorized from the equations or met there in the quadratic form a 2 (u) 2 = a f (u) 1+f (u) 2 , a f (u) = 1+d(u) 2 g(u) , which allows us to escape the double square roots in the equations.
As example, let us represent the investigation of the first equation (A.2). After the mentioned expansion we find, that there are only two following equations, which one has to prove and which can be done after some not so complicated calculations. The main 4 × 4 solutions to YBE obtained heretofore Here we are presenting in the apparent matrix form the already known solutions, observed and classified in [5,13,14] and two-parametric solutions in [17,18].
The eight-vertex solution with elliptic parameterization is [5] R xyz (u) =  All other solutions presented below have free-fermionic property.
The non-homogeneous solution a 1 (u) = a 2 (u), b 1 (u) = −b 2 (u) corresponds to [17] R(u; v) = The non-homogeneous one-parametric solution a 1 (u) = a 2 (u), b 1 (u) = b 2 (u) corresponds to [17] R I (u) = If u 0 = K(k), where K is the complete elliptic integral of the first kind with the module k, then the matrix R I (u) corresponds exactly to the R-matrix, which we have found in [10] for 2d Ising Model. After so called "transformations of the first degree" of the elliptic functions, this solution at u 0 = K(k) corresponds to the XY -model matrix.
The three-parametric (colored) solution in [14,13]. Defining   The functionf (ε) is related to ε andh(ε) by means of the last relation, h 0 and g 0 are arbitrary numbers [18], which can be expressed by the constant parameters f 0 and x f .