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 $({\mathbb 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^2 + y^2=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_1,y_1)$ and $M(x_2,y_2)$ 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.