Two-parameter families of matrix product operator integrals of motion in Heisenberg spin chains

Tw o-parameter families of MPO in tegrals of motion in Heisen b erg spin c hains Vsev olo d I. Y ashin 1, 2 , ∗ 1 Steklov Mathematic al Institute of R ussian A c ademy of Scienc es, Mosc ow 119991, R ussia 2 Russian Quantum Center, Skolkovo, Mosc ow 143025, R ussia (Dated: F ebruary 24, 2026) Recen tly , F endley et al. (2025) revealed a new wa y to demonstrate the integrabilit y of XYZ Heisen b erg mo del by constructing a one-parameter family of integrals of motion in the matrix pro duct op erator (MPO) form. In this short note, I rep ort on the discov ery of tw o-parameter families of MPOs that commute with with the Heisenberg spin chain Hamiltonian in the XXX, XXZ, and XYZ cases. I describ e a sym b olic algebra approach for finding such in tegrals of motion and speculate about possible applications. I. INTR ODUCTION One-dimensional Heisenberg spin chain models [ 1 ] are very p opular exactly solv able mo dels for studying critical p oin ts and phase transitions of magnetic systems in quan tum many-bo dy ph ysics. There is a long and impactful history of mathematical and ph ysical breakthroughs connected with Heisen b erg mo dels. Early dev elopments include Hans Bethe’s solution to isotropic Heisenberg mo del using his famous ansatz [ 2 ] and its generalization to XXZ anisotropic case [ 3 – 6 ]. Later, Sutherland found the connection of the XYZ Heisen b erg mo del with classical eight-v ertex mo del [ 7 ], and Baxter managed to solv e these mo dels [ 8 – 10 ]. The metho ds developed in these works ev olv ed to become the algebraic Bethe ansatz approach [ 11 , 12 ] and w ere applied to wide range of integrable models [ 13 ]. Bethe ansatz metho ds do not directly lead to concrete description of the local conserv ed quantities of the mo del. The question of explicitly constructing lo cal charges (higher Hamiltonians) w as solved for isotropic XXX Heisenberg spin c hain b y Grabowski and Mathieu [ 14 , 15 ]. They used bo ost operator and com binatorics to construct explicit form of lo cal charges in terms of recursiv e relations. In [ 16 ], Y amada and F uk ai enco ded these conserv ed quantities using matrix product op erator (MPO) [ 17 – 19 ] form. Anisotropic XXZ mo del allo ws for description in terms of T emperley- Lieb algebra, Nienhuis and Huijgen ga ve a closed-form description to XXZ lo cal c harges in terms of T emp erley-Lieb algebra elements [ 20 ]. Also, T emperley-Lieb algebra w as used to give closed-form expression to Flo quet lo cal integrals of motion of the XXZ mo del with anisotropy parameter ∆ = i , whic h is a free-fermionic mo del related to logarithmic conformal field theories [ 21 ]. Nozaw a and F uk ai work ed out the complete description via recursion relations to lo cal c harges of the XYZ mo del [ 22 ]. Recen tly , Paul F endley , Sascha Gehrmann, Eric V ernier and F rank V erstraete found a simplified wa y to sho w in tegrability of the XYZ Heisenberg spin chain [ 23 ]. They prop osed to study in tegrals of motion in MPO form and form ulated a simple sufficient condition on lo cal tensors that makes MPO commute with the Hamiltonian. They gav e an explicit MPO solution dep enden t on one parameter with b ond dimension 4 satisfying these conditions. Expanding the MPO charge as a series ov er this parameter, they obtained even lo cal charges of the XYZ mo del. Their work also con tains results ab out the relation b et ween MPO charges and Bethe ansatz transfer matrices, and generalizations to v arious b oundary conditions and defects. In general, it app ears that the scop e of the metho d should b e broad. Even more recently , F uk ai and Y amada [ 24 ] constructed a one-parameter family of MPO integrals of motion that contains the information about all lo cal charges of the XYZ mo del. This solution help ed them to find a b etter explanation to the combinatorial structure of the lo cal charges. The solution uses dual num ber ε such that ε 2 = 0, which can b e rewritten as the usual MPO with b ond dimension higher than 4. In this work, I exploit the idea of [ 23 ] and some symbolic algebra to find tw o-parameter families of MPO integrals of motions with b ond dimension 4 applicable to XXX, XXZ, and XYZ cases. Expanding this solutions into series, one can obtain lo cal c harges of the spin c hains, and I conjecture the solutions might contain information ab out all in tegrals of motion. The parameters turn out to hav e rather conv enien t geometry , for example the XXX and XXZ solutions dep end on p oin ts of a sphere. Also, unlik e Baxter’s eight-v ertex solutions, the MPO do not require using the theory of sp ecial functions. The work w as done indep endently of [ 24 ]. This pap er is organized as follows. After giving concise preliminaries ab out quantum Heisen b erg spin c hains and their their MPO integrals of motion ( Section I ), I state the main result which consists of writing do wn tw o-parameter families of in tegrals of motion for XXX, XXZ, and XYZ mo dels ( Section II I ). Then I disclose t wo basic ideas that help ed me during the search ( Section IV ), and finally I comment on p ossible future researc h and applications ( Section V ). ∗ yashin.vi@mi-ras.ru 2 I I. PRELIMINARIES Let me start by giving minimal preliminaries necessary for formulating the main result. 1. Hamiltonians Supp ose there are L spin- 1 2 sites, each site lab eled by co ordinate j ∈ { 1 , . . . , L } . I denote Pauli op erators acting on j -th site as { X j , Y j , Z j } so that X 2 j = Y 2 j = Z 2 j = 1 , X j Y j Z j = i for all j (here 1 is an identit y op erator and i is an imaginary unit), and the op erators on different sites commute. I will condsider one-dimensional Heisenberg mo dels with p erio dic b oundary conditions. The Hamiltonians of suc h mo dels in isotropic XXX and anisotropic XXZ and XYZ cases are defined as H X X X = 1 2 L X j =1 [ X j X j +1 + Y j Y j +1 + Z j Z j +1 ] , H X X Z = 1 2 L X j =1 [ X j X j +1 + Y j Y j +1 + ∆ Z j Z j +1 ] , H X Y Z = 1 2 L X j =1 [ J x X j X j +1 + J y Y j Y j +1 + J z Z j Z j +1 ] , (1) where ∆ and J x , J y , J z are coupling constants. F or simplicity , I study translationally inv ariant Heisenberg mo dels with p eriodic b oundary conditions, which means L + 1 ≡ 1. 2. Matrix pr o duct op er ators Giv en a Hamiltonian H , an y operator M commuting with Hamiltonian [ H , M ] = 0 remains constant during the ev olution and is called an inte gr al of motion or a char ge . In this work, I search for integrals of motion given in a matrix pr o duct op er ator (MPO) [ 17 – 19 ] form M = tr[ A 1 A 2 · · · A L ] , (2) where each A j is a 4 × 4-matrix with entries of given b y op erators acting on site j , each entry is defined b y tw o indices from { 0 , x, y , z } : A j =    ( A j ) 0 , 0 ( A j ) 0 ,x ( A j ) 0 ,y ( A j ) 0 ,z ( A j ) x, 0 ( A j ) x,x ( A j ) x,y ( A j ) x,z ( A j ) y , 0 ( A j ) y ,x ( A j ) y ,y ( A j ) y ,z ( A j ) z , 0 ( A j ) z ,x ( A j ) z ,y ( A j ) z ,z    . (3) Ro w and column dimensions of such matrices are called b ond dimensions , in the describ ed case they equal 4. Multi- plication of tw o op erator matrices is defined as the usual matrix m ultiplication and the trace tr is a sum of diagonal en tries. Because I consider translationally in v arian t mo dels, the op erator matrix A j remains constant ov er all sites j . Later in text I might leav e lo wer index and write A for visual clarit y . It is customary and illustrative to depict op erator matrices A as tensor diagr ams [ 25 ] ⟨ n | ( A j ) α,β | m ⟩ = n m α β j , (4) where v ertical lines represen t physic al indic es n, m and horizon tal lines represen t virtual indic es α, β . Connecting lines b et w een diagrams corresponds to contracting the tensor index. The tensor diagram of the matrix product operator 3 M with p eriodic b oundary conditions reads M = · · · 1 2 3 L . (5) 3. Main e quation Heisen b erg model Hamiltonians [ Eq. (1) ] consist of nearest-neighbour interactions and eac h Hamiltonian H is decomp osed as H = P j H j,j +1 , where H j,j +1 is lo cal term acting of tw o sites j, j + 1. Recently , Paul F endley , Sascha Gehrmann, Eric V ernier and F rank V erstraete prop osed to study a simple sufficient condition for MPO in tegrals of motion [ 23 ]: for all sites j , for local Hamiltonian term H j,j +1 and op erator matrices A j there should exist op erator matrices E j called err or terms satisfying i H j,j +1 j j + 1 − i H j,j +1 j j + 1 = E j j + 1 − j E j +1 . (6) By contracting the diagrams with tensors A j and summing ov er j , error terms cancel, giving [ M , H ] = 0. Note that the error terms are not unique: if E j satisfies Eq. (6) , then E j + γ A j also satisfies Eq. (6) for any v alue γ . So, in order to find MPO integrals of motion, one migh t search for solutions of Eq. (6) . I II. MAIN RESUL T Let me write down the main result of the work. I giv e a list of solutions to Eq. (6) for XXX, XXZ, and XYZ mo dels. These solutions dep end on tw o-dimensional space of parameters (excluding scaling). In comparison, the solution prop osed in [ 23 ] is one-parameter and only expresses ev en XYZ lo cal c harges; the one-parameter solution found in [ 24 ] expresses all lo cal charges on XYZ mo del but uses dual num b er ε such that ε 2 = 0. One can chec k that the expressions listed b elo w are indeed solutions to Eq. (6) by direct substitution. I give a commen t on how they were found in Section IV . A. XXX Heisenberg mo del F or isotropic XXX Heisenberg spin chain mo del with Hamiltonian H X X X = 1 2 X j ∈ [ L ] [ X j X j +1 + Y j Y j +1 + Z j Z j +1 ] , (7) there is a parametrized solution to Eq. (6) whic h dep ends on three (p ossibly complex) v ariables w , p, r . The solution reads: A =    w 2 + p 2 + r 2 pw X pw Y pw Z pw X w 2 − r w Z r w Y pw Y r w Z w 2 − r w X pw Z − r w Y r w X w 2    (8) with corresp onding error term E =    r w ( γ + 1)( p 2 + r 2 + w 2 ) − 2 r w pr ( γ + 1) X pr ( γ + 1) Y pr ( γ + 1) Z pr ( γ + 1) X r w ( γ − 1) − ( γ r 2 − p 2 ) Z ( γ r 2 − p 2 ) Y pr ( γ + 1) Y ( γ r 2 − p 2 ) Z r w ( γ − 1) − ( γ r 2 − p 2 ) X pr ( γ + 1) Z − ( γ r 2 − p 2 ) Y ( γ r 2 − p 2 ) X r w ( γ − 1)    , (9) where γ is a parameter describing freedom in choosing error terms. 4 When examining this solution, it might b ecome tempting to introduce spherical co ordinates w = ρ sin θ , p = ρ cos θ cos ϕ, r = ρ cos θ sin ϕ. (10) In terms of these co ordinates, the solution reads A = ρ 2     1 sin θ cos θ cos ϕ X sin θ cos θ cos ϕ Y sin θ cos θ cos ϕ Z sin θ cos θ cos ϕ X sin 2 θ − sin θ cos θ sin ϕ Z sin θ cos θ sin ϕ Y sin θ cos θ cos ϕ Y sin θ cos θ sin ϕ Z sin 2 θ − sin θ cos θ sin ϕ X sin θ cos θ cos ϕ Z − sin θ cos θ sin ϕ Y sin θ cos θ sin ϕ X sin 2 θ     . (11) Here ρ b ecomes a scaling parameter, whic h is insignifican t. So, the solution is dep endent on t wo-parameter family A = A ( θ , ϕ ) where ( θ , ϕ ) enco de p oin ts of a sphere. Multiplying L lo cal tensors A j , one constructs a t wo-parameter family of in tegrals of motion M = M ( θ , ϕ ) as in Eq. (2) . At p oin t θ = 0 and ϕ = 0 the charge is a unit op erator M (0 , 0) = 1. Expanding M near the zero θ = 0 and ϕ = 0, the terms of expansion are lo cal (i.e., sums of close-range interactions) conserved quantities of the mo del: M ( θ, ϕ ) = 1 + θ 2 H + θ 3 ϕH (3) + . . . , (12) where H = H X X X is a Hamiltonian and H (3) is a w eigh t-3 lo cal integral of motion; higher terms in the series are also lo cal charges by construction. I conjecture that M ( θ, ϕ ) is in fact complete, meaning that it contains information ab out any integral of motion of the XXX model, and all lo cal charges of the mo del can be obtained using this expansion. (In comparison, the solution found in [ 23 ] only contains information ab out even charges and corresp onds to setting ϕ = 0.) B. XXZ Heisenberg mo del F or anisotropic XXZ Heisenberg spin chan mo del with Hamiltonian H X X X = 1 2 X j ∈ [ L ] [ X j X j +1 + Y j Y j +1 + ∆ Z j Z j +1 ] , (13) there is parametrized solution to Eq. (6) which dep ends on three (p ossibly complex) v ariables w , p, r : A j =     ∆( w 2 + p 2 + r 2 ) √ ∆ pw X j √ ∆ pw Y j √ ∆ pw π ∆ Z j √ ∆ pw X j w 2 − r w Z j r wπ ∆ Y j √ ∆ pw Y j r w Z j w 2 − r wπ ∆ X j √ ∆ pw π ∆ Z j − r wπ ∆ Y j r wπ ∆ X j w 2 ω ∆     , (14) where I denote π ∆ = p ∆ + η ∆ , η ∆ = ∆ 2 − 1 ∆ w 2 p 2 + r 2 , ω ∆ = r 2 + p 2 ∆ 2 ∆( p 2 + r 2 ) . (15) The corresp onding error term reads E =     r w ( γ + 1)( p 2 + r 2 + w 2 )∆ 2 − 2 rw √ ∆ pr [( γ + 1) + η ∆ ] X √ ∆ pr [( γ + 1) + η ∆ ] Y √ ∆ pr ( γ + 1)∆ π ∆ Z √ ∆ pr [( γ + 1) + η ∆ ] X r w ( γ − 1) − [( γ r 2 − p 2 ) + r 2 η ∆ ] Z ( γ r 2 − p 2 )∆ π ∆ Y √ ∆ pr [( γ + 1) + η ∆ ] Y [( γ r 2 − p 2 ) + r 2 η ∆ ] Z r w ( γ − 1) − ( γ r 2 − p 2 )∆ π ∆ X √ ∆ pr ( γ + 1)∆ π ∆ Z − ( γ r 2 − p 2 )∆ π ∆ Y ( γ r 2 − p 2 )∆ π ∆ X [( γ + 1)∆ ω ∆ − 2] rw     , (16) When ∆ → 1, this solution replicates XXX case Eq. (8) . In tro ducing spherical co ordinates Eq. (10) , it b ecomes A = ρ 2         ∆ √ ∆ sin θ cos θ cos ϕ X √ ∆ sin θ cos θ cos ϕ Y √ ∆ sin θ q ∆ − 1 ∆ sin 2 θ cos ϕ Z √ ∆ sin θ cos θ cos ϕ X sin 2 θ − sin θ cos θ sin ϕ Z sin θ q ∆ − 1 ∆ sin 2 θ sin ϕ Y √ ∆ sin θ cos θ cos ϕ Y sin θ cos θ sin ϕ Z sin 2 θ − sin θ q ∆ − 1 ∆ sin 2 θ sin ϕ X √ ∆ sin θ q ∆ − 1 ∆ sin 2 θ cos ϕ Z − sin θ q ∆ − 1 ∆ sin 2 θ sin ϕ Y sin θ q ∆ − 1 ∆ sin 2 θ sin ϕ X sin 2 θ  1 ∆ sin 2 ϕ + ∆ cos 2 ϕ          . (17) This solution leads to an MPO integral of motion M = M ( θ , ϕ ) which dep ends on points of a sphere ( θ , ϕ ). Lo cal in tegrals of motion can b e found by series expansion near θ = 0 and ϕ = 0. In the next subsection ( Section I I I C ) we will also men tion another solution which lacks the inherent symmetry of the XXZ-mo del. 5 C. XYZ Heisenberg mo del There is a family of integrals of motions for XYZ Heisen b erg mo del with Hamiltonian H X Y Z = 1 2 L X j =1 [ J x X j X j +1 + J y Y j Y j +1 + J z Z j Z j +1 ] , (18) that dep ends on (generally complex) parameters π x , π y , π z and is given by A =         1 q π x Π 2 Π 3 X q π y Π 2 Π 3 Y q π z Π 2 Π 3 Z q π x Π 2 Π 3 X J z π y − J y π z J y π y − J z π z − q π x π y Π 1 Π 3 Z q π z π x Π 1 Π 3 Y q π y Π 2 Π 3 Y q π x π y Π 1 Π 3 Z J x π z − J z π x J z π z − J x π x − q π y π z Π 1 Π 3 X q π z Π 2 Π 3 Z − q π z π x Π 1 Π 3 Y q π y π z Π 1 Π 3 X J y π x − J x π y J x π x − J y π y         (19) where abbreviations Π 1 ,Π 2 ,Π 3 descib e the determinants Π 1 = det   1 1 1 J 2 x J 2 y J 2 z J x π x J y π y J z π z   , Π 2 = det   1 1 1 J 2 x J 2 y J 2 z J x π y π z J y π z π x J z π x π y   , Π 3 = det   1 1 1 J x π x J y π y J z π z J 2 x π 2 x J 2 y π 2 y J 2 z π 2 z   = ( J y π y − J z π z )( J z π z − J x π x )( J x π x − J y π y ) . (20) The corresp onding error terms are E =         γ q Π 1 Π 3 √ π x Π 1 Π 2 γ + J y π y + J z π z Π 3 X p π y Π 1 Π 2 γ + γ + J z π z + J x π x Π 3 Y √ π z Π 1 Π 2 γ + J x π x + J y π y Π 3 Z √ π x Π 1 Π 2 γ + J y π y + J z π z Π 3 X q Π 1 Π 3 γ ( J z π y − J y π z ) − 2( J 2 y − J 2 z ) π y π z J y π y − J z π z − √ π x π y e Π 2 +( γ − J z π z )Π 1 Π 3 Z √ π z π x e Π 2 +( γ − J y π y )Π 1 Π 3 Y p π y Π 1 Π 2 γ + J z π z + J x π x Π 3 Y √ π x π y e Π 2 +( γ − J z π z )Π 1 Π 3 Z q Π 1 Π 3 γ ( J x π z − J z π x ) − 2( J 2 z − J 2 x ) π z π x J z π z − J x π x − √ π y π z e Π 2 +( γ − J x π x )Π 1 Π 3 X √ π z Π 1 Π 2 γ + J x π x + J y π y Π 3 Z − √ π z π x e Π 2 +( γ − J y π y )Π 1 Π 3 Y √ π y π z e Π 2 +( γ − J x π x )Π 1 Π 3 X q Π 1 Π 3 γ ( J y π x − J x π y ) − 2( J 2 x − J 2 y ) π x π y J x π x − J y π y         , (21) where e Π 2 = det   1 1 1 J x π x J y π y J z π z J 3 x π x J 3 y π y J 3 z π z   . (22) This solution defines a set of MPO charges M = M ( π x , π y , π z ) whic h is inv ariant under parameter scaling ( π x , π y , π z ) 7→ ( λπ x , λπ y , λπ z ), it is natural to consider it as dep ending on homogeneous co ordinates [ π x : π y : π z ] describing p oin ts of (complex) pro jectiv e plane. 1. XXZ limit In the limit of XXZ mo del J x = 1, J y = 1, J z = ∆ this solution b ecomes A =         1 p π x Π X q π y Π Y q 1 Π Z p π x Π X ∆ π y − 1 ∆ − π y − q − π x π y Π Z p − π x Π Y q π y Π Y q − π x π y Π Z ∆ π x − 1 ∆ − π x − q − π y Π X q 1 Π Z − p − π x Π Y q − π y Π X 1         (23) where I set π z = 1 and denote Π = (∆ − π x )(∆ − π y ) ∆ 2 − 1 . (24) 6 This solution leads to a family of MPO charges M = M ( π x , π y ) where it is natural to treat parameters ( π x , π y ) as deco ding p oints of a (p ossibly complex) plane. The drawbac k of this XXZ solution is that the tensor A is not symmetric under rotations ov er Z , as we will discuss in Section IV A . In the XXX case ∆ → 1, the op erator matrix A tends to iden tity matrix, thus becomes trivial; in the XX case ∆ = 0 it gives a non-trivial solution. 2. Two one-p ar ameter solutions Unfortunately , I was not able find a “sphere-parametered” solution that w ould generalize Eq. (14) and I leav e its existence (or non-existence!) as a conjecture. As a partial result asso ciated with this hypothetic solution, let me note that there are tw o one-parameter solutions. The first is the solution of [ 23 ], which can b e rewritten for a parameter v as A =     v p J x v − J y J z X p J y v − J z J x Y p J z v − J x J y Z p J x v − J y J z X J x 0 0 p J y v − J z J x Y 0 J y 0 p J z v − J x J y Z 0 0 J z     (25) with the error term E =     0 0 0 0 0 0 p J x v − J y J z p J y v − J z J x Z − p J z v − J x J y p J x v − J y J z Y 0 − p J x v − J y J z p J y v − J z J x Z 0 p J y v − J z J x p J z v − J x J y X 0 p J z v − J x J y p J x v − J y J z Y − p J y v − J z J x p J z v − J x J y X 0     + γ A . (26) The other one-parameter solution (dep enden t on parameter ξ ) is A =          0 0 0 0 0 1 J x − √ ξ − J 2 z J z √ J x J y Z √ ξ − J 2 y J y √ J z J x Y 0 √ ξ − J 2 z J z √ J x J y Z 1 J y − √ ξ − J 2 x J x √ J y J z X 0 − √ ξ − J 2 y J y √ J z J x Y √ ξ − J 2 x J x √ J y J z X 1 J z          , (27) with the error term E =          0 0 0 0 0 0 − √ ξ − J 2 x √ ξ − J 2 y J z √ J x J y Z √ ξ − J 2 z √ ξ − J 2 x J y √ J z J x Y 0 √ ξ − J 2 x √ ξ − J 2 y J z √ J x J y Z 0 − √ ξ − J 2 y √ ξ − J 2 z J x √ J y J z X 0 − √ ξ − J 2 z √ ξ − J 2 x J y √ J z J x Y √ ξ − J 2 y √ ξ − J 2 z J x √ J y J z X 0          + γ A . (28) If the “sphere-parametered” solution exists, I think it should generalize these tw o one-parameter solutions. IV. METHODS Let me explain the metho ds that were used during the search. F rist, I comment on symmetry considerations that help to choose a form of A j , then I describ e ho w to reduce the problem to a set of algebraic equations to solve. A. Considerations regarding symmetry The idea of F endley , Gehrmann, V ernier and V erstraete to searc h MPO in tegrals of motion and to solve Eq. (6) should w ork general one-dimensional quantum chain mo dels with v arious b oundary conditions and defects. So, it is 7 in teresting to study some refinements of the metho ds when studying concrete mo dels. Here I discuss ho w one could come up with to the MPO form that I consider (although to b e honest, one is more likely to adopt from [ 23 , 24 ]). Heisen b erg mo dels include v arious symmetries, it is natural to imp ose those symmetries on the integrals of motion. First of all, there is a translational inv ariance that suggests the tensor A j to b e indep endent of the site j . T ranslationally inv ariant op erator matrices A has a gauge fr e e dom preserving MPO charge M : A 7→ S A S − 1 , 7→ S S − 1 , (29) where S is some complex matrix. Note that the symmetry can partly fix the gauge. The XYZ Heisenberg mo del has Z 2 × Z 2 symmetry given b y rotations Z ⊗ L , Z ⊗ L , Z ⊗ L acting on all spins. F or MPO M to b e in v arian t under s uc h symmetries, we impose on A a condition that physical symmetries U translate to gauge symmetries S : U A α,β U † = X µ,ν S α,µ A µ,ν ( S − 1 ) ν,β , U U † = S S − 1 . (30) That means, the group of symmetries acts as a representation on the space of operator matrices A . Cho osing such represen tation also leads to the c hoice of appropriate b ond dimension for A . In the case w e consider, the b ond dimension is 4 and the Z 2 × Z 2 -represen tation is given by O 7→ X O X acts as A 7→    1 1 − 1 − 1    A    1 1 − 1 − 1    , O 7→ Y O Y acts as A 7→    1 − 1 1 − 1    A    1 − 1 1 − 1    , O 7→ Z O Z acts as A 7→    1 − 1 − 1 1    A    1 − 1 − 1 1    . (31) An op erator matrix A resp ects this Z 2 × Z 2 -symmetry if and only if it has form A =    v q x X q y Y q z Z p x X w x − r x,y Z r x,z Y p y Y r y ,x Z w y − r y ,z X p z Z − r z ,x Y r z ,y X w z    (32) where all the app eared eleme n ts v , p • , q • , r • , w • (together with low er indices) are complex v ariables. No w, the gauge symmetry S comm utes with Z 2 × Z 2 -symmetry if and only if it is a diagonal matrix S = diag( s 0 , s x , s y , s z ). The gauge is finally fixed by ensuring that q x = p x , q y = p y , q z = p z . I also require (p erhaps not entirely justifiably) that the elements r • are anti-symmetric, so that the op erator matrix A has form A =    v p x X p y Y p z Z p x X w x − r z Z r y Y p y Y r z Z w y − r x X p z Z − r y Y r x X w z    (33) and is parametrized by 10 v ariables { v , w x , w y , w z , p x , p y , p z , r x , r y , r z } . I also enco de the error term E in the same form: E j =    e v e p x X j e p y Y j e p z Z j e p x X j e w x − e r z Z j e r y Y j e p y Y j e r z Z j e w y − e r x X j e p z Z j − e r y Y j e r x X j e w z    (34) 8 and parametrize it with 10 v ariables { e v , e w x , e w y , e w z , e p x , e p y , e p z , e r x , e r y , e r z } . The U(1) ⋉ Z 2 -symmetry of XXZ mo del means that the tensor A has form A =    v p x X p x Y p z Z p x X w x − r z Z r x Y p x Y r z Z w x − r x X p z Z − r x Y r x X w z    (35) and dep ends on 7 parameters { v , w x , w z , p x , p z , r x , r z } . Note that the solution Eq. (23) do es not follow this symmetry . The SU(2)-symmetry of XXX mo del means that the tensor A is written as A =    v p X p Y p Z p X w − r Z r Y p Y r Z w − r X p Z − r Y r X w    (36) and dep ends on 4 parameters { v , w , p, r } . B. System of equations Ha ving parametrized matrices A and E [see Eqs. (33) and (34) ], one can substitute A and E to main equation Eq. (6) , obtaining a systems of algebraic equations: J z p z r y + J y p y r z = v ˜ p x − ˜ v p x , J y p y p z − J z r y r z = r x ˜ w y − w y ˜ r x , J y p z r y + J z p y r z = w x ˜ p x − p x ˜ w x , J z p y p z − J y r y r z = r x ˜ w z − w z ˜ r x , p x ( v J z − J y w x ) = r y ˜ p z − p z ˜ r y , r x ( J z w y − J y w z ) = r y ˜ r z − r z ˜ r y , p x ( v J y − J z w x ) = r z ˜ p y − p y ˜ r z , r x ( J y w y − J z w z ) = p y ˜ p z − p z ˜ p y , and cyclic p erm utations of x, y , z . (37) This is a set of 24 homogeneous quadratic equations in 20 v ariables. As shown in examples from Section I II , the solution set contains 4-dimensional comp onen ts (tw o meaningful parameters, one scaling parameter and one error term parameter γ ). The v ariables of E fit linearly in the system, and they can b e eliminated for example using substitutions e w x → e p x p x w x − J y r y p z + J z r z p y p x , e r x → e p z p z r x − J x p y w y − v J z p y p z , e p x → e v v p x + J y r y p z + J z r z p y v , and cyclic p erm utations of x, y , z . (38) In principle, all solutions to system Eq. (37) can b e fully characterized using symbolic algebra and Gr¨ obner basis metho ds [ 26 ]. In practice, for XXX and XXZ cases the Gr¨ obner bases are quite easy to obtain, but the XYZ case seems to require to o muc h time and memory resources to b e computed on a laptop. Also, approximate n umerical solutions can be found using dev elop ed methods of semidefinite programming and rank minimization [ 27 ]. I found solutions describ ed in Section I I I with the help of computer algebra systems by sequen tially choosing conv enien t c hanges of v ariables. V. OUTLOOK As announced in the abstract, the scop e of this work is limited to rep orting on the found solutions. Let me discuss the questions left for future research and sp eculate on p ossible applications. First of all, the listed solutions should b e studied more deeply . It is in teresting to chec k if solutions Eqs. (8) and (14) completely describ e all lo cal integrals of motion for XXX and XXZ mo dels. That study migh t as w ell give insights on combinat orial identities b etw een lo cal charges [ 14 , 16 , 22 , 24 ]. 9 Supp ose that the completeness of the family of the integrals of motion is pro ven. This result can help dev eloping new metho ds for studying integrabilit y . Given a quan tum state ρ , its tra jectory in time would b e completely characterized b y the set of av erage v alues ⟨M ( θ , ρ ) ⟩ ρ . It would then b e p ossible to study the dynamical prop erties of the mo del it terms of function ( θ , ϕ ) 7→ ⟨M ( θ , ρ ) ⟩ ρ defined on a sphere. This function can b e easily computed for pro duct states or for matrix pro duct states, sometimes even in the thermo dynamic limit. Other interesting direction would b e to use the approach of generalized Gibbs ensembles [ 28 – 31 ] in application to in tegrable mo dels. In this approach, the known integrals of motion of the system are included into the definition of thermo dynamic ensemble, allowing for describing equilibrium states using generalized temp eratures. Pro vided that w e ha ve an analytical characterization of all c harges in an integrable mo del, p ossibly one can describ e an y prop erties of the system using generalized temp eratures. It w ould b e in teresting to chec k if there are other interesting MPO integrals of motions to Heisenberg mo dels. In the XXX case, it is p ossible to list all b ond dimension 4 solutions and find that there are no other in teresting ones except Eq. (8) . As mentioned in Section I I I C , I was not able to find a “sphere-parametered” XYZ solution with b ond dimension 4 that would generalize XXX and XXZ solutions, and I am not sure if such solution exists. The authors of [ 24 ] managed to find a one-parameter MPO solution to XYZ Heisenberg mo del using dual num bers ε suc h that ε 2 = 0. Despite the great b eaut y and strength of their deriv ations, the dual n um b ers are not to o easy to interpret physically . Using the encoding ε as a 2 × 2 m atrix, their solution becomes a b ond dimension 8 matrix, p ossibly not resp ecting the Z 2 × Z 2 symmetry – so, searching b ond dimension 4 solutions might hav e b enefits. Certainly , the metho ds describ ed in this work can b e applied to more general mo dels of man y-b o dy physics: most imp ortan tly , to describ e the integrals of motion of the SO(4)-inv ariant one-dimensional Hubbard mo del [ 32 , 33 ]. The authors of [ 24 ] announced that they are working in this direction. The metho ds of studying MPO integrals of motion introduced in [ 23 ] app ear to b e quite general. They seem to b e applicable at least to general one-dimensional spin chains with nearest-neigh b our in teractions, with arbitrary lo cal dimensions, b oundary conditions, defects. It would b e nice to find their application for more “disordered” systems than translationally-in v ariant integrable spin chains. As a natural place of application, one might try to searc h for Flo quet charges in in tegrable brick-w all Flo quet proto cols [ 21 , 34 – 37 ]. Moreo ver, the Eq. (6) can b e generalized to study not only the integrals of motion, but also the exact dynamics of observ ables preserving MPO form during the evolution. Supp ose that we are given an MPO M defined by a sequence of op erator matrices A j o ver sites j , and supp ose that there exist op erator matrices L j and R j +1 suc h that i H j,j +1 j j + 1 − i H j,j +1 j j + 1 = L j j + 1 − j R j +1 . (39) In this case, Heisenberg equation guarantees that the MPO form of M = M ( t ) is preserved during the evolution and d dt A j = L j − R j , at least for small times t . It should b e natural to trial this approach on XX-mo del [ 38 ]. Preliminary in vestigation indicates that such MPO-preserving observ ables to exist, they should break the global symmetry of the in tegrable mo del. In Section IV B , I discussed that the solutions to Eq. (6) reduce to a set of quadratic algebraic equations, and suc h equations could b e solved numerically using semidefinite optimization. I assume that the dev elopmen t of suc h metho ds might b e helpful e.g. for numerically searching approximate integrals of motion. Finally , it is p ossible to search for applications to quantum computing, esp ecially for studying structural prop erties of bric k-w all circuits and v ariational ansatzes [ 39 , 40 ]. Note that the w ell-studied class of dual unitary circuits includes lo cal XXZ interactions [ 41 , 42 ], making solution Eq. (14) somewhat more relev an t. A CKNOWLEDGEMENTS I am grateful to Oleg V. Lyc hko vskiy for his interest in the topic and for useful discussions, as well as to Denis V. Kurlov for introducing me to sym b olic algebra methods in many-bo dy physics. The w ork w as supp orted by the F oundation for the Adv ancemen t of Theoretical Physics and Mathematics “Basis” (pro ject no. 25-1-5-45-1). [1] W erner Heisenberg. “Zur theorie des ferromagnetismus”. Zeitschrift f ¨ ur Physik 49 , 619–636 (1928). 10 [2] Hans Bethe. “On the theory of metals. 1. eigen v alues and eigenfunctions for the linear atomic c hain”. Zeitschrift f¨ ur Ph ysik 71 , 205–226 (1931). [3] R. Orbach. “Linear antiferromagnetic chain with anisotropic coupling”. Ph ys. Rev. 112 , 309–316 (1958). [4] C. N. Y ang and C. P . Y ang. “One-dimensional c hain of anisotropic spin-spin interactions. I. pro of of Bethe’s hypothesis for ground state in a finite system”. Phys. Rev. 150 , 321–327 (1966). [5] C. N. Y ang and C. P . Y ang. “One-dimensional chain of anisotropic spin-spin interactions. II. prop erties of the ground-state energy p er lattice site for an infinite system”. Phys. Rev. 150 , 327–339 (1966). [6] C. N. Y ang and C. P . Y ang. “One-dimensional chain of anisotropic spin-spin in teractions. I I I. applications”. Ph ys. Rev. 151 , 258–264 (1966). [7] Bill Sutherland. “Two-dimensional hydrogen b onded crystals without the ice rule”. Journal of Mathematical Physics 11 , 3183–3186 (1970). [8] Ro dney Baxter. “Eight-v ertex mo del in lattice statistics and one-dimensional anisotropic Heisenberg c hain. I. some fun- damen tal eigenv ectors”. Annals of Ph ysics 76 , 1–24 (1973). [9] Ro dney Baxter. “Eight-v ertex mo del in lattice statistics and one-dimensional anisotropic Heisenberg chain. I I. equiv alence to a generalized ice-type lattice model”. Annals of Physics 76 , 25–47 (1973). [10] Ro dney Baxter. “Eigh t-v ertex model in lattice statistics and one-dimensional anisotropic Heisen berg c hain. I I I. eigen v ectors of the transfer matrix and hamiltonian”. Annals of Physics 76 , 48–71 (1973). [11] L. A. T akhtadzh yan and L. D. F adeev. “The quan tum metho d of the inv erse problem and the Heisenberg xyz mo del”. Russian Math. Surveys 34 , 11–68 (1979). [12] N. A. Slavno v. “Algebraic Bethe ansatz” (2019). . [13] R.J. Baxter. “Exactly solved mo dels in statistical mechanics”. Dov er b ooks on physics. Do ver Publications. (2007). url: https://store.doverpublications.com/products/9780486462714 . [14] Marek P . Grab o wski and Pierre Mathieu. “Quantum integrals of motion for the Heisenberg spin chain”. Mo dern Physics Letters A 09 , 2197–2206 (1994). [15] M.P . Grab owski and P . Mathieu. “Structure of the conserv ation la ws in quantum integrable spin chains with short range in teractions”. Annals of Physics 243 , 299–371 (1995). [16] Kyoic hi Y amada and Kouhei F uk ai. “Matrix pro duct op erator representations for the lo cal conserv ed quantities of the Heisen b erg chain”. SciP ost Physics Core 6 (2023). [17] F. V erstraete, J. J. Garc ´ ıa-Rip oll, and J. I. Cirac. “Matrix pro duct density op erators: Simulation of finite-temp erature and dissipative systems”. Physical Review Letters 93 (2004). [18] Michael Zwolak and Guifr ´ e Vidal. “Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent sup eroperator renormalization algorithm”. Physical Review Letters 93 (2004). [19] B Pirvu, V Murg, J I Cirac, and F V erstraete. “Matrix pro duct op erator representations”. New Journal of Physics 12 , 025012 (2010). [20] Bernard Nienhuis and Onno E Huijgen. “The lo cal conserved quan tities of the closed XXZ chain”. Journal of Physics A: Mathematical and Theoretical 54 , 304001 (2021). [21] Vsevolod I. Y ashin, Denis V. Kurlov, Aleksey K. F edorov, and Vladimir Gritsev. “In tegrable Flo quet systems related to logarithmic conformal field theory”. SciPost Physics 14 (2023). [22] Y uji Nozaw a and Kouhei F uk ai. “Explicit construction of local conserved quantities in the XYZ spin-1 / 2 chain”. Ph ysical Review Letters 125 (2020). [23] Paul F endley , Sascha Gehrmann, Eric V ernier, and F rank V erstraete. “XYZ integrabilit y the easy wa y” (2025). arXiv:2511.04674 . [24] Kohei F uk ai and Kyoic hi Y amada. “Matrix pro duct op erator representations for the lo cal conserved quantities of the spin-1 / 2 XYZ chain” (2026). . [25] Bram V ancraeynest-De Cuip er, W eronik a Wiesiolek, and F rank V erstraete. “Les Houches lecture notes on tensor net- w orks” (2026). . [26] David A. Cox, John Little, and O’Shea Donal. “Ideals, v arieties, and algorithms: An introduction to computational algebraic geometry and commutativ e algebra”. Undergraduate T exts in Mathematics . Springer Cham. (2025). 5 edition. [27] Alex Lemon, Anthon y Man-Cho So, and Yin yu Y e. “Low-rank semidefinite programming: Theory and applications”. F ound. T rends Optim. 2 , 1–156 (2016). [28] Lev Vidmar and Marcos Rigol. “Generalized Gibbs ensem ble in integrable lattice mo dels”. Journal of Statistical Mec hanics: Theory and Exp erimen t 2016 , 064007 (2016). [29] Bal´ azs Pozsga y . “The generalized Gibbs ensemble for Heisenberg spin chains”. Journal of Statistical Mechanics: Theory and Exp erimen t 2013 , P07003 (2013). [30] Marius Costeniuc, Richard S. Ellis, Hugo T ouc hette, and Bruce T urkington. “The generalized canonical ensemble and its univ ersal equiv alence with the micro canonical ensemble”. Journal of Statistical Physics 119 , 1283–1329 (2005). [31] M. Costeniuc, R. S. Ellis, H. T ouc hette, and B. T urkington. “Generalized canonical ensembles and ensem ble equiv alence”. Ph ysical Review E 73 (2006). [32] F abian HL Essler, Holger F rahm, F rank G¨ ohmann, Andreas Kl ¨ umper, and Vladimir E Korepin. “The one-dimensional Hubbard mo del”. Cambridge Univ ersit y Press . (2005). [33] Kohei F uk ai. “All lo cal conserved quantities of the one-dimensional Hubbard mo del”. Physical Review Letters 131 (2023). [34] Vladimir Gritsev and Anatoli Polk ovnik o v. “Integrable Flo quet dynamics”. SciPost Ph ysics 2 (2017). [35] Matthieu V anicat, Lenart Zadnik, and T oma ˇ z Prosen. “Integrable Trotterization: Lo cal conserv ation laws and boundary driving” (2017). . 11 [36] Y uan Miao, Vladimir Gritsev, and Denis V. Kurlov. “The Flo quet baxterisation”. SciPost Physics 16 (2024). [37] Pietro Ric helli, Kareljan Schoutens, and Alb erto Zorzato. “Bric k wall quantum circuits with global fermionic symmetry”. SciP ost Physics 17 (2024). [38] Alexander T eretenk o v and Oleg Lychk ovskiy . “Exact dynamics of quantum dissipative XX mo dels: Wannier-Stark lo cal- ization in the fragmented operator space”. Physical Review B 109 (2024). [39] Matthew P .A. Fisher, V edik a Khemani, Adam Nah um, and Sagar Vijay . “Random quan tum circuits”. Annual Review of Condensed Matter Physics 14 , 335–379 (2023). [40] Romain V asseur. “Les Houches lectures on random quantum circuits and monitored quantum dynamics” (2026). arXiv:2602.17258 . [41] Bruno Bertini, P a vel Kos, and T omaˇ z Prosen. “Exact correlation functions for dual-unitary lattice mo dels in 1 + 1 dimensions”. Physical Review Letters 123 (2019). [42] Lorenzo Piroli, Bruno Bertini, J. Ignacio Cirac, and T omaˇ z Prosen. “Exact dynamics in dual-unitary quantum circuits”. Ph ysical Review B 101 (2020).