Matrix product operator representations for the local conserved quantities of the Heisenberg chain

We present the explicit expressions for the matrix product operator (MPO) representation for the local conserved quantities of the Heisenberg chain. The bond dimension of the MPO grows linearly with the locality of the charges. The MPO has more simple form than the local charges themselves, and their Catalan tree patterns naturally emerge from the matrix products. The MPO representation of local conserved quantities is generalized to the integrable $\mathrm{SU}(N)$ invariant spin chain.


Introduction
Quantum integrable models are special many-body systems that allow for exact solutions [1,2].The Bethe Ansatz, with its origins tracing back to Hans Bethe's seminal work on the exact solution of the spin-1/2 Heisenberg spin chain [3], enables the exact calculation of energy spectrum and physical observables.Through various generalizations, the Bethe Ansatz has become the most renowned method for solving integrable models.
The defining feature of quantum integrable systems is the presence of an extensive number of local conserved quantities, denoted as {Q k } k=2,3,4,... .They are local in the sense that they are a linear combination of operators that act on a finite range of local sites.In our notation here, Q k is the site translation sum of the local operator that acts on the adjacent k sites.The existence of these local conserved quantities has been well established by the quantum inverse scattering method [4]: the local conserved quantities can be derived from the expansion of the transfer matrix T (λ) with respect to the spectral parameter λ, given by log T (λ) ∼ k≥2 λ k−1 Q k , where Q 2 is usually the Hamiltonian itself.The commutativity of the transfer matrix, [T (λ), T (µ)] = 0, ensures the mutual commutativity of the local conserved quantities, [Q k , Q l ] = 0. Another way to obtain Q k is the usage of the Boost operator, denoted by B, if it exists.Q k can be calculated recursively by Although the formal methodology for generating local conserved quantities Q k , through the expansion of the transfer matrix and using the Boost operator, has been known, determining their general expressions in practice remains a formidable challenge.This difficulty stems not only from the excessive computational expense associated with higher-order charges but also from finding a general pattern within the vast data sets generated by these calculations.The general forms of local conserved quantities were firstly found for the spin-1/2 Heisenberg chain(XXX chain) independently by [5] and [6] and subsequently found for its SU(N ) generalization [7].For these models, the structure of the local charges has the Catalan tree pattern [8]: Q k is constructed from the linear combination of the polynomial of spin operators with the coefficient of generalized Catalan number.More recently, the general forms of the local charges were found in the spin-1/2 XYZ chain [9], the Temperly-Lieb models [10], which include the spin-1/2 XXZ chain, and the one-dimensional Hubbard model [11].Nonetheless, their expressions are still slightly complicated, even for these models where the general forms of the local charges are now known.A universal, more simple description of the general form of local conserved quantities is highly desirable.
Over the past two decades, there has been a growing trend to introduce tensor network techniques [12][13][14] to reformulate quantum integrability.The matrix product state(MPS) representation of the Bethe eigenstate of quantum integrable systems has been studied [15][16][17][18][19][20][21][22].The symmetry of integrable systems was also investigated in terms of matrix product operator(MPO).In [23], the MPO commuting with the Hamiltonian of the Heisenberg chain was constructed, which was proved to be the product of two transfer matrices with different spectral parameters.More recently, the non-commutative symmetry of models with fragmented Hilbert space was obtained in the form of MPO, even for non-integrable cases [24].The hidden symmetry of an integrable Lindblad system, which leads to multiple non-equilibrium steady states, was exactly given by MPO [25].However, despite the fact that the transfer matrix, which is the source of Q k , is defined by the MPO constructed from the Lax operator, there has yet to be an exploration of the MPO representation for local conserved quantities themselves.
For example, the transfer matrix of the SU(N) invariant chain is the MPO with bond dimension N .It is worth noting that the fact that the transfer matrix is expressed as an MPO had been known even before the term "MPO" was coined.
In this work, we present the MPO representation of the local conserved quantities for the spin-1/2 Heisenberg spin chains.The local conserved quantity Q k can be represented by the MPO whose bond dimension is 3k−1.We found the Catalan tree pattern of the local conserved quantities [8] naturally emerges from the product of the MPO by using the identity of the generalized Catalan number.Unlike the local conserved quantities themselves, their MPO representation is only involved with the usual Catalan number, implying that the complexity of these expressions is folded within the product of the MPO.This MPO representation of the local conserved quantities can be immediately generalized for the integrable SU(N ) invariant spin chains of fundamental representation.To our knowledge, this is the first study investigating the MPO representation for local conserved quantities of quantum integrable systems.
This paper consists of the following Sections: in Section 2, we review the local charges for the Heisenberg chain and the SU(N ) integrable spin chain.In Section 3, we present the main result of the MPO representation for the local charges in the Heisenberg chain.In section 4, we give the demonstration that the MPO introduced in section 3 actually reproduces the local charges of [5][6][7].Section 5 contains the summary of our work and future outlooks.

Local conserved quantities of the Heisenberg chain
In this section, we review the result of the local conserved quantities of the Heisenberg chain [5,6] and its SU(N ) generalization [7].

spin-1/2 Heisenberg chain
The Hamiltonian of the spin-1/2 Heisenberg chain is given by: where σ i = (X i , Y i , Z i ) stands for the vector of the usual Pauli matrices acting non-trivially on i-th site, and L is the system size.We assume the periodic boundary condition: σ i+L = σ i .The Hamiltonian (1) is integrable [3] and has an extensive number of local conserved quantities {Q k } k=2,3,4,... .The key sign of its integrability is the mutual commutativity To represent the expression of Q k , we introduce some notations.A sequence of n sites C = {i 1 , i 2 , . . ., i n } with i 1 < i 2 < . . .< i n , will be called a cluster of order n.A cluster C can be further classified by hole, defined by i n − i 1 + 1 − n, which is the number of the sites between i 1 and i n that are not included in C. For example, the cluster C = {1, 3, 4, 7} is the cluster of order 4, and whose hole is 3.
For a cluster C = {i 1 , i 2 , . . ., i n } of order n, we define the nested products of Pauli matrices f n (C) by: Then we define components of the local conserved quantities: where C (n,m) denotes the set of clusters of order n with m holes, satisfying 1 ≤ i 1 ≤ L.
The general expression of Q k is given by [5,6]: where is the generalized Catalan number 1 .

SU(N ) generalization
The structure of the local conserved quantities of the isotropic SU(N ) version of the spin-1/2 Heisenberg chain in the fundamental representation is the same as that of the spin-1/2 Heisenberg chain [7].The Hamiltonian of SU(N ) invariant chain is given by [26]: where t a i , a = 1, . . ., N 2 − 1 are the su(N ) generators in the fundamental representation.For the N = 2 case, (5) reduces to the Hamiltonian of the spin-1/2 Heisenberg chain.We choose the normalization of the generators so that t a 's are the su(N ) Gell-Mann matrix, satisfying the following algebra: where f a bc is the structure constant of su(N ), and d a bc is a completely symmetric tensor, which is non-trivial for N > 2.
For the SU(N ) case, f n (C) for a cluster C = {i 1 , i 2 , . . ., i n } is defined by: where t i = t 1 i , . . ., t N 2 −1 i stands for the vector of the su(N ) Gell-Mann matrices acting nontrivially on i-th site.The outer product of the vector A, B with N 2 − 1 elements is defined by The local conserved quantities of the SU(N ) invariant chain are expressed in the same form as (4), with f n (C) defined in (8).

Doubling-product representation for SU(2) case
For the SU(2) case, the local conserved quantities can be represented using the doublingproduct notation, which was initially introduced in the proof of the non-integrability of the spin-1/2 XYZ chain in a magnetic field [27], and subsequently employed to construct the local conserved quantities of the spin-1/2 XYZ chain without a magnetic field [9].
A doubling-product is a notation for the product of the Pauli matrices, defined by: 1 The coefficients of local charges of Eq. (4.2) in [7] are slightly complicated; thus, we used more simple notation here.The relation between our C k,n and α k,l of Eq. (4.2) in [7] is C k,n = α k+1,k+1−n .This can be easily seen by using Eq.(4.5) in [7]: Note that the definition of the symbol C n,m employed in [7] is distinct from the generalized Catalan number in this work.
where A i ∈ {X , Y, Z} and A l A l+1 is the product of A l and A l+1 .(•) i denotes the operator acting on the i-th site.The hole is equal to the number of j that satisfies A j = A j+1 .We define the support of an operator as the range of sites on which it acts.The support of the doublingproduct of ( 9) is n + 1.
With the doubling-product notation, components of local conserved quantities for the spin-1/2 Heisenberg chain (3) can be rewritten by: A, (10) where S l,m is the set of all doubling-products with (support, hole) = (l, m).The local conserved quantities constructed with (10) differ from those constructed with (3) by the factor i n .

MPO representation for the local conserved quantities
In this section, we present the matrix product operator (MPO) representation for the local conserved quantities of the spin-1/2 Heisenberg chain and its SU(N ) generalizations.Given an operator as a sum of finite-range interactions, it is possible to construct the MPO representation for such a local operator using entirely upper (or lower) triangular matrices [28].We show the local conserved quantities can be represented by an open boundary MPO constructed with the upper triangular matrix Γ i k .

MPO for the Hamiltonian
MPO component for the Hamiltonian of the spin-1/2 Heisenberg chain (1) has been known as the matrix with bond dimension 5, the element of which are local operators acting on the physical Hilbert space [29]: where ) is treated as a row vector and σ ⊤ i is its transpose and I is the identity operator.O denotes the zero-operator, and O m,n is an m × n matrix whose entries are all O.
The Hamiltonian is reproduced by the MPO constructed from Γ i 2 : where The bulk term can be written simply with the boundary vectors 〈L| = (1, 0, 0, 0, 0) and |R〉 = (0, 0, 0, 0, 1) ⊤ : The Hamiltonian under the periodic boundary condition ( 1) is reproduced by the sum of the bulk term and boundary term: In the following, we generalize this result to the higher-order local conserved quantities Q k .

Building blocks of MPO
We introduce a 3 × 3 matrix M i , components of MPO which are local operators acting on the physical Hilbert space.M i serves as the building block of the MPO for the local conserved quantities of the spin-1/2 Heisenberg chain.M i is defined by: where the off-diagonal elements are the Pauli matrices acting on the i-th site.Note that M i has the form of the so(3) generator assigned with the Pauli matrices.The nested product of ( 2) is represented with these building blocks by: where in the second line, ) is treated as a row vector, and σ ⊤ i is its transpose.For the more general SU(N ) invariant chain, the building block for the MPO becomes the N 2 − 1 by N 2 − 1 matrix defined by: where the indices run a, c = 1, 2, . . ., N 2 −1.We note that M The nested product (8) for the SU(N ) case is represented by: where ) is treated as a row vector, and t ⊤ i is its transpose.We introduce another building block for the MPO, the 3 by 3 diagonal matrix e for the SU(2) case and the N 2 − 1 by N 2 − 1 diagonal matrix e (N ) for the general SU(N ) cases, the diagonal elements of which are the identity operators: where the non-diagonal elements are all O.We note that e (2) = e.We give the expressions for the local charges in the L = 4 case using the building blocks in appendix D.

Explicit expressions of MPO
Next, we construct the MPO for the local conserved quantities from the building block defined above.While we focus on the spin-1/2 Heisenberg chain (SU(2) case) here, the scenario for the more general SU(N ) case is similar as well: simply replace M i , e, σ i with M (N ) i , e (N ) , t i , respectively.
MPO components for k-th local conserved quantity Q k for k > 2 are the upper triangle square matrices with bond dimension 3k − 1: where M i k and Θ k are the square matrices of size 3(k − 2), and the blank blocks are entirely filled with zero-operators.M i k is defined by: where the diagonal 3 by 3 block elements are all M i and the non-diagonal elements are all O 3,3 .Θ k is defined by: One can calculate the k-th local conserved quantity Q k by taking the product of Γ i k over all the sites and doing some boundary treatment.To be specific, for L ≥ k one can write the explicit form of the matrix product as: where m ≡ 3k − 2 and u k,L = (u 1 k,L , . . ., u ) is the row vector of 3(k − 1) dimension with its elements being the operator localized at the boundary.
has non-trivial action on the sites from the (L − k + 1)-th site to the L-th site, and v j k,1 has non-trivial action from the 1-st site to the k-th site.Q c k is the bulk term of the Q k , which can be more simply written with boundary vector: where 〈L| = (1, 0, . . ., 0) and |R〉 = (0, . . ., 0, 1) ⊤ are the boundary vectors of the same size with Γ i k .The local conserved quantity Q k under the periodic boundary condition is given with the boundary term The boundary term q B k is constructed from the operators that jump over the boundary, i.e. the linear combination of the operators that act across the boundary from the L-th site to the 1-st site.In the periodic boundary case considered here, the structure of the boundary terms is the same as that of the bulk term.
For the case that the Hamiltonian is defined on the infinite chain, i.e. the summations over sites L i=1 are all replaced with i∈ where is the set of all integers, the boundary term q B k becomes irrelevant and it is enough to consider only the bulk term.The MPO representation for the local charge Q ∞ k in the infinite chain is simply given by: where for the ordering of the operators in the product, we choose a convention that operators with lower site indices act first.When taking products of Γ i k , we can treat each building block as if it was a non-commutative scalar because the building blocks appearing in the MPO are conformable for multiplication.Therefore, in the following, we denote the zero-operator block matrix O n,m just simply as O.
We give the expressions up to Γ i 8 below: , where the lower triangular elements are all zero-operators and the O's in these expressions represent zero-operator block matrices, suitably shaped for the position of their respective blocks.
In the following section, we treat the building block elements of Γ i k such as M i and σ and e and O as a symbol, and the indices indicate the positions of these building block elements.For example, Γ i k 1,2 = σ i , and Γ i k 2,3 = M i for k ≥ 3, and Γ i 8 2,8 = C 2 e.We note that Γ i k is only involved with the usual Catalan number, while the expressions for the local charges themselves (4) are also involved with the generalized Catalan number.In this point, the MPO representation (18) is more simplified than the traditional bare expressions (4).
We briefly explain why the local charges can be represented simply by the MPO.The point is that basic components of the local charges (2) or ( 8) are the nested product of the Pauli matrices, and consequently, we can rewrite them by the products of M i like (14) or (16).We can then factorize them into the form of an MPO, as shown in (21).However, this does not apply straightforwardly to the anisotropic generalizations: basic components of the local charges for the XXZ or XYZ chains do not have this nested product structure [9,10], and they cannot be easily expressed in a matrix product form similar to ( 14) and (16).We leave the anisotropic generalization of our result to the future problem.

Boost operator
One may think we can obtain the recursive relation between Γ i k and Γ i k+1 using the Boost operator [30].In this subsection, we investigate the recursive way to obtain Γ i k .For simplicity, we consider the infinite chain here.Within this subsection, we denote the local charges in the infinite chain Q ∞ k simply by Q k .The boost operator in the infinite chain is defined by: which generate the recursive relation between Q k and Q k+1 : Strictly speaking, the local charge obtained from the boost operation is different from our Q k defined in (4); local charges have the freedom to add the lower-order charges, and the boostderived charge and our Q k exhibit distinct linear combinations of lower-order charges.Thus, Q k in ( 26) is slightly different from its original definition (4), and the corresponding MPO representation is also slightly different from (18).However, the bond dimension of the MPO for boost-derived charges seems to be still the same as that of (18).Within this subsection, we denote boost-derived charges by Q k and its MPO by Γ j k .The MPO representation for the boost operator is given by: where the bond dimension is 5, which is the same as that of the Hamiltonian, and the boost operator is written as: where we denote 〈L n | = (1, 0, . . ., 0) and |R〉 = (0, . . ., 0, 1) ⊤ with 0, . . ., 0 representing the sequence of 3n − 2 zeros.
We denote the local physical space corresponding to j-th site by h j ∼ = 2 , and denote the auxiualy space of Γ j k , corresponding to the bond dimension 3k − 1, by V a 1 ∼ = 3k−1 .Similary, we denote the auxiualy space of K j , corresponding to the bond dimension 5, by V a 2 ∼ = 5 .We define the Hilbert space W ≡ j∈ h j ⊗ V a 1 ⊗ V a 2 .We can treat Γ j k as an operator on W, acting non-trivially on h j and V a 1 , and K j as an operator on W, acting non-trivially on h j and V a 2 .Given this context, we denote Γ j k and K j as Γ j,a 1 k and K j,a 2 , respectively.With these notations and (24), we have: where , and a vector with the indices a 1(2) denotes a vector in V a 1 (2) .We can represent the boost recuresion relation (26) as: where we denote L ≡ ( a 〈L| , a 〈L| ), and R ≡ ( |R〉 a , − |R〉 a ) ⊤ , which are vectors of dimension 10(3k − 1), and where the elements of two-by-two matrix Γ j k+1 is the operator on W, and O is zero-operator on W, and the bond dimension of Γ j k+1 is 10(3k − 1).In this way, we can derive the MPO representation for Q k+1 as Γ j k+1 from the boost recursion relation (26).However, Γ j k+1 is different from Γ j k+1 : the bond dimension of Γ j k+1 is 10(3k − 1) whereas Γ j k+1 possesses a bond dimension of 3(k + 1) − 1, which is smaller than that of Γ j k+1 .Thus, Γ j k+1 is simpler than Γ j k+1 .Though it might be feasible to derive the expression for Γ j k+1 from Γ j k+1 through the compression of the MPO, addressing this is beyond our current study.

Proof of the MPO
In this section, we give the demonstration that the MPO introduced in the previous section actually reproduces the local conserved quantities of the spin-1/2 Heisenberg chain.

Catalan number identity
There is the useful identity of generalized Catalan numbers C n+m,n and usual Catalan numbers C n for the proof of the MPO: One can prove this recurrence relation (32) by relating C n+m,n to the numbers of monotone lattice paths from (x, y) = (0, 0) to (n, n + m) that never crosses the line y = x + m.
Figure 1: Graphical explanation of the proof of the recurrence relations of generalized Catalan numbers eq.(32).n = 11, m = 3 case is shown.C n+m,n equals the number of possible lattice paths from (0, 0) to (n + m, n) that never crosses the line y = x +m.If one divides a path P into (m+1) pieces P 0 , P 1 , . . ., P m by the point where P first reaches the line y = x + j, each P j can be related one by one to a possible lattice path from (0, 0) to (x j+1 − x j , x j+1 − x j ) that never reaches the line y = x.
Let Ω n+m,n represent the set of monotone lattice paths that meet this condition.For each P ∈ Ω n+m,n , let x j (P) be the smallest x coordinate of the point on the path where the path reaches the line y = x + j for the first time for 1 ≤ j ≤ m.We define x 0 (P) ≡ 0 and x m+1 (P) ≡ n.We define the subset of Ω n+m,n as follows: Ω n+m,n can be represented by the direct sum of Ω n 0 ,...,n m n+m,n : Let us denote P j as the part of a path P from (x j , x j ) to (x j+1 , x j+1 ), corresponding to the element of Ω n j ,n j , that is, to a monotone lattice path from (0, 0) to (n j , n j ) that never crosses the line y = x (see Fig. 1).Because the number of elements in Ω n j ,n j equals C n j , one gets where #[•] denotes the number of elements in a (finite) set.Using (34), we have Thus, we have proved (32).
After the completion of this work, we became aware of the conference proceedings [31], where the identity (32) is discussed.While they derive the identity analytically, we derive the identity combinatorially with the lattice path explanation.A formula similar to but distinct from ours has also been studied in [32], where the summation of indices traverses n j ≥ 1.
where the indices indicate the position of the building blocks, and γ p a p ,a p+1 represents the (a p , a p+1 )-th block element of Γ p 12 , and the second equality holds from the fact that Γ p 12 is an upper triangular matrix.In the summation in (37), σ 1 M 2 M 4 σ ⊤ 6 is included in the restricted summation of (37) as follows: where the other variables in (37) is fixed as a p = 13 for p ≥ 7, and therefore γ p a p ,a p+1 = I for p ≥ 7.
Each term in (38) corresponds to the "path" from "start" to "end" traversing through the matrix in Figure 2, which has the same structure as Γ i 12 .A "path" is defined as follows: first, we pick up the element in the first row, where only σ is non-zero2 .Thus we have to pick σ first, corresponding to γ 1  1,2 = σ 1 .Every time we pick an element, we then vertically descend within the same column till we reach the diagonal line, after which we move horizontally to the right within the same row and pick up an element on the row.The element chosen in our p-th pick corresponds to γ p a p ,a p+1 in equation (38).At the final (sixth) pick, we have to pick σ ⊤ , corresponding to γ 6 a 6 ,a 7 = σ ⊤  6 , and finish the procedure for this path.By considering all possible paths, we can compute equation (38).For calculating the coefficient of 12 , we have to pick up the point of σ, M , C n 1 e, M , C n 2 e, σ ⊤ in the matrix in Figure 2 in that order.There are three possible paths that generate σ 1 M 2 M 4 σ ⊤ 6 in Q c 12 , which are depicted by the teal, the bold blue, and the purple line in Figure 2, corresponding to (n 1 , n 2 ) = (0, 2), (1, 1), and (2, 0) respectively.These paths all satisfy n 1 + n 2 = 2.All the contribution to the coefficient is where we used the Catalan number identity (32).Thus, we have demonstrated that the coefficients of σ 1 M 2 M 4 σ ⊤ 6 in the MPO is actually C 3,2 , which is the same as the case of the local conserved quantity Q 12 .
Let us consider a more general situation, leaving the rigorous proof to the appendix A. We calculate the coefficient of the operator have to pick up the elements of C n 1 e, . . ., C n m e, corresponding to the m-holes.Every time we pick up C n e on the paths, there is a horizontal move of (2n + 2) columns to the right, and every 12 , which is one of components of F 4,2 .The three paths depicted by the teal, the bold blue, and the purple line contribute to the coefficients.For example, the bold blue line represents the path that picks σ, M , C 1 e, M , C 1 e, σ ⊤ in this order, and the contribution is To obtain the coefficients of F k−2(n+m),m in more general Q c k , we must take all the possible paths including C n 1 e, . . ., C n m e in this order that satisfy time we pick up M , there is a horizontal move of one column to the right.The number of the columns between σ and σ ⊤ of the matrix corresponding to Γ i k is k − 2, thus we have the relation and all the contributions to the coefficients become where we used the Catalan number identity (32), and we can see the coefficient of ( 4) is reproduced by the MPO.

Summary and Outlook
In this work, we show the matrix product operator(MPO) representation of the local conserved quantities for the spin-1/2 Heisenberg chain and its SU(N ) generalization.In terms of our MPO representation, the local conserved quantities are more simply written: the pattern in the expression of the MPO is more straightforward than that of the local conserved quantities themselves.Especially the coefficients appearing in the MPO are simpler than those in the bare expressions for the local charges (4), and the Catalan tree pattern of the local conserved quantities naturally appears from the product of the MPOs.This means the complexity of the local conserved quantities is folded into the product of the MPO.
The generalization of our result to the local charges in the spin-1/2 XYZ chain [9], which is the anisotropic generalization of the Heisenberg chain considered here is an interesting topic.It would also be worthwhile to verify how the mutual commutativity of the local conserved quantities can be explained with our MPO representation.
Our strategy may be useful to find the general expressions for the local charges in more general integrable quantum spin chains.To obtain the general expression for the local charges, one must identify the operator basis that constructs them and discern the regularity of their coefficients.When looking at the bare local charges, these tasks are very hard for interacting integrable spin chains.Hence, once one represents the lower-order charges with MPO, the patterns may become clearer and we may predict the general forms for the local charges via MPO representation.The local charges considered in this work have already been found for decades before.Therefore, it is interesting to investigate the MPO representation for local charges whose general forms have not been found.In this spirit, the investigation of local charges in open boundary cases is also interesting.In the periodic boundary condition, the boundary term of the local conserved quantity is trivial: they have the same structure as the bulk term, as discussed in Appendix B. However, the boundary term in the open boundary condition is non-trivial, and the expressions have yet to be elucidated even for the most famous spin-1/2 Heisenberg chain [33].Using the MPO representation, we might infer the pattern of the boundary terms in the open boundary case.

A Rigorous proof for bulk term
In appendix A, we give the rigorous proof of the bulk part of (23).
In the following, the matrix indices indicate the position of the building blocks, i.e., treating the elements of Γ i k such as M i and σ and e and O in Γ i k as a symbol.Calculating the product of the MPO explicitly, we have: where the second equality holds from the fact that Γ i k is an upper triangular matrix, and we denote the matrix element of Considering the first row and the last column of Γ p k is written by: and γ p a,a = O for 2 ≤ a ≤ k, the summation in (41) can be decomposed as: where the second summation on the first row is taken over a i+2 , . . ., a j−1 and the other a p 's are fixed by a 1 , . . ., The first term of (43) becomes First term of (43) = 1≤i≤L−k+1 where C n+m,m bulk denotes the set of clusters C = {i 1 , . . ., i n } of order n with m holes, satisfying 1 ≤ i 1 < i n ≤ L, and we can see F bulk k,0 is the bulk term of F k,0 .We denote the number of p's that satisfies b p > 1 in the variable b i+1 , . . ., b j−1 of the second summation of (43) by m.Then the second term of (43) becomes: Second term of (43) In the third equality, we used the following relation: , and e behave as identity operator: eM p = M p e = M p and σ i e = σ i .In the fourth equality, we replace the variable by c r = b i ′ r .If c r is odd for some r, we can see γ c r = 0. Thus the non-zero contribution comes from the case that all c r is even, and in the sixth equality, we rewrite c r by c r = 2d r .For k − w + m being an even integer, w and m have to satisfy the relation w = k − 2n − m with ∃n ∈ .In the seventh equality, we change the variable by d r − 1 → d r , and in the ninth equality, we used the Catalan number identity (32).F bulk k−2(n+m),m is the bulk term of F k−2(n+m),m .Therefore, with (45) and (46), we have proved Q c k becomes: and we can see this is actually the bulk term of Q k .

B Boundary treatment
We prove the boundary terms of the local conserved quantities Q k are given by q B k ≡ u k,L v ⊤ k,1 .We write the k-th matrix product operator from i-th site to j-th site ( j − i ≥ k) by: where u k, j and v k,i are the row vectors of dimension 3(k −1), whose elements are independent of i and j, respectively, and are constructed from the local operators that act on at most k adjacent sites between the ( j−k+1)-th site and j-th site and between the i-th site and (i+k−1)th site, respectively.u k, j (v k,i ) is independent of i( j).Q c k,[i: j] is constructed from the local operators that act on at most k adjacent sites between the i-th site and the j-th site.
We can decompose the products in the MPO for the k-th local conserved quantities as: and we have Now, we can see the term in Q c k, [1:L+N ] whcih acts across the L-site and the L + 1-site is u k,L v ⊤ k,L+1 .Consequently, the boundary term q B k in system size L for periodic boundary condition is given by u k,L v ⊤ k,1 .

C Recursion equation for Θ k
Θ k can be obtained from the following recursion equation: where E k is a square matrix of order 3(k − 2), defined by: where a and b (1 ≤ a, b ≤ k − 2) indicate the coordinate for the 3 by 3 block matrices in Θ k , and = {0, 1, 2, . ..} is the set of all natural numbers and C n = 2n n − 2n n−1 is n-th Catalan number, which is the special case of the generalized Catalan number C n,n = C n .Θ k can also be obtained from a simple recursion equation, which will be shown in the appendix C.

where γ b p is
defined by γ b p ≡ C b p /2−1 (0) for b p is even(odd), and the factor on the left hand side as a c-number, and factored out on the right hand side because γ is proportional to e for b i ′ r >1 ) a,b := e (b − a = 1) O 3,3 (otherwise) , (52) where a and b (1 ≤ a, b ≤ k − 2) indicate the coordinate for the 3 by 3 block matrices in E k .
and in the second equality we define b p ≡ a p+1 − a p for i < p < j and γ only on the difference a p+1 − a p .In the last equality, we decompose the summation to the case that all b p = 1 for all p and the other cases that at least one b p is greater than one.The explicit form of γ depends