Yang-Baxter and the Boost: splitting the difference

In this paper we continue our classification of regular solutions of the Yang-Baxter equation using the method based on the spin chain boost operator developed in \cite{deLeeuw:2019zsi}. We provide details on how to find all non-difference form solutions and apply our method to spin chains with local Hilbert space of dimensions two, three and four. We classify all $16\times 16$ solutions which exhibit $\mathfrak{su}(2)\oplus \mathfrak{su}(2)$ symmetry, which include the one-dimensional Hubbard model and the $S$-matrix of the ${\rm AdS}_5 \times {\rm S}^5$ superstring sigma model. In all cases we find interesting novel solutions of the Yang-Baxter equation.


Introduction
The Yang-Baxter equation (YBE) appears in many different areas of physics [2]. It signals the presence of integrability which implies the existence of higher conservation laws. The equation emerges in some form in virtually every area of physics, including condensed matter, statistical physics, (quantum) field theory, string theory and even quantum information theory [3]. The Heisenberg spin chain and the Hubbard model [4] are just some of the famous integrable models and were important for our understanding of low-dimensional statistical and condensed matter systems and, similarly, over the last few years, exceptional progress has been made in understanding the AdS/CFT correspondence [5] due to the discovery of integrable structures [6]. Given the clear ubiquity of the Yang-Baxter equation throughout theoretical physics it is clear that understanding and classifying its solutions is a highly interesting and non-trivial task.
The presence of quantum integrability in a given physical model with Hilbert space C n is dictated by the existence of a solution R(u, v) ∈ End(C n ⊗ C n ), dubbed R-matrix, of the Yang-Baxter equation, which reads on C n ⊗ C n ⊗ C n and the subscripts denote which of the three spaces R acts on. The parameters u, v, w are known as spectral parameters with one associated to each of the three spaces. Once R is known one can construct the so-called transfer matrix t(u, θ) for a spin chain of length L as t(u, θ) = tr a (R aL (u, θ) . . . R a1 (u, θ)) , (1.2) which generates an infinite tower of conserved charges (Q i , i = 1, . . . , ∞) via the expansion log t(u, θ) = Q 1 (θ) + (u − θ)Q 2 (θ) + 1 2 (u − θ) 2 Q 3 (θ) + . . . . (1.3) The parameter u is an auxiliary spectral parameter, whereas the parameter θ is a physical parameter such as the rapidity of a particle in a scattering process. An immediate property of the YBE is that the charges Q r commute: [Q r (θ), Q s (θ)] = 0 (1.4) which is the cornerstone of integrability. A particularly interesting class of R-matrices are the so-called regular solutions which are those R-matrices R(u, v) which satisfy the regularity condition R 12 (u, u) = P 12 where P 12 denotes the permutation operator on the two copies of C n . The significance of such solutions is that for the corresponding integrable system, momentum is a conserved charge, or more precisely the tower of conserved charges commutes with the operator of cyclic permutations which is a prevalent feature of many integrable models such as the Hubbard model. In this case the conserved charges Q r (θ) are a sum of densities of range r, meaning each density acts on r-adjacent spin chain sites. For example, the Hamiltonian Q 2 (θ) can be written as a sum of nearest-neighbour (range 2) densities H j,j+1 (θ) as Q 2 (θ) = L j=1 H j,j+1 (θ) (1.5) and periodic boundary conditions are imposed, that is H L,L+1 (θ) = H L,1 (θ). This density is itself related to the R-matrix in a very simple way: H 12 (θ) = P 12 ∂ u R 12 (u, θ)| u→θ . (1.6) Hence, the moment one knows the R-matrix one knows the Hamiltonian and the dynamics of the system. Throughout the history of quantum integrable systems numerous different approaches have been developed for finding solutions of the Yang-Baxter equation. In the early days a very fruitful approach has been through requiring the solutions to have certain symmetries [7][8][9]. For example, if we wish for the Hamiltonian Q 2 to commute with the generators a of some Lie algebra g then one should impose that [R(u, v), a ⊗ 1 + 1 ⊗ a] = 0. More generally given some bialgebra A we require that ∆ op (a)R(u, v) = R(u, v)∆(a) where ∆ and ∆ op denote the coproduct and opposite coproduct related by conjugation on A, respectively. In many cases this is enough to completely fix R up to a small number of functions, drastically simplifying the construction, as was demonstrated in the case of AdS/CFT integrable systems [10][11][12][13], see also [14] for recent developments using this approach. Of course, this approach first requires one to know what the corresponding symmetry is and there are R-matrices which may have no such symmetry at all. Still within the realm of algebra, a more abstract approach is that of Baxterisation which initially appeared in the realm of knot theory [15] and consists of constructing solutions of the YBE as representations of certain algebras, for example Hecke algebras and Temperly-Lieb algebras. Numerous different R-matrices have been obtained in this way [16] and further advancements have also been achieved recently [17].
A more hands-on approach is to simply try and solve the Yang-Baxter equation directly. The upside to this is that in principle one can obtain all solutions in this way, but this is contrasted with the enormous difficulty of solving cubic functional equations. This approach is usually supplemented with differentiating 1 the YBE and reducing the cubic functional equations to a system of coupled partial differential equations. This approach has recently been used to provide a full classification of R-matrix of size 4 × 4 so-called 8-and-lower-vertex models [18] obeying the difference property R(u, v) = R(u − v) and to obtain certain 9 × 9 models [19] whose R-matrix satisfies the so-called ice rule but it quickly becomes unwieldy as the size of the R-matrix increases.
In this paper we follow a bottom-up approach which we have developed in a series of recent papers [1,20,21] and further develop here. In our approach, instead of starting with the R-matrix and using it to find the Hamiltonian and the corresponding dynamics, we start with the Hamiltonian and use it to obtain the R-matrix. The mechanism for carrying out this procedure hinges on the so-called boost automorphism [22,23] which is an alternative, yet equivalent, way to generate the tower of conserved charges for regular integrable models without the need to construct the transfer matrix and expand it. The boost automorphism B[Q 2 ], or simply the boost operator, is defined by (1.7) The infinite sum should be interpreted in a formal sense but what we are interested is not the boost operator itself but rather its commutator with the tower of conserved charges which is perfectly well-defined even for finite chains. In fact, it can be shown that, see Hence, by knowing just the Hamiltonian density H 12 (θ) we can construct the full tower of commuting conserved charges directly and our approach is based on exploiting this observation. Namely, instead of starting with a solution of the YBE we will start with a generic operator on C n ⊗ C n which we identify as a Hamiltonian density H 12 (θ) and construct the corresponding global charge Q 2 (θ). We will then use the boost operator to construct Q 3 (θ) by the relation (1.8). A priori there is no reason for the two constructed operators Q 2 and Q 3 to commute with each other, but if we impose this it will place a number of constraints on the entries of the density H 12 in the form of a system of ODEs. We then solve the set of constraints and show that the resulting Hamiltonian defines an integrable system, meaning it can be obtained from a solution of the YBE. In order to do this we use the so-called Sutherland equations which are obtained from the YBE and read where in each of the above Sutherland equations R ij := R ij (u, v) andṘ and R denote the derivatives of R with respect to the first and second variable respectively. The Sutherland equations constitute two sets of ODEs for the entries of the R-matrix and the boundary conditions are fixed by the requirement of regularity R(u, u) = P and the fact that the Hamiltonian density can be obtained from the R-matrix by the expansion Hence, in effect, solving the condition [Q 2 (θ), Q 3 (θ)] = 0 singles initial conditions for the Sutherland equations which have a chance to be consistent with the existence of an Rmatrix. What is remarkable is that all initial conditions obtained in this way lead to an R-matrix, at least for the cases discussed in this paper and in [1,20,21]. Let us remark that our approach is actually not the first which uses the Hamiltonian as a starting point for constructing integrable systems and R-matrices. In a series of papers [24,25] the authors present a method for determining if a given nearest-neighbour Hamiltonian system is solvable by the coordinate Bethe ansatz [26] and also produce Rmatrices for some of these systems. There is also the earlier work [27] where an iterative procedure for reconstructing the R-matrix from the Hamiltonian was developed for models where the R-matrix satisfies the difference property R(u, v) = R(u−v), as well as the work [28]. In the case of 1 + 1-dimensional integrable field theories such R-matrices correspond to S-matrices which are Poincaré invariant and include integrable systems such as the XYZ spin chain and its derivatives and Zamolodchikov's O(N ) sigma model. When one restricts to this case, our procedure described above for constructing R-matrices from Hamiltonians simplifies enormously. In particular the conserved charges Q r (θ) become independent of θ and so the set of ODEs arising from the condition [Q 2 (θ), Q 3 (θ)] = 0 reduces to a set of coupled cubic polynomial equations. This simplification was exploited in the papers [1,20] in order to find a plethora of new integrable systems with a range of interesting physical properties. In this paper, in order to demonstrate the full power of our approach we will not impose the difference property and consider the most general R-matrices. One of the main results of this paper is that the single consistency condition [Q 2 (θ), Q 3 (θ)] = 0 on the Hamiltonian is enough in order to completely determine the R-matrix even in the absence of the difference property, which goes back to a conjecture of [29]. The analysis of the most general possible R-matrices using the boost approach was initiated in [21] and here we continue that analysis. For the higher-rank case however, that is beyond 4 × 4 R-matrices, we will impose that our R-matrices have certain symmetries in order to render the calculations tractable. For 9×9 R-matrices we will impose that our R-matrices commute with the Cartan subalgebra of su(3) and for 16 × 16 R-matrices we will impose su(2) ⊕ su(2) symmetry which appears in various interesting models such as the su(4) Heisenberg XXX spin chain [8], the Hubbard model [30] and the AdS/CFT S-matrix [10] and the related Shastry R-matrix [31]. Furthermore, for 16 × 16 models we determine all possible integrable models which preserve fermion number and a Hamiltonian based on a generalisation of the usual Hubbard model. Let us point out that such restrictions are not strictly necessary to implement our approach, but due to the fact that our method produces huge numbers of integrable systems, many of which are trivially related as will be explained in the main text, providing a full classification of all possible 9×9 and 16×16 R-matrices is highly difficult and so we limit ourselves to a subset of models which are physically interesting.
This paper is organised as follows. In Section (2) we review the procedure developed in [1,20,21] and outlined here in the introduction and explain in detail how to go from a generic Hamiltonian to an integrable Hamiltonian and subsequently find a solution of the Yang-Baxter equation. In Section (3) we will apply our procedure to models with two-dimensional local Hilbert space which produces R-matrices of size 4 × 4 which were initially presented in our letter [21]. In Sections (4) and (5) we discuss models with three-and four-dimensional local Hilbert spaces with certain symmetries imposed, namely u(1)⊕u(1)⊕u(1) ⊂ su(3) and su(2)⊕su(2) respectively as well as the generalised Hubbard models mentioned earlier. Finally, we discuss some further directions for research. In Appendix A we review the construction of the charges Q r (θ) using the boost operator.
We have attached a Mathematica notebook to the arxiv submission of this paper which contains all of the R-matrices obtained in this paper as well as those in [1,20,21] together with the corresponding Hamiltonians and commands to check various properties. In Appendix B we provide some details on this notebook.

Set-up and method
In this section we will give more details on our method and discuss an explicit example to illustrate the procedure that we follow.

Method
As described in the introduction our starting point is a nearest-neighbour Hamiltonian density H 12 (θ) on C n ⊗ C n . Using the density we construct the full Hamiltonian H(θ) = Q 2 (θ) on a spin chain of length 4 and we identify sites 4 + j := j, j = 1, . . . , 4. The reason for our restriction to length 4 will be explained below. From now on, for shortness, we will sometimes omit the θ dependence on the density Hamiltonian. Using the boost operator (1.7) we construct Q 3 (θ), which in this case is given explicitly by (2.2) As was described in the introduction Q 3 (θ) can be written as a sum of range 3 densities Q j,j+1,j+2 (θ): Next, we impose the condition [Q 2 (θ), Q 3 (θ)] = 0 which is a necessary condition for the model to be integrable and solve the resulting set of ODEs for the entries of H 12 .
Then we plug the obtained density H 12 (θ) into the Sutherland equations described in the introduction solve the resulting set of differential equations for R(u, v) and fix this uniquely by using the boundary conditions Finally, we check that the Yang-Baxter equation is indeed satisfied. This procedure can be conveniently outlined in the diagram of Figure 1. At this point we would also like to highlight that in case the R-matrix is of difference form (that is satisfy the difference form property), the Hamiltonian will not depend on the spectral parameter and hence all derivative terms drop out. This means that the integrability condition simply reduces to a set of polynomial equations and that the Sutherland reduces to ordinary (non-linear) differential equations since R effectively only depends on one variable.
Finally a small comment is due on why we restrict to spin chains of length 4. Since Q 2 is a sum of range 2 densities and Q 3 is a sum of range 3 densities the non-vanishing terms in their commutator [Q 2 , Q 3 ] is a sum of densities of range 2 + 3 − 1 = 4. Hence, if we restrict to a spin chain of length 3 say, then these non-zero commutators will effectively wrap around the spin chain producing cancellations which do not happen in general. Hence for our construction we must consider spin chains of at least length 4 in order to avoid this happening. Alternatively, we can also derive the same system of equations by simply looking at the densities. In case one needs to consider the commutation relations between higher conserved charges, the length of the spin chain needs to be adjusted accordingly -if one wants to consider the commutator [Q n , Q m ] then a spin chain of length L = n + m − 1 should be considered.

Identifications
As we explained in [1,20,21], our approach of solving the YBE by imposing the condition [Q 2 , Q 3 ] = 0 leads to quite a large redundancy in solutions. Specifically, for a given  integrable Hamiltonian density (or R-matrix) there are various transformations one can do which preserve this condition and preserve regularity. Hence in what follows we will only concern ourselves with a single representative of equivalence classes of Hamiltonians and R-matrices. We now describe the transformations which lead to equivalent solutions.
Local basis transformation If R(u, v) is a solution of the Yang-Baxter equation and V an invertible matrix depending on one spectral parameter and with the same size of the R-matrix, then we can generate another solution by defining This new solution is trivially compatible with regularity and just corresponds to a change of basis on each site. On the level of the Hamiltonian it gives rise to a new integrable Hamiltonian which takes the form where everything is evaluated at θ and I is the identity matrix. In particular, we see that terms of the form A ⊗ I − I ⊗ A in the Hamiltonian can be removed by performing the basis transformation (2.8) with the matrix V (u) satisfyingV = AV which can be solved by means of a path-ordered exponential.
Reparameterization If R(u, v) is a solution, then R(g(u), g(v)) clearly is a solution of the YBE as well. This transformation affects the normalization of the Hamiltonian since by the chain rule the logarithmic derivative of R will give an extra factorġ, so that H(u) →ġH(g(u)). (2.9) Notice furthermore that this will similarly affect the derivative term in the boost operator. We are also free to reparameterize any other functions and constants in both the R-matrix and Hamiltonian. For instance the R-matrices from [32] can be obtained by using a reparameterization of the usual XXX R-matrix.
Normalization We can normalize the R-matrix in any way we want since multiplying any solution R of the YBE by an arbitrary function g is clearly allowed. On the level of the Hamiltonian this corresponds to a simple shift of the Hamiltonian where I is the identity matrix. We have imposed g(θ, θ) = 1 in order to preserve R(θ, θ) = P .
Discrete transformations It is straightforward to see that for any solution R(u, v) of the Yang-Baxter equation, P R(u, v)P, R T (u, v) and P R T (u, v)P are solutions as well.
All the above transformations are universal and hold for any integrable model. Moreover, they have a trivial effect on the spectrum, which means that they basically describe the same physical model. Additionally, there are some transformations called twists that we can use for identifications that are model dependent. Twists generically change the spectrum and more generally the physical properties of the integrable model in a nontrivial way. However, on the level of the R-matrix a twist is a simple transformation.
is a solution of the YBE provided R is. Note that much more general transformations which preserve the YBE can be obtained by combining (2.11) together with other transformations. For example, if both U and V are constant invertible matrices satisfying which can be obtained by applying (2.11) together with a similarity transformation and applying (2.11) again. We will refer to any transformation obtained by combining (2.11) with the other transformations mentioned above as a twist. Under the transformation (2.11) the Hamiltonian density H 12 transforms as and the analogue of the condition [U (u)⊗U (v), R 12 (u, v)] = 0 for the Hamiltonian density can be easily worked out to be (2.14) Alternatively this relation may be derived by plugging the twisted R-matrix (2.11) and Hamiltonian (2.13) into the Sutherland equations (1.9) and sending v → u, which is not surprising given the striking similarity between (2.14) and the Sutherland equations.
Finally, there can be other, model dependent, twists such as Drinfeld twists [33]. Moreover, the condition [U (u) ⊗ U (v), R 12 (u, v)] = 0 can also be extended to, for instance, depend on two twists U, V or the spectral dependence of the twist can be modified. We will usually only use standard twists (2.11) unless stated otherwise.

Example
As a demonstration of our method let us work out an example in full detail. From here on we will use the following notation: Hamiltonian Let us classify all regular solutions of the YBE whose Hamiltonian densities have the following form From the boost operator construction we find that the corresponding charge Q 3 has density and is quadratic in the components h i (θ) of the Hamiltonian density H. We have suppressed the θ dependence. The next step is to impose [Q 2 (θ), Q 3 (θ)] = 0 which gives the equationṡ These are solved by for some constants c 3,4 . Thus we find that if H 12 is to be obtained from an R-matrix it must have the form R-matrix We make an ansatz for our R-matrix of the following form (2.20) We will first solve the Sutherland equations using brute force before using identifications to greatly simplify the process. The Sutherland equations (1.9) give the following independent set of PDEs From this we see that Since we need to impose regularity R(u, u) = P , we find that A = 1 and B = c 3 /c 4 . Next, we derive that We are then left with three unsolved PDEṡ In order to solve these we redefine so that the last equation becomes This equation can now be most conveniently solved by substituting cylindrical coordinates, so that we find for some function φ to be determined by the remaining two differential equations. Notice that this is an overdetermined system. Plugging (2.28) then back into the remaining Sutherland equations gives the followinġ which is easily solved upon using the boundary condition that φ(u, u) = 0. Setting and combining everything we are left with the following Rmatrix after choosing the overall normalisation r 3 to correctly reproduce the Hamiltonian. Owing to the dependence on both H + and H − , this R-matrix is manifestly of non-difference form.
It is straightforward to check that R indeed satisfies the Yang-Baxter equation and that its logarithmic derivative gives the density Hamiltonian (2.19).
Using identifications The above method of finding the R-matrix can be greatly simplified if we use some identifications that relate various solutions of the Yang-Baxter equation that we discussed in the previous section. We start from (2.19) and use a local basis transformation to set h 1 = h 2 . This is achieved using the matrix V (θ) with together with the transformation law (2.8).
Next, we use reparameterization symmetry to set h 1 = h 2 = 1. Thus, it follows that all the entries of the Hamiltonian are constant and the resulting Hamiltonian density has the form (2.32) Moreover, we can use a twist and set c 3 = c 4 = c. Indeed, it is trivial to check that the twist condition (2.14) is satisfied for any constant invertible diagonal matrix U and the matrix can be used to bring the Hamiltonian density to the form after applying H 12 → U 1 H 12 U −1 1 . The Sutherland equations are now also easily solved since all the coefficients of the Hamiltonian are simply constants. As a consequence, the R-matrix is of difference form and is given by the usual XXZ solution. Putting ω 2 = c 2 − 1 we find In order to see that this solution is equivalent to the solution (2.30), let us undo the identifications that we performed to make the Hamiltonian constant. First we undo the twist and apply R 12 → U −1 2 R 12 U 1 to (2.35) and put c = √ c 3 √ c 4 so that we arrive at the R-matrix for the Hamiltonian (2.32). Next we reparameterize and finally we apply the inverse of the local basis transformation (2.31), immediately obtaining (2.30).
Difference vs. Non-difference After using all the identifications, we see that (2.30) is actually just an R-matrix of difference form in disguise. The non-difference nature of the rapidity dependence of the R-matrix only resides in local basis transformations, a rescaling and a reparameterization. These can obviously be applied to any solution of difference form to generate a non-difference form solution. In the remainder of this work we will also encounter models which are genuinely of non-difference form, but it is easy to see already at the level of the Hamiltonian if this is the case. More precisely, after solving the integrability condition [Q 2 (θ), Q 3 (θ)] = 0 our Hamiltonian will depend on a number of free functions. One will usually correspond to a shift, one can be absorbed in a reparameterization of the spectral parameter and then remains a number that can be absorbed by local basis transformations and potentially twists. The exact number of the latter will depend on the set-up. Thus in case of (2.19), we count 2 free functions h 1 , h 2 and we could have already at that point concluded that the underlying model was actually of difference form.

Two-dimensional local Hilbert space
We will now apply our approach to the classification of various different integrable systems. The first case we will consider will be the case where the local Hilbert space has dimension two and, consequently, the R-matrix is of size 4 × 4. In [21] we classified all solutions of the Yang-Baxter equation of 8-and-lower-vertex type and we review this case now for completeness. Further details can be found in [21].

8-and-lower-vertex models
8-and-lower-vertex type models are solutions of the form and, consequently, the corresponding Hamiltonian densities are of the form We will briefly recap the solutions of the Yang-Baxter equation of this type that we found. After identifying solutions, we found only four different types of integrable 4 × 4 Hamiltonians that solve the integrability condition [Q 2 (θ), Q 3 (θ)] = 0: Let us discuss the models in more detail.

6-vertex A
where c 3,4 are constants. The Hamiltonian is actually equivalent to that of the XXZ spin chain. Indeed, by applying a local basis transformation, twist, reparameterization and normalization we can bring the Hamiltonian density to the form which is precisely the Hamiltonian density (2.34), and so its R-matrix is given by (2.35).
Notice that this solution also contains the most general diagonal Hamiltonian since only the off-diagonal elements h 3,4 are restricted by the integrability condition. We now introduce a reparameterization of the spectral parameter which kills the non-derivative term in the Riccati equation and removes the explicit dependence on h 3 . It is then straightforward to solve our system of differential equations to find where again H i (x, y) = x y h i . It is instructive to write the R-matrix as We see that h 5 gives rise to the non-difference nature of this solution. In particular, when h 5 is constant the R-matrix reduces to an R-matrix of XXZ type. It is easy to show that it satisfies the Yang-Baxter equation and the correct boundary conditions. This model can be mapped by a twist into the solution A of the pure colored Yang-Baxter equation considered in [34].

8-vertex A
In the case h 6 = 0, the integrability constraint gives that where c i are constants. The resulting Hamiltonian is that of the XYZ spin chain [8,18] under our identifications.

8-vertex B
In the case when h 6 = 0, we find the following differential equationṡ 14) We use a local basis transformation to set h 2 = 0 and then these equations are solved by .
By using a local basis transformation we can set c 8 = c 7 and after applying further identifications the remaining functions can be brought to the following form where η is some free function. This further results in h 7 = h 8 = 2c 7 := k, which all together imply that r 5 = r 6 = 1 and r 7 = r 8 for the R-matrix. The remaining functions are easily determined from the Sutherland equations and we find where sn, cn, dn are the usual Jacobi elliptic functions with modulus k 2 and where η ± = η(u)±η(v) 2 and all the Jacobi elliptic functions depend on the difference u − v, i.e. sn = sn(u − v, k 2 ). This solution indeed satisfies the Yang-Baxter equation and has the correct boundary conditions. Moreover, it is easy to see that in the case where η is constant, it becomes of difference form and reduces to the well-known solution found in [18,35]. Furthermore, in the limit k → ∞ the R-matrix reduces to that of the AdS 2 integrable system [12,21]. We would also like to remark that in [21] we presented the Hamiltonian with a different parameterisation than used here as well as two solutions of the differential equations and hence the R-matrix -the two solutions are actually related by a twist thanks to the symmetry [R 12 , σ z ⊗ σ z ] = 0.
Off-diagonal model As can be seen from (3.14)-(3.16), the cases where h 5 = 0 and h 3 = −h 4 need special attention due to possible singularities. In particular it is easy to see that by setting h 5 = 0 it follows that the Hamiltonian is constant unless h 3 = −h 4 . And, indeed, in our final expression the limit h 5 = 0 corresponds to setting η(x) = π/2. However, the case h 3 = −h 4 warrants special attention. In this case, the entries of the Hamiltonian are We see that the Hamiltonian for this model only has off-diagonal entries. It can be shown that it is possible to recover this model, starting from the Hamiltonian of 8-vertex B.
Since the procedure is highly non-trivial, we explain the steps of this identification. In order to recover (3.27) we followed the following steps: In order to makeH 8V B verify the integrability condition [Q 2 , Q 3 ] = 0, we fixed one entry of the twist a → s 1 √ c 8 b, with s 1 = ±1, ±i and we had to impose a constraint on the entries of the Hamiltonian with α 3 some θ-dependent function. Notice that this twist is non-standard as is does not satisfy (2.14).
2. We apply a diagonal local basis transformation V (θ). In particular by using (2.8), we first fixV V −1 to eliminate the elements in the (2,2) and (3,3) positions of the Hamiltonian. Then by solving the differential equations, we fixed the matrix V (θ).

3.
We get an off-diagonal Hamiltonian density and we checked that the sum of the elements at position 2,3 and 3,2 is zero if s 1 (defined in step 1) is ±i. Moreover the ratio between elements in 1,4 and 4,1 is constant.
In this way we have recovered model (3.27) from H 8V B . Since the twist that we used is non-standard, it is unclear how to easily lift it to the level of the R-matrix. Nevertheless, it is easy to solve the Sutherland equations for this model directly and we obtain

(3.30)
We see that it is of quasi-difference form, meaning all of the dependence on the spectral parameters is of the form H 3 (u) − H 3 (v) and H 7 (u) − H 7 (v).

Hermitian solutions
We postpone the classification of all regular 4 × 4 solutions of the Yang-Baxter equation to future work due to the complexity of the equations and their solutions. However, there is one interesting physical subcase which we can fully classify. We can classify 4 × 4 solutions that give a Hermitian spin chain Hamiltonian. Hence, let's assume we have a Hamiltonian density of the form where h (r) and h (i) are the real and imaginary parts of the entries of the Hamiltonian. All the functions are now real-valued and imposing hermiticity leaves us with 16 independent real functions. Hence in solving the integrability condition we can set both real and imaginary parts to 0. Moreover, we can discard all solutions that have complex numbers in them. This greatly simplifies our computation and, remarkably, we find that all the solutions of this type can be brought into 8 vertex form under our identifications. Note that this does not mean that the corresponding 8 vertex models are Hermitian. There are non-Hermitian 8 vertex models which, after a non-diagonal basis transformations, become Hermitian but are then no longer of 8 vertex type.

Three-dimensional local Hilbert space
Next we apply our method to R-matrices of size 9 × 9 corresponding to local Hilbert spaces of dimension 3. In the literature there are many examples of such models including [9,19,24,25,27,28,[36][37][38][39]. We consider models whose R-matrix and Hamiltonian density commute with the Cartan subalgebra of su(3) which are usually referred to as 15-vertex models. These models are a special case of models satisfying the so-called ice rule [27] which states that for an R-matrix with components R µα νβ in the standard basis we have the constraint We complete the classification of fifteen-vertex models developed in [19,37]. As a result of the Cartan symmetry we consider a Hamiltonian density of the form where H := H(θ) and h ij := h ij (θ). The corresponding R-matrix is of the form where R := R(u, v) and r ij := r ij (u, v).
Applying the procedure described in the previous sections we obtain ten independent models of non-difference form. We had four models for which all h ij in (4.2) and all r ij in (4.3) were nonzero. However, after applying the identifications presented in section 2.2 we learned that these models were actually difference form models disguised by twists, local basis transformations and reparameterizations and corresponded exactly to the models obtained in [19] and [37] 2 .
The six remaining models are fundamentally of non-difference form, i.e. they cannot be brought to difference form by applying the transformations described in section 2.2, and to our knowledge are new models. An interesting fact is that for the non-difference form 9 × 9 cases with Cartan symmetry su(3), none of the Hamiltonians are Hermitian. This is in contrast with the 4 × 4 case where one can construct models which commute with the Cartan subalgebra of su (2) and are still of non-difference form and Hermitian, for example the model (3.1) under a special choice of its free functions.
A curious fact, is that all the cases that are fundamentally of non-difference form for the fifteen vertex models, were the ones coming from singular cases, i.e. cases where we started with some extra zeros in the Hamiltonian (4.2) from the beginning of the procedure.
Below we present the new models of non-difference form. They are divided in two classes depending on whether p(u, v) in equation (4.7) is equal to the permutation operator P (Class 2) or not (Class 1).

Hamiltonian densities
The nonzero matrix elements of the Hamiltonian density (4.2) for each of the six new models are presented in the tables below.

Class 1
The matrix elements of the density Hamiltonian for models in Class 1 are given in Table  1 Model

R-matrices
By exploiting identifications we found that we could bring all of the corresponding Rmatrices to a form closely resembling the XXX spin chain which we remind the reader is of the form R(u) = uI + P (4.6) where I denotes the identity operator and P is the permutation operator. We found that we can always write the R-matrices for this section in the form where f (u, v) is some function, d(u, v) is a diagonal 9 × 9 matrix and p(u, v) is a matrix with the same entries being non-zero as the usual permutation operator in the standard basis. In fact for models 5 and 6 this operator is exactly the permutation operator.

Class 1
The models in Class 1 have in common the fact that they all possess the same p i given by where E 86 is a matrix with 1 in position (8, 6) and 0 everywhere else. All of the diagonal matrices for these models can be written in the following form where E kk denotes the matrix with 1 in position (k, k) and 0 everywhere else.

Class 2
This class has two models and they both have p equal to permutation, i.e., p = P. (4.10) Model 5 It is described by and where Model 6 It is given by where G ± = G(u) − G(v), I ± = I(u) ± I(v) and I(u) is defined in (4.4) and (4.5). Notice that model 6 actually defines two separate models due to the choice ± of signs. They are independent due to the nontrivial dependence of I(u) on j(u) (see (4.4) and (4.5)) and cannot be mapped to each other using any of the transformations in section 2.2.
In order to check the YBE for model 6 in Mathematica one needs to be careful with the choice of branch in j(u). The best approach is to substitute R(u, v) in the YBE without specifying j(u), then simplify as most as possible and only then substitute j(u) 2 as in equation (4.5). By doing in this way, one never actually needs to choose a branch and the YBE is immediately satisfied. This was already implemented in our Mathematica notebook.

Four-dimensional local Hilbert space
We now apply our method to the case where the local Hilbert space is of dimension 4. In order to have a manageable set-up we restrict ourselves to models which have su(2) ⊕ su(2) symmetry. This class of solutions of the YBE contains important R-matrices which correspond to the Hubbard model and AdS/CFT. Moreover, in [20] we also discovered some interesting new models whose R-matrix was of difference form. Following it, we see that there are two classes of models with su(2)⊕su(2) symmetry. The first class are models where both su(2) transform in a four-dimensional representation. These models are of so(4) type via the isomorphism so(4) ∼ su(2) ⊕ su (2). The second class are models where the su(2) are represented two-dimensionally. The Hubbard model falls into this category. Finally, we discuss some generalisations of the Hubbard model obtained by taking the Hamiltonian of the free Hubbard model and including the most general possible potential term which preserves fermion number, which allows to interpret the model as electrons moving on a one-dimensional lattice or conduction band. The different form analogue of this setting was discussed in [20].

so(4) type models
The most general Hamiltonian density underlying this symmetry takes the form where P is the permutation operator, I is the 16 × 16 identity matrix and K = E ij ⊗ E ij with (E ij ) αβ = δ i,α δ j,β . Summation over repeated indices is assumed and i, j, k, l = 1, ..., 4. We found only one possible integrable model of non-difference form with the following Hamiltonian This is the non-difference form model corresponding to the usual so(4) spin chain, but here the constant coefficients become functions of the spectral parameter. The R-matrix corresponding to (5.2) is given by where again H i (u, v) = u v h i . Notice that this model is indeed manifestly of non-difference form. One function can be absorbed into a normalization (H 1 ) and one can be used in a reparameterization, so we are left with one additional free function. To be more precise, the R-matrix (5.3) is of quasi-difference form, in fact the spectral parameters always
Here φ 1,2 and ψ 1,2 span the two independent su(2) fundamental representations. Explicitly in matrix form, the Hamiltonian density is given by where the h i 's are dependent on the spectral parameter θ.
Similarly, for the R-matrix, we write
Integrable Hamiltonians Upon imposing our integrability constraint we find a total of eight integrable models of non-difference form, in particular six of them have h 3 = 0 and are listed in Table 4.
To our best knowledge all of these models are new. They have some interesting properties. These models at most either exhibit electron pair formation or electron pair splitting, not both. Model 1 can only be made Hermitian if c = 0 and the dependence on the spectral parameter drops out, models 2 to 6 are Hermitian if we impose some conditions on the functions, see Table 5. More generally, we can relate the Hamiltonians of the models 1-6 and their Hermitian conjugate by a unitary transformation if we impose the conditions 3 on Table 6.
Model 5 is a quadruple embedding of model 6-vertex B of 3.1. This can be seen by applying a constant local basis transformation 4 to the Hamiltonian of model 6 V B after the redefinition h 5 → 1 2 h 4 h 5 , and by making the identifications

Model
Reality conditions Table 5: Conditions on models 1-6 to make them Hermitian. We put the superscript (r) to identify the real part of the functions and constants.

Model
Unitarity conditions  Model 7 This model is the most general of non-difference form with h 3 = 0. In fact, we will show that model 8 can be obtained from this one by performing a double limit. In order to solve this model we fixed the normalization of the Hamiltonian such that h 10 = 1 5 . We set h 4 = h 6 = 0 by using the identifications described in 2.2, after which we get the following set of coupled differential equations , (5.14) Summing the first two equations of (5.15) and taking into account the third one, substituting 6 we find that where Ξ i = ξ i , σ = ±1 and c 1 is a constant. To find ξ 1 we made the substitution (5.16) in the differential equation for h 7 and we get For general c 1 the equation (5.18) can be solved by performing the substitutions Ξ 1 (u) → w(u), ξ 1 (u) →ẇ(u) andξ 1 (u) →ẅ(u) and we find that the corresponding differential equation is solved by elliptic functions ξ 1 (u) = i 8 c 2 1 cs(z|m)ds(z|m)ns(z|m), (5.19) where z = i 2 c 1 (u + c 2 ) and m = 8c 3 c 2 1 , c 2,3 are constants. To summarize we then get (5.14) together with 20) The Hamiltonian found does not depend on any free functions, so it should be equivalent to the Hamiltonian of AdS/CFT. In order to prove this, we compared the Hamiltonian that we found with the one derived from requiring centrally extended su(2|2) symmetry [10]. After using an appropriate normalization and shift, the entries of the Hamiltonian of AdS/CFT are 5 We can notice that, if h 10 cannot be normalized to 1, we should impose h 10 = 0 from the beginning and we found a HamiltonianH. This model is equivalent to model 1 under the transformation PH P T . 6 For simplicity we will omit the dependence of ξ 1 , ξ 2 , Ξ 1 and Ξ 2 on the spectral parameter.
where α is a free constant and x + and x − the Zhukovksy variables. x ± can be conveniently parametrized using elliptic functions [40] as 7 we indeed see that the two Hamiltonian densities are the same under The other choice of σ = −1 is not independent, in fact it can be related to the previous one by a twist 8 . So, remarkably, by using our method we are naturally lead to the elliptic parameterization of the AdS/CFT R-matrix.
Model 8 Model 8 can be obtained as a special limit of Model 7. However, since the limit is somewhat singular, let us spell out this case explicitly. If we solve (5.18) for c 1 = 0, we get ξ 1 (u) = c 2 e c 3 u and so It is interesting to notice that the limit c 1 → 0 in the Hamiltonian of model 7 is not well defined because some of the Jacobi functions are divergent in this limit. In order to find the correct results one should take the result for general c 1 and then follow the steps below  (5.26) Comparison between difference and non-difference form models In order to have a complete classification of the models with su(2) ⊕ su(2) symmetry, we compare the models in Table 1 of [20] with the ones in Table 4 using the allowed identifications, i.e. normalization, shift, rescaling and twists. With this, one can see which non-difference form models constructed here reduce to the difference form given in [20]. By doing this comparison we found the correspondence listed in Table 7. We should mention that to verify that model 3 of difference form can be obtained from model 4 of non-difference form one needs to perform a limit because the solution is found for c 2 = 0 (pole for h 5 ) and c 1 = −2.
We furthemore notice that models 2, 4 and 7 of non-difference form generate more than one independent difference form models and that models 1 and 3 do not have a difference form version.
Finally, models 9, 10 and 11 of [20] cannot be obtained from any non-difference form version, so if we want a complete classification of 16 × 16 matrices with su(2) ⊕ su(2) symmetry we should add those three models.
Difference form Non-difference form  In the following c, c 1 , c 2 are constants, F ± = F (u) ± F (v) and similarly for G and H, r i = r i (u, v) and σ = ±1.
Explicitly, the entries of the R-matrices that we got are Model 1  We gave here the full classification of integrable models with su(2) ⊕ su(2) symmetry.
In particular we showed that model 7 corresponds to the AdS 5 × S 5 integrable system derived in [10,41]. This model, as shown in [42] contains the one-dimensional Hubbard model. Our next goal is to address the question of finding new Hubbard-type solutions with a more general form.

Generalised Hubbard model
The Hubbard model is an integrable spin chain with 4 dimensional local Hilbert space identified with two bosons and two fermions. The R-matrix was constructed by Shastry [31] and is notably of non-difference form. It was known early on that the model possesses su(2) ⊕ su(2) symmetry, but this is not enough to uniquely fix the R-matrix. It was later found that the su(2) ⊕ su(2) could be embedded into the centrally extended su(2|2) superalgebra which arises naturally with the worldsheet S-matrix of the AdS 5 ×S 5 integrable system [42] and is equivalent to the Shastry R-matrix. Hence, the Hubbard model is of particular significance in the context of the AdS/CFT correspondence and the corresponding integrable structures.
In this respect we are motivated to look for new integrable models beyond the conventional ones relevant for AdS/CFT integrability potentially associated with new quantum algebras [43,44]. In order to make progress in this direction we implement an ansatz which preserves fermion number but can violate spin conservation as was done in [20] H = H Kin + K Flip + K Pair + V. (5.47) where H Kin is the kinetic term of the free Hubbard model and K Pair and K Flip are further kinetic terms which describe the hopping of a pair of electrons and a term which flips the spins of the electrons on neighbouring sites, respectively, and V is a general potential term, and the precise form of these operators can be found in [20]. The total space of solutions is very large, but we can single out solutions that have the maximal amount of non-zero entries. In particular, demanding that all entries are non-zero, we only find one independent solution. Remarkably it reduces to an R-matrix of difference form. The Hamiltonian is given by

Discussion and conclusions
In this paper we have classified various types of integrable systems and found a plethora of new models generalizing a method based on the boost operator initially put forward in [1]. By starting with a generic Hamiltonian we constrained it to potentially belong to an integrable model by imposing that it commutes with the first higher conserved charge generated by the boost operator. In all cases we showed that this condition is sufficient and we were able to subsequently derive the corresponding R-matrices, guaranteeing integrability. For 4 × 4 models we reviewed the classification of 8-and-lower vertex models originally presented in [21]. We also proved that any Hermitian integrable Hamiltonian can be reconducted (using a local basis transformation) to be 8 V type. Next, we examined 9 × 9 models and completely classified all 15-vertex models satisfying the ice-rule. Finally, for 16 × 16 R-matrices we classified all models with su(2) ⊕ su(2) symmetry. There are various interesting avenues for future research. We discuss some of them here. One natural direction involves applications to holography and the integrable systems which appear in that context. Generalised Shastry-type models provide a base for a search of new types of solutions that are relevant for AdS 4,5 integrable models. In particular it would be interesting to search for new deformations of the AdS 4,5 S-matrix and establish potential contact with q-deformations of the underlying twisted Hopf algebra as for ηdeformed AdS 5 × S 5 [45,46] or for λ-deformed systems which appear to be non-ultralocal [47]. We already discussed in [21] that the AdS 2 S-matrix could be embedded into the 4×4 model 8VB and admitted a one-parameter deformation. Similarly we showed that the S-matrices of AdS 3 governing the scattering of particles with the same chirality could be embedded into both 6VB and 8VB, and both had tunable parameters which correspond to sources of deformations. It would be highly interesting to find a physical interpretation for these parameters and to determine the symmetry algebras of the resulting S-matrices and, in the case of AdS 3 , to check if the deformations of the same chirality S-matrices induce a corresponding deformation of the opposite chirality S-matrices. As well as this, in this paper we demonstrated that there are no integrable deformations of the AdS 5 /CF T 4 Smatrix compatible with su(2) ⊕ su(2) symmetry. However, integrable deformations which do not have this symmetry do exist, for example the q-deformed model [45] and associated quantum algebra [44,48], and it may still be possible to construct integrable deformations which do not have this symmetry by using a perturbative approach. For example, one can start from a given integrable Hamiltonian H 12 . By imposing the commutativity of the first two conserved charges, one will obtain a set of linear equations for the coefficients of H 12 . After factoring out trivial deformations arising from identifications, the resulting set of equations should be quite tractable.
It would be also very interesting to start from Hamiltonians with more vertices and less symmetry. For a Hilbert space of dimension three for example, one could consider 19-vertex Hamiltonians and do a full classification of these models. It is likely that many new non-difference form models could be found. It would be interesting to verify if some of the models we found in the 15-vertex case are actually special cases of more general 19-vertex models.
Furthermore, for Hilbert spaces of dimension 4 we obtained a correspondence between some the non-difference form models found and the difference one of [20]. In particular we found that models 9, 10 and 11 of [20] cannot be obtained from any of the nondifference form models presented in this paper. It would be very interesting to check if these models can be obtained as special limits from other non-difference form models with less symmetry then the su(2) ⊕ su(2) considered here.
Another question that could be addressed is the construction of finite length open spin chains for all the new models constructed in this paper. In order to do that, the first step would be the construction of all possible integrable boundary conditions, meaning all solutions of the Boundary Yang-Baxter equation [49] for each of the R-matrices introduced here.
Remarkably, all of our solutions of the Yang-Baxter equation can be characterized by the integrability condition [Q 2 , Q 3 ] = 0. It is unclear to us why this is the case. Indeed, all the reverse lines in the flowchart, Figure 1, can be shown to hold. The reverse arrows that we exploit here, however, appear to be valid as well and it seems to indicate an equivalence relation. It would be very important to understand and prove these relations.
There are also interesting related mathematical questions to be asked. In the case of difference form models the condition [Q 2 , Q 3 ] = 0 results in a set of cubic polynomial equations for the Hamiltonian entries which seems to be fully equivalent to the Yang-Baxter equation. It would be highly interesting to construct a proof of this claim and in doing so perhaps obtain a closed form expression for the R-matrix in terms of the Hamiltonian entries. In this paper we have relied on a brute force approach to solving the constraint [Q 2 , Q 3 ] = 0 and to a large extent have exhausted the cases where such an approach is applicable.
In order to make more progress it could be important to make use of the extensive toolbox of algebraic geometry. Indeed, [Q 2 , Q 3 ] = 0 describes an algebraic variety in projective space described by a set of coupled, cubic polynomials. For instance, in the 4 × 4 case the integrable models will correspond to algebraic varieties in CP 16 . It would be very interesting to exactly understand what the algebraic varieties are that describe integrable models and how exactly they can be characterized.

(A.3)
Now take (A.2) and multiply from the left with the product of R-matrices . . . R a,k+2 and from the right with R a,k−1 . . . . We then multiply the resulting equation by k and sum over k from −∞ to ∞. The two terms on the right hand side of (A.2) telescopically cancel and we are left with

B Notebook
We have provided a Mathematica notebook with the arxiv submission of this paper. The notebook provides a database of all R-matrices presented in this publication along with those presented in [1,20,21]. A model is defined by three parameters, which we denote in Mathematica as spec, dimHS and model. spec is a list containing the set of spectral parameters of the R-matrix. For a non-difference form R-matrix R(u, v) we have spec = {u,v} and for a differenceform R-matrix R(u) we have spec = {u}.
The parameter dimHS can take the values 2, 3, 4 and specifies the dimension of the local spin chain Hilbert space. If dimHS = n then the corresponding R-matrix is of size n 2 × n 2 .
Finally, model can take the values 0, 1, 2, . . . and specifies, together with dimHS and spec which of the models in the papers [1,20,21] and in this publication we are referring to. The precise map between the value of model and R-matrices is specified in the notebook.

R-matrix
The R-matrix corresponding to a given model as explained above is obtained by running the command where spec={u} for non-difference form and spec={} for difference form. where spec = {u,v,w} for a non-difference form model and spec = {u,v} for difference form. This command evaluates R 12 (u, v)R 13 (u, w)R 23 (v, w) − R 23 (v, w)R 13 (u, w)R 12 (u, v) for non-difference form and R 12 (u−v)R 13 (u)R 23 (v)−R 23 (v)R 13 (u)R 12 (u−v) for difference form. If the Yang-Baxter equation is satisfied the output is {0}.

Yang-Baxter equation
Regularity Regularity is the condition that