On integrable matrix product operators with bond dimension D = 4

We construct and study a two-parameter family of matrix product operators of bond dimension D = 4. The operators M(x, y) act on ({\Bbb{C}}_2)^{\otimes N} , i.e. the space of states of a spin-1/2 chain of length N. For the particular values of the parameters: x = 1/3 and y=1/\sqrt{3} , the operator turns out to be proportional to the square root of the reduced density matrix of the valence-bond-solid state on a hexagonal ladder. We show that M(x, y) has several interesting properties when (x, y) lies on the unit circle centered at the origin: x² + y² = 1. In this case, we find that M(x, y) commutes with the Hamiltonian and all the conserved charges of the isotropic spin-1/2 Heisenberg chain. Moreover, M(x₁, y₁) and M(x₂, y₂) are mutually commuting if x^2_i + y^2_i=1 for both i = 1 and 2. These remarkable properties of M(x, y) are proved as a consequence of the Yang-Baxter equation.