Solving and classifying the solutions of the Yang-Baxter equation through a differential approach. Two-state systems

The formal derivatives of the Yang-Baxter equation with respect to its spectral parameters, evaluated at some fixed point of these parameters, provide us with two systems of differential equations. The derivatives of the $R$ matrix elements, however, can be regarded as independent variables and eliminated from the systems, after which two systems of polynomial equations are obtained in place. In general, these polynomial systems have a non-zero Hilbert dimension, which means that not all elements of the $R$ matrix can be fixed through them. Nonetheless, the remaining unknowns can be found by solving a few number of simple differential equations that arise as consistency conditions of the method. The branches of the solutions can also be easily analyzed by this method, which ensures the uniqueness and generality of the solutions. In this work we considered the Yang-Baxter equation for two-state systems, up to the eight-vertex model. This differential approach allowed us to solve the Yang-Baxter equation in a systematic way and also to completely classify its regular solutions.


The Yang-Baxter equation
The Yang-Baxter equation (ybe) is one of the most important equations of contemporary mathematical-physics. It originally emerged in two different contexts of theoretical physics: in quantum field theory, the ybe appeared as a sufficient condition for the many-body scattering amplitudes to factor into the product of pairwise scattering amplitudes [1,2,3]; in statistical mechanics it represented a sufficient condition for the transfer matrix of a given statistical model to commute for different values of the spectral parameters [4,5].
Since the pioneer works in quantum integrable systems -see [6,7,8] for a historical background -, the ybe has become a cornerstone in several fields of physics and mathematics: it is most known for its fundamental role in the quantum inverse scattering method and in the algebraic Bethe Ansatz [9,10,11], although it also revealed to be important in the formulation of Hopf algebras and quantum groups [12,13,14,15,16], in knot theory [17], in quantum computation [18], in AdS-CFT correspondence [19,20] and, more recently, in gauge theory [21,22,23].
The ybe can be seen as a matrix relation defined in End (V ⊗ V ⊗ V ), where V is an nth dimensional complex vector space. In the most general case, it reads: R 12 (x, y)R 13 (x, z)R 23 (y, z) = R 23 (y, z)R 13 (x, z)R 12 (x, y), (1) where the arguments x, y and z, called spectral parameters, have values in C. The solution of the ybe is an R matrix defined in End (V ⊗ V ). The indexed matrices R ij appearing in (1) are defined in End (V ⊗ V ⊗ V ) through the formulas R 12 = R ⊗ I, R 23 = I ⊗ R, and R 13 = P 23 R 12 P 23 = P 12 R 23 P 12 , where I ∈ End (V ) is the identity matrix, P ∈ End (V ⊗ V ) is the permutator matrix (defined by the relation P (A ⊗ B) P = B ⊗ A for ∀ A, B ∈ End(V )) and P 12 = P ⊗ I, P 23 = I ⊗ P .
For physical reasons, a particular form of the ybe is usually considered: this consists in assuming that the R matrix depends only on the difference of the spectral parameters x, y, and z. In this case, the ybe assumes the simpler form: where the new spectral parameters are related to the older ones by u = x − y and v = y − z. In this work, we shall consider only the "additive" ybe (3).
For each solution of the ybe, a given integrable system can be associated. In fact, in statistical mechanics, the R matrix represents the Boltzmann weights of a given statistical model while, in quantum field theory, the R matrix is associated with factorizable scattering amplitudes between relativistic particles. From the ybe we can prove that systems described by an R matrix possess infinitely many conserved quantities in involution -the Hamiltonian being one of them -, the reason why they are called integrable [24].
We say that a given solution R(u) of the ybe (3) is regular if R(0) = P . Regular solutions of the ybe have several important properties. We list below some of them [8]: Two solutions that differ from each other only by these transformations are said to be equivalent. Regular solutions of the ybe can also have several properties or symmetries; the most common are the following [8]: -Unitarity (U) symmetry: R(u)P R(−u)P = ρ(u)I; -Permutation (P) symmetry: P R(u)P = R(u); -Transposition (T) symmetry: R(u) t = R(u); -Permutation-Transposition (PT) symmetry: P R(u)P = R(u) t ; -Crossing (C) symmetry: R(u) t1 R(−u − 2ζ) t1 = σ(u)I.
Here, t denotes transposition in End(V ⊗ V ), t 1 and t 2 denote the partial transposition in the first and second vector spaces, respectively, ζ is called the crossing parameter of the R matrix, I is the identity matrix and, finally, ρ(u) and σ(u) are two complex functions specific to each model. We highlight that we shall not impose any of these symmetries in our search for solutions of the ybe (3). Nevertheless, we indicate in Appendix B which symmetries each solution presents.

The differential Yang-Baxter equations
The ybe corresponds to a system of non-linear functional equations. Several particular solutions of the ybe are known [6,7,8]. The first solutions were found by a direct inspection of the functional equations, which are in fact very simple because the R matrix is assumed to have many symmetries. Nevertheless, there are other more advanced methods for solving the ybe: we can cite, for instance, the Baxterization of braid relations [25], the use of Lie algebras and superalgebras [26,27,28], the construction of Hopf algebras and quantum groups [13,14,15], and also techniques relying on algebraic geometry [29].
The methods mentioned above usually require that the R matrix presents one or more symmetries from the very start. From a mathematical point of view, would be desirable to develop a method that requires in principle as few as possible symmetries and, at the same time, that is powerful enough in order to find and classify the solutions of the ybe. This paper is concerned with the development and extensive use of such a method, which is based on a differential approach. To be more precise, this method consists mainly of the following: if we take the formal derivatives of the ybe (3) with respect to the spectral parameters u and v and then evaluate the derivatives at some fixed point of those variables (say at zero), then we shall get two systems of ordinary non-linear differential equations for the elements of the R matrix. The derivatives of the R matrix elements, however, can be regarded as independent variables, so that, after they are eliminated, two systems of polynomial equations are obtained in place. Thus, these polynomial systems can be analyzed -for instance, through techniques of the computational algebraic geometry [30] and eventually completely solved. It happens, however, that these polynomial systems usually have a positive Hilbert dimension, which means that the systems are satisfied even when some of the R matrix elements are still arbitrary. The remaining unknowns, nonetheless, can be found by solving a few number of differential equations that arise from the expressions for the derivatives we had eliminated before. These auxiliary differential equations, therefore, can be thought as consistency conditions of the method.
For example, if we take the formal derivative of (3) with respect to v and then evaluate the result at the point v = 0, then we shall get the equation, E := R 12 (u)D 13 (u)P 23 + R 12 (u)R 13 (u)H 23 = H 23 R 13 (u)R 12 (u) + P 23 D 13 (u)R 12 (u), (4) while, on the other hand, if we take the derivative of (3) with respect to u and then evaluate the result at u = 0, we shall get, In (4) and (5), we introduced the quantities: where P is the permutator matrix 1 so that R(0) = P and H = D(0). We highlight that H = H L P , where H L is nothing but the local Hamiltonian associated with the model -see, for instance, [6,8].
The idea of transforming a functional equation into a differential one goes back to the works of the Niels Henrik Abel, who solved several functional equations in this way [31]. Abel's method presents many advantages when compared with other methods of solving functional equations. For instance, it consists in a general method that can be applied to a huge class of functional equations; it establishes the generality and uniqueness of the solutions (which would be difficult, if not impossible, to establish in other ways) by reducing the problem to the theory of differential equations and so on -see [32] for more. Notice moreover that although Abel's method requires the solutions to be differentiable (there can be non-differentiable solutions of some functional equations), this restriction is not a problem when dealing with the ybe, as its solutions are always assumed to be differentiable because of the connection between the R matrix and the corresponding local Hamiltonian.
Concerning the theory of integrable systems, the differential method is perhaps most known in connection with boundary ybe [33,34]: This equation -which is also known as the reflection equation -is a generalization of the ybe for non-periodic boundary conditions. The fundamental unknown of the reflection equation is the K matrix -also known as the reflection matrix -, while the R matrix is assumed to be given. Notice that the K matrices figuring in (7) always depend only on a single variable, which is why the differential method transforms (7) into a system of linear algebraic equations instead of a non-linear differential system. This particularity makes the differential method as simple as powerful to study the reflection equation (7). In fact, this method was extensively used by Lima-Santos and collaborators in a series of papers, where solutions of (7) associated with non-exceptional Lie algebras and superalgebras were found and classified [35,36,37,38,39,40,41,42]. Differential methods were also employed to solve the periodic ybe (3). In fact the first solutions of the ybe were found precisely in this way [6]. What distinguishes the present approach from the previous ones is that here we regard the derivatives of the R matrix elements as independent variables, so that the system of equations (4) and (5) can be solved in an algebraic way -it is only at the end of the calculations that a few number of differential equations arise as consistency conditions. This approach allowed us to make a systematic study of the ybe; it provided a simpler and comprehensive analysis of the possible branches of the solutions which enabled us to make a complete classification of the regular R matrices for two-state quantum systems, up to the eight-vertex model 2 . Our results agree with early classifications proposed by Sogo & al. in [44] and by Khachatryan and Sedrakyan in [45]; in fact, many of the solutions derived here are equivalent to the ones presented in [44] and [45], although a few of them seems to be new -which is the case, for instance, of some solutions for six-vertex models with unusual shapes, among others.
The eight-vertex model is the most general two-state vertex model satisfying the Z 2 -symmetry [43,46], which means that the non-null elements r j1,j2 i1,i2 (u) of the R matrix must satisfy the relation i 1 − j 1 + i 2 − j 2 ≡ 0 mod 2, where the indexes i 1 , i 2 , j 1 and j 2 can assume only the values 0 or 1. Therefore, let us write the most general R matrix associated with the eight-vertex model as follows: Besides, let us denote the matrices D(u) = R ′ (u) and H = D(0) respectively by In this work, we shall look for regular solutions of the ybe (3) which means that R(0) = P , where P is the permutator matrix, explicitly given in this case by In the next section we shall present a detailed analysis of the differential ybe's (4) and (5), from which we derive the possible solutions of the ybe (3) for two-state systems, up to eight-vertex model. Other important properties of these solutions will be presented in the Appendixes.

Solutions for the four-vertex model
The most simple regular solution of the ybe occurs when the R matrix has the same shape as the permutator matrix. (8) and we are left with a four-vertex model. In this case the R matrix is,

This means that
2 A priori, the most general R matrix associated with a two-state system would be a four-by-four matrix with no zero elements -in this case, the ybe (3) would represent a set of 64 functional equations for 16 unknowns. The Hamiltonian of such a sixteenvertex model would describe a completely anisotropic Heisenberg chain in the presence of external fields and with arbitrary ionized configurations [43]. It is not known, however, if the sixteen-vertex model is integrable: its R matrix, if exists, would have no symmetry at all -not even the unitarity one. For this reason (and because the problem of solving the ybe for the sixteen-vertex model is insurmountable at the present), we shall restrict ourselves to the eight-vertex model, whose Hamiltonian is related to a Heisenberg chain in the presence of external fields but with ionized configurations occurring only in pairs [43]. 3 For sake of clarity, we shall often hide the dependence of the R matrix elements on the spectral parameter u.
so that the general R matrix of the four-vertex model is: The solution depends on four free-parameters (namely, α 1 , α 2 , γ 1 and γ 2 ). Notice, however, that one of these parameters can be removed due to the multiplicative property of the ybe, which means that the four-vertex R matrix above presents only three bare free-parameters 4 .

Solutions for the usual six-vertex model
Let us consider now the usual six-vertex model. In this case we require just that d 1 = d 2 = 0 in (8) and, consequently, the most general six-vertex R matrix becomes: We can verify that in this case the systems E and F -respectively, the equations (4) and (5) -become different each from the other, so that we get two complementary systems at our disposal. In fact, several simple relations immediately follow from these equations: for instance, subtracting E 2,5 from F 2,5 we are led to the relation, 4 In general, one or more free-parameters can be removed from the solutions thanks to the equivalence properties of regular R matrices (see Section 1). For example, by multiplying the R matrix through a given regular function, redefining the spectral parameter, performing a similarity transformation or by grouping together some of the parameters and renaming them. The remaining parameters that cannot be removed anymore from the R matrix through these equivalence transformations will be called bare free-parameters of the solutions.
On the other hand, from E 3,3 or F 3,3 it follows as well that, This means that we can write: Now we can solve the equations E 2,5 , E 4,7 , E 4,6 and E 3,5 for a ′ 1 , a ′ 2 , b ′ 1 and b ′ 2 , respectively, which provide us with, After these derivatives are eliminated, we can verify that the equations E 2,3 and E 6,7 become, respectively, On the other hand, the equation F 2,3 gives us a relation between a 2 and a 1 : Then, using (17) and (22), it follows from (21a) and (21b), assuming that b 1 = 0 and b 2 = 0, the following constraint: This means that the solutions of the ybe for the six-vertex model admit two main branches.
4.1 The first case a 2 = a 1 Let us first assume that α 2 = α 1 which, according to (22), implies a 2 = a 1 . Then we look for a solution with the following properties: Notice that, because b 2 is fixed through (17) and c 1 , c 2 are given by (19), it only remains to find a 1 and b 1 . From (20) it follows that, This system of linear differential equations can be easily solved with the initial conditions a 1 (0) = 1 and b 2 (0) = 0.
The solution is: where, At this point, all equations of the systems E and F are satisfied. Moreover, we indeed have b ′ , so that the solution is consistent with the differential method. Therefore we get the following solution: which depends on five free-parameters, namely, α 1 , β 1 , β 2 , γ 1 and γ 2 . This solution can also be written in a more convenient form by introducing the parameter η defined by the relation so that we get 5 sinh (ωη) = ǫω/ √ β 1 β 2 , where ǫ = sign (2α 1 − γ 1 − γ 2 ). Thus, in terms of this new parameter, the solution becomes: In order to count the number of bare free-parameters of this solution we can proceed as follows: first, we can simplify the R matrix by dividing all its elements by e 1 2 (γ1+γ2)u (thanks to the multiplicative property of the ybe this gives another equivalent solution). After that, we may notice that only the ratio between β 1 and β 2 appears in the solution, so that we can set β 1 /β 2 → β. In the same fashion, η always appears multiplied by ω so that we can let ωη → λ. Finally, we can redefine the spectral parameter u through ωu → u, after which γ 1 and γ 2 will appear only in the combination (γ 1 − γ 2 ) /ω, which we may call γ. This means that the solution has actually only three bare free-parameters, namely, β, λ and γ. It follows therefore that the R matrix above is equivalent to the solution named 6V(I) by Sogo & al. in [44] and that one given by equation (2.15) in the work of Khachatryan & Sedrakyan in [45]. Now, let us suppose that Taking (22) into account, this means that we are looking for a solution with the properties: In this case, the derivatives (20) become, which can be easily solved as we impose the initial conditions a 1 (0) = a 2 (0) = 1 and b 1 (0) = b 2 (0) = 0. The solution is: where, here, 5 In the whole paper, we shall use the notation ǫ = ±1.
We can verify that all the equations are satisfied when the constraint α 1 + α 2 = γ 1 + γ 2 is taken into account. Therefore, together with the expressions for c 1 = e γ1u and c 2 = e γ2u given by (19) we got a solution with five free-parameters. This solution can also be written in a simpler form after we make the transformation so that sinh (ωη) = ǫω/ √ β 1 β 2 , where ǫ = sign (α 1 − α 2 ). In fact, using (31), the solution above becomes, (37) After simplifying this solution through the equivalence transformations of the ybe (see Section 1), we can verify that it presents three bare free-parameters. This R matrix is therefore equivalent to the solution named 6V(II) in [44] and that given by equation (2.15) in [45].

Solutions for unusual six-vertex models
In the previous section, the six-vertex R matrix (16) was obtained from the most general eight-vertex R matrix (8) by zeroing the elements d 1 and d 2 . This, however, is not the only possibility for constructing six-vertex R matrices that are compatible with the regularity condition. Indeed, we might have zeroed any two of the elements b 1 , b 2 , d 1 and d 2 instead. For instance, if we vanish the elements b 1 and b 2 then we would be led to a six-vertex model whose R matrix has the following unusual shape: In the other cases, we would get the following unusual six-vertex models: In this section we shall show that the ybe (3) admits solutions for such unusual six-vertex R matrices. Let us first consider the R matrix given by (38). Our starting point here is again the analysis of equations E 3,3 and F 3,3 . As in the previous cases, these equations fix the ratio between c 2 and c 1 through the simple relation Here, however, equations E 1,1 and F 1,1 are not null, and difference of them provide us with a nice relation between d 2 and d 1 : Returning to equations E 3,3 and F 3,3 , we realize that c 2 must equal c 1 . Thanks to the multiplicative property of the ybe (see Section 1), this means that we can write 6 : Thus, we can go on by eliminating the derivatives a ′ 1 , a ′ 2 , d ′ 1 and d ′ 2 through the equations E 2,5 , E 4,7 , E 1,4 and E 2,2 , respectively, from which we get the relations: Now, from E 1,7 , it follows that after which E 2,8 becomes a quadratic equation for a 2 1 : Assuming δ 1 positive, the only solution of this equation that satisfies the initial conditions a 1 (0) = a 2 (0) = 1 and d 1 (0) = 0 is The other solutions differ from the above one by negative signs in front of the square roots, but they do not satisfy the required initial conditions for δ 1 positive 7 . After a 1 is fixed by (46), we can verify that the remaining equations imply the constraint α 2 2 = α 2 1 , so that the solutions branches into two ways. Let us consider first the branch in which α 2 = α 1 . In this case, from (44) and (46) we get that Therefore, the solution is characterized by 6 Notice that this is the same as setting γ 2 = γ 1 = 0. If we wish, we can recover the parameter γ 1 by renormalyzing the solution (e.g., by multiplying the R matrix by e γ 1 u ). With that the ybe will still be satisfied, although differential ybe's will be satisfied only if we redefine the other parameters of the solution so that the renormalized R matrix satisfies the consistency condition R ′ (0) = H. The number of bare free-parameters of the solutions are not altered by this choice, of course. 7 In general, for some complexes values of δ 1 , the initial conditions can also be satisfied when the signs in front the square roots are negative. A detailed study of these cases shows, however, that this leads to the same solution as the one presented in the text. This can be explained from the fact that δ 1 is a free-parameter of the solution, so that the negative signs can be absorbed into its definition. For this reason, in similar situations, we shall assume that the free-parameters of the solutions are positive, although the corresponding solution may be valid as well for negative, or even complex, values of these parameters.
As d 2 is already fixed by (41), it remains to find d 1 . This follows from the third equation in (43), which provides the following non-linear differential equation: This is a Riccati differential equation with constant coefficients (see [47]). To solve it, let us rewrite it in the form (we have written y in place of d 1 for convenience): where A = δ 2 , B = 2α 1 , C = δ 1 and ξ 1 , ξ 2 are the roots of the quadratic equation Ay 2 + By + C = 0. Let us first assume that ξ 2 = ξ 1 (the degenerated case ξ 2 = ξ 1 will be presented in Appendix A). In this case, the differential equation (50) can be reduced to an integral, which, by its turn, can be easily solved through partial fractions method: where c is the constant of integration. Inverting this relation, we get the desired general solution of (50): Thereby, the corresponding solution of (49) can be found by imposing the initial condition d 1 (0) = y(0) = 0, which implies the value c = log (ξ 2 /ξ 1 ) / [A(ξ 1 − ξ 2 )] for the constant of integration. Thus, after we replace back the values of A, ξ 1 and ξ 2 , we shall get that, from what follows the desired R matrix: which depends on three free-parameters (namely, α 1 , δ 1 and δ 2 ). Introducing the parameter η through the relation coth (ωη) = α 1 /ω so that sinh (ωη) = ω/ √ δ 1 δ 2 , the solution above can be rewritten as from which we see that the number of bare free-parameters is two (ω can be removed by redefining u and η). Therefore, it can be verified that this solution is equivalent to that given by equation (5.25) in [45].

The subcase a 2 = a 1
Now, let us consider the second possibility in which α 2 = −α 1 . From (44) and (46) it follows that we have, in this case, Therefore, this branch is characterized by As in the previous case, d 1 can be found through the differential equation provided by the third equation in (43). Here, however, we get the following non-linear differential equation 8 : where In the what follows, we shall show that the general solution of (58) is given in terms of Jacobian elliptic functions. To see why, notice first that (58) can be reduced to the following integral: where we have written y in place of d 1 for convenience and ξ 2 1 and ξ 2 2 are the two roots of the quadratic equation Ay 2 + By 2 + C = 0 with A, B and C given by (59). Now, assuming ξ 2 = ξ 1 (the degenerated case ξ 2 = ξ 1 will be presented in Appendix A) and making the change of variable y → ξ 1 sin φ, it follows that (60) can be rewritten as: where k = ξ 1 /ξ 2 . The integral above is the definition of the trigonometric elliptic integral of first kind of modulus k (see [47]), denoted here by F (φ|k), so that we get, where c is the constant of integration. Now, the inverse of F (φ|k) is the Jacobi amplitude function, φ(u) = am(u|k), from which we get the general solution for the function φ: 8 We highlight that the non-linear differential equation (dy/dx) 2 = Ay 4 +By 2 +C appears in several fields of mathematics and physics. In fact, in the last decade it becomes a cornerstone in some expansion methods applied to non-linear wave equations, providing in this way new elliptic solutions for several non-linear partial differential equations. For instance, this expansion method provided a countless number of solutions for the Klein-Gordon-Schrödinger [48], Kronig-Penny [49], Boussinesq [50], Korteweg-de-Vries [51], Burgers [52], Ostrovsky [53] equations and their generalizations -only to cite a few. It is interesting to notice that all the elliptic solutions of the ybe derived here follows from the solutions of this differential equation.
The solution for d 1 follows from the relation d 1 = ξ 1 sin φ, after we use the identity sin(am(u, k)) = sn(u, k) and impose the correct initial conditions. In fact, the condition d 1 (0) = 0 implies c = 0, and from A = δ 2 2 , we get that 9 , where ǫ = cosign (ξ 2 δ 2 ) . Therefore, we have determined all elements of the R matrix. After using the following well-known identities for the elliptic functions (see [47]), to simplify the elements of the R matrix, we shall get the solution, where This solution can be simplified further by introducing a parameter η through the relation iǫ δ 2 /δ 1 sn (ωη|k) = 1/ξ 1 . Thus, we can verify from the relations (59) and the identities (65) that cn (ωη|k) dn (ωη|k) = 2α 1 ξ 1 /δ 1 . Then, from the addition formula of the elliptic sinus (see [47]), it follows that (66) can be rewritten as where Λ (u, η) = 1 − k 2 sn 2 (ωu|k) sn 2 (ωη|k). Alternatively, from the identity (see [47]), we can also rewrite (66) as follows 10 : This solution is characterized by two bare free-parameters (e.g., δ 1 /δ 2 and k, as ωη is a function of k). We did not find such R matrices in the literature; we mention, however, that an elliptic R matrix with this same shape was already discussed in [45], which corresponds to a specific limit of two asymmetric eight-vertex R matrices (which are equivalent to the solutions discussed in Section 6.2 of this paper), after several elliptic transformations are performed. We were not able to verify if the R matrix (70) can be mapped to that one reported in [45] by performing some combinations of elliptic transformations and identities (e.g., Jacobi's imaginary modulus transformation, the doubleargument identities or others, see [47]). Nonetheless, we did check that these solutions have essentially the same trigonometric limits -see Appendix A.

The remaining cases
For the remaining cases of the unusual six-vertex R matrices given by (39), the functional equations are quite simple, reason for which we shall report only the final results here. They are: and 6 Solutions for the eight-vertex model The most general R matrix considered in this work belongs to the eight-vertex model and it has the following shape: (73) 10 Notice that (65) and (67) imply the remarkable identity: . This provides another way of written the R matrix (66).
Here again, we can verify that the equations E 1,1 , F 1,1 , E 3,3 and F 3,3 imply the relations Therefore, without loss, we can assume henceforward that The difference of E 2,5 with F 2,5 also provides the relation: Now we can eliminate the derivatives a ′ 1 , a ′ 2 , b ′ 1 and b ′ 2 from the equations E 2,5 , E 4,7 , E 4,6 and E 3,5 , respectively, which become, after simplification, Besides, equation F 2,3 also provides, in the same fashion as in the usual six-vertex model.
At this point, if we take the difference between E 1,4 and F 2,8 and assume that β 1 = 0, then we shall get the which evince how the solution branches.
6.1 The first case a 2 = a 1 The first branch occurs when α 2 = α 1 which, according to (78), implies as well that Now, multiplying E 2,3 by β 2 and E 6,7 by β 1 and taking the difference of them, we are led to the equation, for β 1 , δ 1 and δ 2 different from zero. It seems, therefore, that we have other three possible branches for the solutions.
Let us consider first the possibility β 2 = β 1 , which implies b 2 = b 1 . Therefore, this solution is characterized by the ratios: Now, from equation E 1,4 we can eliminate d ′ 1 : Besides, from E 1,6 and E 1,7 we can fix a 1 and d 1 : where we assumed β 1 positive (see Footnote 7) and a 2 1 = b 2 1 , since otherwise the solution would not be regular. At this point, we can verify that all the equations of the systems E and F are satisfied. It remains, however, to find b 1 . This can be achieved from the third equation in (77), which provides the following non-linear differential equation: with The differential equation (85) has the same form as (58) has same form as that for d 1 . Therefore, the desired solution for b 1 is where ξ 2 1 and ξ 2 2 are the roots of the quadratic equation Ax 2 + Bx + C = 0, with A, B and C given by (86) and we assumed ξ 2 2 = ξ 2 1 (for the degenerated case ξ 2 2 = ξ 2 1 , see Appendix A). From the identities (65), we can simplify all the elements of the R matrix, from which we obtain the following solution: where ω = ξ 2 2 δ 1 δ 2 = β 2 1 /ξ 2 1 and k = ξ 1 /ξ 2 is the modulus of the elliptic functions. Notice that ω, k, ξ 1 and ξ 2 are functions of the parameters α 1 , β 1 , δ 1 and δ 2 .
This solution (88) can be simplified further if we introduce a new parameter η through the relation sn (ωη|k) = 1/ξ 1 , so that cn (ωη|k) dn (ωη|k) = ǫα 1 /β 1 . Then, from the addition formula of the elliptic sinus (67) and using the relations (86), it follows that (88) can be rewritten as This solutions corresponds to a generalization of the eight-vertex R matrix found originally by Baxter in [4,5] and by Zamolodchikov in [54]. It contains three bare free-parameters (e.g., k, ωη and δ 1 /δ 2 ) and it is equivalent to the solution named 8V(I) in [44] and that given by equation (3.7) in [45].
Let us now to consider the possibility β 2 = −β 1 in (81), which implies b 2 = −b 1 . This means that we are looking for a solution with the following properties: As in the previous case, we can eliminate d ′ 1 from the equation E 1,4 : Then, multiplying E 2,8 by b 1 and E 3,8 by a 1 and taking the difference of them, we are led to the equation The only possibility for a 1 , b 1 and d 1 different from zero is α 1 = 0, which, of course, does not means that a 1 = 0. Now, multiplying E 1,6 by b 1 and F 1,7 by a 1 and taking again the difference, we get the equation Clearly a 2 if the solution is regular, hence, the second factor in the equation above should vanish. Solving (93) for d 1 we get, At this point, all equations of the systems E and F are satisfied. It remains however to find a 1 and b 1 . These remained unknowns can be found through the respective differential equations in (77), which become, Here we remark that a 2 1 −b 2 1 −1 = 0 implies d 1 = 0, which would lead us to a particular solution of the usual six-vertex model. Therefore, let us assume that a 2 1 − b 2 1 − 1 = 0. In this case, the system of differential equations above, with the initial conditions a 1 (0) = 1 and b 1 (0) = 0, has the solution: Now, taking the derivative of b 1 and b 2 at u = 0 we get that b ′ 1 (0) = ǫβ 1 and b ′ 2 (0) = ǫβ 2 . This means that we must set ǫ = 1 for consistency with the differential method 11 . Therefore, after simplify the previous expressions, we find the solution we were looking for: This R matrix depends on three free-parameters (β 1 , δ 1 and δ 2 ) but one of them can be removed by redefining u, so that the number of bare free-parameters is two (say, β 1 / √ δ 1 δ 2 and δ 1 /δ 2 ). It is equivalent to the solution named 8V(III) in [44] and that given by equation (3.11) in [45].
Let us back now to equation (79) and assume α 2 = α 1 . Thus, we are looking for a solution in which Before solving (79), let us eliminate d ′ 1 from E 1,4 : Then we can solve (79), say for a 1 , which provides Besides, after d ′ 1 is eliminated, we can verify that F 1,7 reduces to It seems therefore that we have two additional branches to consider, depending on whether β 2 = β 1 or β 2 = −β 1 . The second possibility, however, implies α 2 = α 1 , which lead us to the previous case discussed above (in fact, assuming β 2 = −β 1 we can verify that the difference between E 1,6 and E 8,3 gives the equation (α 1 − α 2 ) b 1 d 1 = 0). Therefore, let us consider the case β 2 = β 1 , that is, We continue in this way by solving E 1,6 for d 1 : where we assumed β 1 positive (see Footnote 7). Now, E 2,3 reduces to The first possibility α 2 = α 1 lead us again to the previous considered case in which a 2 = a 1 . Therefore, let us consider that Now all the functional equations are satisfied. It remains, however, to found b 1 . As in the previous cases, b 1 can be determined through the differential equation which is provided by the third equation in (77). In this case, however, Therefore, the desired solution of the differential equation above satisfying the initial value b 1 (0) = 0 is, where ξ 2 1 and ξ 2 2 are the roots of the quadratic equation Ax 2 + Bx + C = 0 with A, B and C given by (107) and we have assumed ξ 2 1 = ξ 2 2 (the case ξ 2 1 = ξ 2 2 will be discussed in Appendix A). Therefore, from the equations (98), (100) and (103), we can write down all the elements of the R matrix, which are: ) and we have used the identities (65) to simplify the square root in (103).
This solution can still be simplified by introducing the parameter η through the relation sn(ωη|k) = 1/ξ 1 and noticing that, remarkably, the expressions for ξ 2 1 and ξ 2 2 are quite simple in this case. In fact, we can set either or, conversely, The first case implies the identity cn (ωη|k) = ǫα 1 /β 1 where ǫ = cosign (α 1 /β 1 ), which leads to the following R matrix: where ω = β 2 ). The second case implies dn (ωη|k) = ǫα 1 /β 1 where ǫ = cosign (α 1 /β 1 ), from which we get the solution where, now, ω = √ δ 1 δ 2 and k = (β 2 To count the number of bare free-parameters of the solutions (112) and (113), we should use the relations (107) and (110) or (111). In this way, we can verify that β 1 can be removed and the solution will depend only on δ 1 /δ 2 , k and ωη, so that these solutions are characterized by three bare free-parameters. The R matrices (112) and (113) are related to the solutions originally found by Felderhof in [55] and by Bazhanov & Stroganov in [56] and they are also equivalent to the solutions reported in [44] (solution named 8V(II) in Table I) and in [45] (eqs. (3.23) and (3.19), respectively). We remark that the R matrices (112) and (113) are related to each other by the inversion of the modulus k -see Footnote 9.

Solutions for the five-vertex models
In the previous sections, we obtained the general regular solutions of the ybe for the case where the R matrix had an even number of non-null elements. There exist, however, regular R matrices with an odd number of non-null elements, which correspond to the five and seven-vertex models. In this section we shall concern with the solutions for the five-vertex models. The solutions for the seven-vertex models will be treated in the next section.
Five-vertex models can be obtained in four different ways. Two types of five-vertex models can be obtained by zeroing one of the elements b 1 or b 2 in the six-vertex R matrix (16): These five-vertex models are related to the so-called Totally Asymmetric Simple Exclusion Process (tasep) -see [57,58].
Other two types of five-vertex models can be obtained as well by zeroing the elements d 1 or d 2 in the unusual six-vertex R matrix (38): Because the main steps to solve the ybe for the five-vertex models are similar to the previous cases, in the what follows we shall only comment on the possible branches and present the final results.
This is the case corresponding to the R matrices given in (114). The solutions branch into two classes regarding on whether α 2 = α 1 or γ 1 + γ 2 = α 1 + α 2 . In the first case where α 2 = α 1 , we obtain the solutions, If we set γ 2 = γ 1 in the solutions above and make the change of variable u → log (u) /(α 1 − γ 1 ) then we shall obtain the same R matrices presented by Motegi & Sakai in [58], which are related to the tasep models.
For the second case where γ 2 = α 1 + α 2 − γ 1 , we get the solutions, These solutions present four free-parameters (α 1 , γ 1 , γ 2 and β 1 or β 1 ); however, one of them can be removed by renormalization and another one by redefining the spectral parameter. This means that the solutions above contain two bare free-parameters only. Finally, we remark that these R matrices can also be obtained from the six-vertex R matrices (30) and (37), respectively, by taking the limits β 1 → 0 or β 2 → 0 (it is necessary to eliminate η before taking the limit to avoid discontinuities).

Solutions for the seven-vertex models
Finally, let us consider the case in which the R matrix has seven non-null elements. In principle, a seven-vertex model can be obtained by zeroing one of the elements b 1 , b 2 , d 1 or d 2 in the eight-vertex R matrix (73). This would provide four possible initial shapes for the seven-vertex R matrices. It happens, however, that β 1 = 0 implies β 2 = 0 and vice-versa whenever d 1 and d 2 are both different from zero, which is a consequence of the constraint β 2 2 = β 2 1 remaining from the eight-vertex model case. Therefore, we must regard both b 1 and b 2 different from zero here, which means that there are only two possibilities for the initial shapes of the seven-vertex R matrices, namely, The ybe for these seven-vertex models can be solved in the same fashion as the usual six-vertex models. In fact, after the elements a 1 , a 2 , b 1 , b 2 , c 1 and c 2 are found by the same equations as before, other simple equations fix d 1 or d 2 and their derivatives, while the remaining equations provide some constraints and determine the branches of the solutions. Relation (23) still holds true here, which means that we have two main branches depending on whether we set α 2 = α 1 or α 1 + α 2 = γ 1 + γ 2 . In the first case, we get a solution in which a 2 = a 1 ; in the second case, we have a 2 = a 1 . We shall present these solutions in the what follows.

The first case a 2 = a 1
In this case, we have α 2 = α 1 , which implies a 2 = a 1 . The solution can be found by following the same steps presented Section 4.1. Some other simple equations determine d 1 (or d 2 ). However, differently from what happens in the sixvertex case, some of the remaining equations implies the constraint β 2 2 = β 2 1 , which means that we have two subcases to work on.

The subcase
Now, let us consider the case α 2 = α 1 and β 2 = −β 1 . Here γ 1 and γ 2 remain arbitrary and we are led to the following solutions: and The R matrices (125) and (126) depend on four free-parameters, namely, β 1 , γ 1 , γ 2 and δ 1 or δ 2 . However, β 1 can be removed by redefining the spectral parameter and, if we divide everything by e γ1u , then we can set γ =(γ 2 − γ 1 ) /β 1 as a bare free-parameter, so that the solutions have actually only two bare free-parameters. These solutions are, therefore, equivalent to the seven-vertex R matrices presented as solution 7V(III) in [44] and that given by equation (5.36) in [45]. We remark as well that the limits β 1 → 0 or β 2 → 0 of (125) and (126) reproduces the five-vertex R matrices given in (118).

The second case a 2 = a 1
Now, let us consider the case in which α 1 + α 2 = γ 1 + γ 2 . The solution can be found by following the same steps presented in Section 4.2, among with other simple equations that determine d 1 (or d 2 ). After all elements are eliminated from the systems of equations, we shall come across with following constraint: We have therefore three cases to consider. In the first case we get a solution in which b 2 = b 1 ; in the other two cases we get solutions in which b 2 = b 1 (we kindly thank the anonymous referees of this paper for drawing our attention for this possibility).

Conclusions and generalizations
In this work we developed a differential method for solving the ybe and to classify its solutions. This method allowed a systematic analysis of the functional equations arising from the ybe for two-state systems: it revealed in a simple way how the solutions branch and allowed a complete classification of their regular R matrices. In total we found thirty-one families of solutions that are associated with the four, five, six, seven and eight-vertex models -see Table  1.
In the Appendices, we also report interesting reduced solutions, which are obtained from the general ones by fixing some of the free-parameters of the R matrices. In this way, trigonometric solutions are obtained from the elliptic R matrices as well as rational solutions are derived from the trigonometric ones. The symmetries of the solutions, their geometric invariants and manifolds, the corresponding Hamiltonians and the respective classical limits (when they exist) are also discussed in the Appendices.
This work can be generalized in several ways. The most obvious generalization is the classification of the R matrices associated with three-state systems, which includes important models as the fifteen-vertex R matrix of Table 1 Classification of the solutions of the ybe for two-state systems according to the ratios of the R matrix elements. We also indicate whether or not the respective solution is of the free-fermion type (i.e., if the quantity Φ = a 1 a 2 + b 1 b 2 − c 1 c 2 − d 1 d 2 is zero or not), and we give as well the number of bare free-parameters of the solutions. In total we found thirty-one families of solutions.
Cherednik and the nineteen-vertex R matrices of Zamolodchikov-Fateev and Izergin-Korepin. Such classification is already in preparation and it will be communicated in the future. Other generalizations may include the computation of the reflection matrices (solutions of the boundary ybe) associated with the R matrices presented here, as well as the study of the statistical mechanics of the respective integrable models. It would be interesting, for instance, to present the Bethe Ansatz of these integrable models. Finally, we believe that the differential method can also be useful in the study of the non-additive ybe, the tetrahedral equation and also to find and classify the solutions of the classical ybe without make reference to the quantum ybe -this analysis could provide classical r matrices that fall outside the Belavin-Drinfel'd classification given in [60].
Acknowledgements The author kindly thanks Professor A. Lima-Santos for his comments and also the referees for their valuable remarks and suggestions. This work was fomented by Coordination for the Improvement of Higher Level Personnel (CAPES).

A Reduced Solutions
The general solutions presented in this work contain several free-parameters. When giving particular values to these parameters, particular solutions are obtained. For instance, if we set ω = 1 and γ 2 = γ 1 = 0 in the six-vertex R matrices (30) and (37), respectively, then we shall obtain the well-known simplest trigonometric R matrices of the six-vertex model [46], namely, In a similar way, the well-known R matrices of the eight-vertex model, for instance Baxter's R matrix [4,5,46], can be obtained from (89) after we set ω = 1 and δ 1 = δ 2 , so that we get ξ 1 ξ 2 = β 1 /δ 1 with k = ξ 1 /ξ 2 . Another important example was suggested by one of the referees of this paper. It corresponds to a reduced solution obtained from the eight-vertex R matrices (112) and (113) by fixing the value of η according to the expression ωη = iF √ 1 − k 2 , where F (k) denotes the complete elliptic integral of first kind whose modulus is k. In this case, we get that sn (ωη|k) → ∞, while cn (ωη|k) /sn (ωη|k) → −i and dn (ωη|k) /sn (ωη|k) → −ik. Noticing further that in this case we get β 1 → 0 but β 1 sn (ωη|k) → (ǫ/k) √ δ 1 δ 2 for (112) and β 1 sn (ωη|k) → ǫ √ δ 1 δ 2 for (113), we can verify that this value for ωη leads, respectively, to the following unusual six-vertex R matrices: and This limit was already discussed in [45] and it is related to an unusual elliptic six-vertex R matrix found by the authors of that work. The R matrices above can be compared with the elliptic R matrix (70), which has the same shape as them (see Section 5).
Other elliptic reduced solutions can be obtained by giving special values for the elliptic modulus k -or, which is the same, by considering solutions of the differential equation (dy/dx) 2 = Ay 4 + By 2 + C for particular values of the coefficients A, B and C. Regarding on these possibilities, we mention the work [53] in which the authors had presented a table with fifty-two particular solutions of this differential equation -each one of them would give place to a particular elliptic R matrix.
The most interesting way of deriving reduced solutions from the general ones is, however, to take special limits for the elliptic modulus k, so that the elliptic functions degenerate into trigonometric ones. Indeed, trigonometric R matrices can be derived from the elliptic ones by taking one of the following well-known limits of the elliptic functions -see [47]: (care should be taken in evaluating these these limits, however, because the modulus k of the elliptic functions appearing in the R matrices generally depends on complicated expressions of the solutions parameters). In the same fashion, rational R matrices can be obtained from the trigonometric ones by taking some limits, usually letting ω → 0, which degenerate the trigonometric functions into rational ones. These degenerated R matrices will be reported in the what follows (many of these degenerated solutions were already reported in [45]).

A.1 From elliptic R matrices to trigonometric ones
Let us begin our analysis with the elliptic R matrix of the unusual six-vertex model given by (70). Here the limit k → 0 is achieved by setting either δ 1 = 0 or δ 2 = 0. In any case, we get that η → ∞ and ω → 2ǫiα 1 where ǫ = sign(α 1 ); whence we obtain the following reduced solutions: and The limit k → 1, on the other hand, can be achieved by making either α 1 = 0 or α 1 = ǫi √ δ 1 δ 2 . In the first case we get the R matrix, while, in the second case, we get, We remark that the limits of (134) and (135) for k → 0 (which requires δ 1 = 0 or δ 2 = 0 in the first case and α 1 = 0 in the second) correspond to the same R matrices as given by (137) and (138), respectively, while their limits for k → 1 (which in any case is achieved as iα 1 = √ δ 1 δ 2 ) are both the same and they are given by (139) with a 1 and a 2 interchanged (this swapping between a 1 and a 2 can be explained by the symmetry of the solutions (134) and (135) regarding the inversion of the elliptic modulus k -see also Footnote 9). We can conclude, therefore, that all the three elliptic unusual six-vertex R matrices (70), (134) and (135) have essentially the same trigonometric limits. Now, let us analyze the trigonometric reductions of the eight-vertex elliptic R matrices. For the R matrix (89) in which a 2 = a 1 , the limit k → 0 can be achieved by letting either δ 1 → 0 or δ 2 → 0. However, we can verify that these limits lead to the seven-vertex R matrices (123) and (124), respectively (with iω in place of ω), so that anything new is obtained here. The limit k → 1, on the other hand, can be achieved as α 1 → ǫ β 1 − √ δ 1 δ 2 and it provides the following trigonometric eight-vertex R matrix: where ω = β 1 (β 1 − ǫα 1 ) = β 1 √ δ 1 δ 2 and tanh (ωη) = ǫω/β 1 so that sech (ωη) = ǫα 1 /β 1 . Besides, for the R matrices given by (112) and (113), in which a 2 = a 1 , we have the following: for k → 0, the limit of (112) if found as we make either δ 1 = 0 or δ 2 = 0; however, we can verify that the two seven-vertex R matrices obtained in this way are equivalent to the R matrices (128) and (129), after we replace ω by iω and u by ǫu, so nothing new is obtained here again. The limit of (113) for k → 0, on the other hand, is achieved when β 1 = ǫα 1 and it leads to the following trigonometric eight-vertex R matrix: where ω = √ δ 1 δ 2 and sin (ωη) = ω/β 1 . Finally, it remains to consider the limit k → 1 of the R matrices (112) and (113). Here, in both cases the limit k → 1 is achieved by imposing the additional constraint δ 1 δ 2 = β 2 1 − α 2 1 and it leads to the same trigonometric R matrix, namely, where ǫ = sign(α 1 ), ω = √ δ 1 δ 2 = β 2 1 − α 2 1 and tanh (ωη) = ω/β 1 .

A.2 From trigonometric R matrices to rational ones
Rational R matrices can be derived from the trigonometric ones from special limits of its free-parameters, usually through the limit ω → 0. As we shall see, all trigonometric R matrices have a non-trivial rational limit with the only exception being the four-vertex R matrix (15), whose rational limit is R(u) = P . Let us begin our analysis with the usual six-vertex R matrices (30) and (37). Setting γ 2 = γ 1 = 0 to eliminate the exponentials and then taking limit ω → 0, which implies α 1 = √ β 1 β 2 , we shall get respectively the following rational R matrices: Now, let us take the unusual six-vertex R matrix (54). Here the limit ω → 0 is achieved by imposing the additional constraint α 1 = ǫ √ δ 1 δ 2 . This provides the following R rational matrix: This solution also corresponds to the case ξ 1 = ξ 2 in the differential equation (50). For the remaining unusual six-vertex R matrices given by (66), (71) and (72), we have the following rational limits: For the eight-vertex R matrices, we have the following: the rational limit of (89) is found as we take the limit ω → 0, which implies that β 1 = ǫα 1 and either δ 1 = 0 or δ 2 = 0. This gives us two rational R matrices with seven non-null entries: Besides, the rational limit of the R matrix (89) is achieved as β 1 → 0 and δ 1 → 0 or as β 1 → 0 and δ 1 → 0; these two cases lead respectively to the same R matrices given at (145). Finally, for the R matrices (112) and (113) the rational limit is obtained as ω → 0, which implies β 1 = ǫα 1 and, either, δ 1 = 0 or δ 2 = 0. This provides us with two other rational R matrices with seven non-null entries: For the five-vertex R matrices the rational limit is found by taking simultaneously the limits α 2 → α 1 → 0 and γ 2 → γ 1 → 0. In this way, the R matrices given by (116) and (117) reduce to the following rational R matrices: while the R matrices presented in (118), (119), (120) and in (121) reduce respectively to the same R matrices given by (145). Finally, it can be verified that the seven-vertex R matrices (123) and (124) reduce to same rational R matrices given at (146), while the seven-vertex R matrices (125) and (126) reduce to the same rational R matrix given at (145). The seven-vertex R matrices (128), (129), on the other hand, reduce respectively to the R matrices given at (147), while the R matrices (130) and (131) reduce to the following R matrices:

B Symmetries and properties of the solutions
Regular solutions of the ybe can satisfy several symmetries. These symmetries can be discovered either by an algebraic approach -analyzing, for instance, the shape of the R matrix and the ratios between its elements -or through a geometric analysisstudying, in this case, the properties of the manifolds associated with the solutions. The unitarity, permutation, transposition, permutation-transposition and crossing symmetries that we commented in Section 1 depend mainly on the ratios of the R matrix elements and they can be said to be of the algebraic type. The solutions we found in this work, however, generally do not enjoy the majority of these symmetries. This, of course, is due to the fact that we assumed a priori none of these symmetries, which led us to quite asymmetric R matrices. Nonetheless, the R matrices can in general satisfy a given required symmetry if we impose the necessary restrictions on the free-parameters of the solutions. In fact, from the definitions given to these symmetries in Section 1, we can verify that the permutation symmetry requires b 2 = b 1 and c 2 = c 1 , while the transposition symmetry requires c 2 = c 1 and d 2 = d 1 . Besides, to the solution satisfy the transpositionpermutation symmetry, it is necessary that b 2 = b 1 and d 2 = d 1 and, finally, other more complicated relations are necessary for the solution to have the unitary and crossing symmetries. Notwithstanding, we highlight that all the solutions we found in this work remarkably satisfy the unitarity condition regardless of any additional constraint (i.e., in their most general form). In Table 1 we have classified the solutions into thirty-one families, according to the different ratios of the R matrix elements. The algebraic symmetries satisfied by these R matrices (within the constraints necessary for the solutions to have a given symmetry, if necessary) are presented in Table 2.
The geometric symmetries of the solutions, on the other hand, can be achieved through the analysis of their invariants, that is, the relations among the R matrix elements that are independent on the spectral parameters u and v. These invariants fix the manifolds associated with the solutions and they are a consequence of the ybe. Here again, the differential forms of the ybe given by equations (4) and (5) provide a more direct way of deriving these invariants. To see why, we can proceed as follows: first, we can eliminate all the derivatives from the systems E and F , after which they reduce to two systems of polynomial equations only. Then, methods of algebraic geometry and commutative algebra -e.g., Gröbner basis, Hilbert series, multivariate resultants, and so on -can be used to study the ideals generated by these polynomial equations [30]. In fact, we can show in this way that the Hilbert dimension corresponding to each affine variety is usually positive, which means that the system of equations are satisfied before all the unknowns are fixed. This explains why we need to solve one or more differential equations at the end of the calculations to find the remaining elements of the R matrix.
As we shall see below, the following quantity plays a key role in the analysis of the eight-vertex model invariants and its descendants [61,62]: In fact, the eight-vertex models are usually classified into two main classes according to whether Φ is zero or not: the case Φ = 0 corresponds to the so-called free-fermion models (or Felderhof-type models, after [55]), while the case Φ = 0 corresponds to the non-free-fermion models (or Baxter-type models, after [46]). If we take the derivative of (150) with respect to u and evaluate the result at u = 0 we shall get the quantity Thus, ϕ = 0 provides a necessary condition (often sufficient) for a model to be of the free-fermion type. Due to its importance, we indicate if a given vertex model is or not of the free-fermion type in Table 1.
In the sequel, we shall analyze in more details the invariants and manifolds associated with the six and eight-vertex models. A similar analysis can be done for the four, five and seven-vertex models but it will be concealed. Table 2 The unitarity (U ), permutation (P ), transposition (T ), permutation-transposition (P T ) and crossing (C) symmetries of the solutions, among with the relation defining the crossing parameter (ζ), when available. The symbol ✓ indicates that the solution has the corresponding symmetry without any constraint. If some conditions are necessary for the solution to have a given symmetry (in such a way that the shape of the R matrix remains the same), we write them down. Finally, the symbol ✗ means that the solution cannot have the given symmetry at all.

B.1 Invariants of the usual six-vertex model
Let us begin with the usual six-vertex model, in which case we must take d 1 = d 2 = 0. After the derivatives are eliminated from the systems (4) and (5), only a few equations do not vanish. Some of them fix the ratios of the R matrix elements, providing the already derived relations (see Section 4): Among the remaining equations, E 2,3 and E 6,7 are of particular importance. After simplification, they become: From these two equations, the following invariants directly follow: from which we can deduce as well that The possible manifolds associated with the six-vertex models are found after we multiply E 2,3 by β 2 and E 6,7 by β 1 and take the sum or the difference of them. In fact, this provide us with the expressions: Therefore, we can see that there are two possibilities, namely, either ϕ = α 1 + α 2 − γ 1 − γ 2 = 0 or a 1 b 1 β 2 − a 2 b 2 β 1 = 0. The first case corresponds to a free-fermion solution which was named 6V2 in Table 1. In the second case, we get a manifold characterized by which actually means that a 2 = a 1 after we use (152). Assuming ϕ = 0 we get a manifold belonging to the non-free-fermion solution that corresponds to the R matrix named 6V1 in Table 1 (notice that the factor a 2 b 2 β 1 + a 1 b 1 β 2 cannot be zero because this would imply a 2 = −a 1 in view of (152), which incompatible with the regularity condition).
B.2 Invariants of the unusual six-vertex model where b 1 = 0 and b 2 = 0 In the case of the elliptic unusual six-vertex model presented in Section 5, we have b 1 = b 2 = 0 but d 1 = 0 and d 2 = 0. After the derivatives are eliminated from the systems E and F , some equations can be used to fix the ratios between the R matrix elements, which are: In fact, we get from E 1,1 and F 1,1 for instance, the invariant, which actually means c 2 = c 1 and d 2 = (δ 2 /δ 1 ) d 2 , as stated above. Moreover, from E 1,4 , E 4,1 , E 5,8 and E 8,5 we get as well the invariants: from which we obtain the ratio between a 2 and a 1 given by the first equation in (158).
The possible manifolds associated with this unusual six-vertex model can be found from the analysis of the equations E 2,8 and F 1,7 . In fact, taking the difference between E 2,8 and E 1,7 , we shall get, after we use (158), (164) and (165), the following relation: ( Therefore, we have two possibilities: if α 2 = α 1 we get the solution given by equation (54), which is a non-free-fermion solution (solution 6U1 in Table 1). In this case, the following invariants hold: The other possibility, on the other hand, means that Φ = 0, that is, it belongs to a free-fermion solution. This case corresponds to the solution 6U2 of Table 1.
We shall not discuss the invariants of the unusual six-vertex R matrices given by (71) and (72) because the analysis is similar to the previous cases.

B.3 Invariants of the eight-vertex model
In the case of the eight-vertex model, several other equations survive after the derivatives are eliminated from the systems E and F . Some of them provide the ratios derived in Section 6, namely, As in the previous case, from E 1,1 and F 1,1 we get the invariant, and, from E 1,4 , E 4,1 , E 5,8 and E 8,5 we also get that, where we used the relation b 2 2 = b 2 1 , which always holds for regular solutions of the eight-vertex models. Here, the equations E 2,3 and E 6,7 are again of particular interest. In this case, they are: Taking the sum and the difference of them, we get, From this, we can analyze the possible manifolds associated with the eight-vertex model. Considering first the free-fermion case, where Φ = 0, we can see from (167) that there is only two possibilities here -in the same fashion as in the usual six-vertex model -, namely, we should have either ϕ = 0 or a 1 b 1 − a 2 b 2 = 0. For the first case, it is more convenient to consider the subcases β 2 = β 1 and β 2 = −β 1 separately. For β 2 = −β 1 we obtain the invariant, which means that a 2 = a 1 and b 2 = −b 1 ; therefore, this manifold corresponds to the solution named 8V2 in Table 1. For the other possibility, β 2 = β 1 , it follows that so that this manifold is associated with to the solutions 8V3 and 8V4 of Table 1. The last possibility a 1 b 1 − a 2 b 2 = 0 leads us to the invariant which actually means that b 2 = b 1 and a 2 = a 1 . In this case, however, other equations imply that ϕ = 0 as well, so that we lie into a particular case of the solutions above. Notice that for the free-fermion case, the following invariants always hold: Now, let us consider the non-free-fermion case in which Φ = 0. Here we should assume that ϕ = 0. Besides, the possibility β 2 = −β 1 leads to d 1 = 0 if the solution is regular, so that there is only one possibility here, namely, β 2 = β 1 . By its turn, this lead us to the invariant, which implies that a 2 = a 1 and b 2 = b 1 . Therefore, this manifold is associated with the solution 8V1 of Table 1. Finally, notice that the following invariants always hold for the non-free-fermion solutions: as we can see from the equations E 1,6 and E 1,7 , for instance, after we make use of the relations (163), (164) and (165). Several other invariants can be found for the eight-vertex model as well. We shall not go further on this subject, however, because these relations do not provide much more than we already have.  Table 3 Hamiltonian coefficients for each vertex model and the corresponding spin chains. Notice that the free-fermion models are characterized by Jz = 0.
or, conversely, In Table 3 we present the values of the J's coefficients for each vertex model listed in Table 1. We can see from the Table 3 that the four-vertex R matrix corresponds to the Ising spin chain, in which the interaction occurs only in the z direction. Besides, all the free-fermion solutions are associated with either the xx or xy Heisenberg spin chains. Finally, the non-free-fermion solutions are related to either the xxz or xyz Heisenberg spin chains (the xxx spin chain is related to some of the reduced solutions).

D Classical limits
We finish this paper by presenting the classical limits of the R matrices classified in Table 1. As commented in Section 1, the ybe plays a fundamental role in theory of quantum integrable systems. For classical integrable systems, an analogous role is played by the so-called classical Yang-Baxter equation: a first-order approximation of the (quantum) ybe (3) that reads: The solutions of the classical ybe are called classical r matrices. Given a solution R(u, h) of the (quantum) ybe (3) that also depends smoothly on a certain parameter h, besides the spectral parameter u, we say that this R matrix has a classical limit if the following condition holds: lim where I is the identity matrix and f (u) is any differentiable complex function. In this case, the classical r matrix is defined by the formula: and we can verify that it satisfies the classical ybe (180). In fact, (180) follows after we differentiate twice the quantum ybe (3) with respect to h, evaluate the result at h = 0 and use the property (181). The classical ybe was introduced by Sklyanin in [63] and Belavin & Drinfel'd in [60] have classified the non-degenerated solutions of the classical ybe for finite dimensional simple Lie algebras (a non degenerated solution of (180) is such that det r(u) = 0). After that, Jimbo [13] and Drinfel'd [14,15] independently introduced the concept of quantum groups as deformations of Lie algebras that allowed the construction of R matrices from the classical ones [26]. Since then, the classical ybe has appeared in connection with many important topics of theoretical physics and mathematics -see, for instance, [7] and references therein.
In this section we shall present the classical limit of the R matrices associated with two-state systems. For the cases where this limit exists, we have identified the parameter h with η and normalized 12 the R matrices so that we have f (u) = 1. We remark, however, that the majority of the R matrices reported in this work does not have a classical limit as well. This is because the condition (181) is too restrictive 13 . In fact, any R matrix whose elements b 1 or b 2 is zero cannot have a classical limit because (181) can never be satisfied. Thus, the four-vertex R matrix, all the five-vertex R matrices and the unusual six-vertex R matrices cannot have such a classical limit. Besides, the condition (181) is satisfied only if a 2 (u, 0) = a 1 (u, 0), b 2 (u, 0) = b 1 (u, 0) etc., which implies that the asymmetric six-vertex R matrix (37), all the seven-vertex R matrices except (123) and (124), and, finally, the eight-vertex R matrices (97), (112) and (113) cannot have such a classical limit. In short, the only R matrices that admit a classical limit are the six-vertex R matrix (30), the seven-vertex R matrices (123) and (124) and the eight-vertex R matrix (89). The classical limits of these R matrices will be presented below.
For the six-vertex R matrix (30), we have the following classical r matrix: where we have made β 2 = β 1 so that (181) is satisfied (see Footnote 13).
For the seven-vertex R matrices (123) and (124), we have respectively the following classical r matrices: and For the eight-vertex R matrix (89) we have the following classical r matrix: Finally, we can compare our results with the Belavin-Drinfel'd classification [60]. For two-states systems, the r matrices correspond to the sl(2) Lie algebra and for this case Belavin & Drinfel'd have shown that, up to certain equivalences, there is 12 Different normalizations lead to different r matrices that differ from each other only by the addition of a term proportional to the identity. These r matrices, however, can be regarded as equivalent thanks to the following additive property of the classical ybe: if r(u) is a solution of (180), then the matrixr(u) = r(u) + f (u)I is also a solution of it. 13 Some solutions can have a classical limit if the condition (181) is weakened. For example, it can be verified that the classical limit of the asymmetric six-vertex R matrix (37) is given by the same r matrix (183), although the condition (181) is not satisfied in this case. Also, for the six-vertex R matrices (30) and (37), the classical ybe is still satisfied without the requirement β 2 = β 1 . These cases, however, fall outside the Belavin-Drinfel'd classification so that we shall not analyze these possibilities further.
only two trigonometric solutions and only one elliptic solution, besides their reduced rational solutions. The two trigonometric solutions are associated, respectively, with the six and seven-vertex models and the corresponding r matrices are indeed equivalent to the r matrices we found here -namely, that ones given by (183), (184) and (185), respectively. The elliptic solution, by its turn, is associated with Baxter's eight-vertex model and the corresponding classical r matrix is indeed equivalent to the r matrix (186). We conclude therefore that our classification completely agrees with the Belavin-Drinfel'd analysis.

E The differential Yang-Baxter equations
We list below the non-null functional equations corresponding to the the ybe (3) for the general eight-vertex model. The functional equations for the four, five, six and seven-vertex models are contained in this system and can be obtained by zeroing the corresponding elements.
The two sets of differential equations derived from the ybe (3) are presented below. The E system corresponds to equation (4) while the F system corresponds to (5).