{'fraction_non_alphanumeric': 0.03965496294496416, 'fraction_numerical': 0.011954805005467136, 'mean_word_length': 4.17035175879397, 'pattern_counts': {'":': 2, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Conserved charges in the quantum simulation of integrable spin chains

When simulating the time evolution of quantum many-body systems on a digital quantum computer, one faces the challenges of quantum noise and of the Trotter error due to time discretization. The Trotter error in integrable spin chains can be under control if the discrete time evolution preserves integrability. In this work we implement, on a real quantum computer and on classical simulators, the integrable Trotterization of the spin-1/2 Heisenberg XXX spin chain. We study how quantum noise affects the time evolution of several conserved charges, and observe the decay of the expectation values. We in addition study the early time behaviors of the time evolution, which can potentially be used to benchmark quantum devices and algorithms in the future. We also provide an efficient method to generate the conserved charges at higher orders.

When simulating the time evolution of quantum many-body systems on a digital quantum computer, one faces the challenges of quantum noise and of the Trotter error due to time discretization. The Trotter error in integrable spin chains can be under control if the discrete time evolution preserves integrability. In this work we implement, on a real quantum computer and on classical simulators, the integrable Trotterization of the spin-1/2 Heisenberg XXX spin chain. We study how quantum noise affects the time evolution of several conserved charges, and observe the decay of the expectation values. We in addition study the early time behaviors of the time evolution, which can potentially be used to benchmark quantum devices and algorithms in the future. We also provide an efficient method to generate the conserved charges at higher orders. One of the most promising applications of quantum computing is the simulation of quantum-mechanical systems in Nature with many degrees of freedom [1]. Future quantum computers are expected to provide powerful resources for simulations, which cannot be matched by their classical counterparts. Ultimately one would like to simulate the time evolution of realistic quantum field theories such as the Standard Model of particle physics, quantum chromodynamics and condensed-matter systems. At present, however, noise is one of the most significant obstacles in quantum computation. It is important to understand and reduce its effects on the digital quantum simulation of quantum many-body systems.
When simulating the time evolution of a quantum mechanical system with Hamiltonian H on a digital quantum computer, one often discretizes time t and approximate e −itH by U d , namely d repeated applications of a product of unitary operators U composed of elementary gates [2]. This approximation, called the Lie-Trotter-Suzuki decomposition or Trotterization, introduces a Trotter error. For given t one can make the Trotter error small by taking the depth d large. However, increasing d introduces more errors due to quantum noise. Thus, there is a competition between noise and the Trotter error.
Integrable spin chains can be Trotterized in a way that preserves integrability, called integrable Trotterization [3,4]. In this prescription, the original conserved charges including Hamiltonian are slightly deformed, but remain mutually commuting. Therefore the new conserved charges are free of the Trotter error, and exactly conserved in the discrete time evolution.
The study of integrable spin chains is important in its own right. They are among the most well-studied quantum many-body systems where many exact results are known, and often serve as useful stepping stones for studying more general non-integrable models. Integrable spin chains may also be regarded as simple lattice quantum field theories in 1+1 spacetime dimensions, where both time and space are discrete.
In this work, we implement the integrable Trotterization of the spin-1/2 Heisenberg XXX spin chain [3] with periodic boundary conditions, on real quantum computers and on classical simulators. We devise protocols to simulate the discrete time evolution and measure the expectation values of the conserved charges. While the Trotter error is under perfect control, we find that the effects of quantum noise are significant, forcing the charges to decay. We study the behaviors of the decay using classical simulators and the theory of quantum error channels.
We also suggest that the early-time behaviors of the charges may be used to benchmark quantum devices and algorithms in the future when the noise exists but has smaller effects. Compared with other benchmarking protocols such as randomized benchmarking [5][6][7][8][9], the (non-)conservation of charges will provide intuitive benchmarking measures in the setting directly related to the simulation of quantum field theories.
As a by-product, we obtain explicit expressions for low-degree higher conserved charges of the integrably-Trotterized XXX spin chain. It is known that these conserved charges can be recursively generated by the so-called boost relation. Since the conserved charges become quickly complicated as the degree increases, we invent an efficient way to solve the recursion relations and implement it as a classical computer (Mathematica) code.
The simulation of integrable spin chains on digital quantum computers has been intensively studied these days [10][11][12][13][14][15][16], and even Lindblad superoperators can be integrable [17,18]. 1 To the best of our knowledge our work is the first study of integrable Trotterization on a real quantum device. The rest of this paper is organized as follows. In section II we review the Trotterized XXX spin chain and discuss recursion related for the conserved charges. In section III we propose quantum circuits for measuring the expectation values of the conserved charges. In section IV we present the results of our simulations (both on simulators and on a real device of IBM Quantum). In section V we comment on the implications of the simulation results. We conclude in section VI with summary and discussions. We also include several appendices for technical materials. Some of our codes and the data presented in this paper are publicly available, as explained in Appendix F.

II. CONSERVED CHARGES IN TROTTERIZED XXX SPIN CHAINS
A. Integrable Trotterization of the XXX spin chain As a prototypical example of an integrable spin chain, we consider the standard spin-1/2 XXX spin chain defined by the Hamiltonian where σ j = (X j , Y j , Z j ) denotes the vector of Pauli operators (σ x , σ y , σ z ) on site j. We take N to be even and impose the periodic boundary condition σ N +1 = σ 1 . This model can be exactly solved by the algebraic Bethe ansatz, which can be formulated in terms of the R-matrix and the transfer matrix (See [24,25] for reviews).
For quantum simulations, we need to discretize the time evolution operator exp(−iHt) into a sequence of discrete time evolutions. In [3], a deformation of the XXX model was introduced, which we call the Trotterized XXX model. The Trotterized XXX model is a special case of the inhomogeneous XXX model, where the inhomogeneities are controlled by the deformation parameter δ. In this model, we can discretize the time evolution while preserving integrability.
Let us write the R-matrix of the original XXX model as R ij (λ) = P ijŘij (λ), wherě and P ij = (1 + σ i · σ j )/2 acts as a permutation of the qubits at sites i and j. We introduced a formal variable λ, called the spectral parameter. The discrete time evolution is generated by d (∼depth) repeated actions of the unitary operator where the deformation parameter δ is regarded as a Trotterized time step, so that the total time for the evolution is t = −(δ/J)d. Note that the ordering matters in (3): theŘ-matrices within each bracket mutually commute, but those from the two different brackets may not commute with each other. The continuous time evolution of the original XXX model can be studied by the relation for large d and fixed J, t. The second term represents the Trotter error. Below we take δ a general real parameter. This unitary operator can be decomposed as 2 where the transfer matrix is defined by Here the product is taken in the ascending order from right to left, and the trace is taken in an auxiliary space (ancilla). The Trotterized XXX model is integrable because the transfer matrices with arbitrary and different values of the spectral parameter mutually commute: [T (λ), T (µ)] = 0.

B. Recursion relations for conserved charges
The mutual commutativity of the family of transfer matrices T (λ) implies conservation of the charges under the Trotterized time evolution: Throughout the paper, the length N of the spin chain must be greater than 2n + 1 when measuring Q ± n , because the conserved charges above this bound are not generated by the simple recursion relation given below.
Let us introduce the precise expressions of the charges Q ± n , which will be used in the quantum simulation. For our purposes it is useful to start with the explicit formulas for Q ± 1 and Q ± 2 as follows. We first define the charges in terms of the densities: which contain Pauli matrices acting on at most 2n + 1 sites. As shown in [3], q and where σ 1 · (σ 2 × σ 3 × · · · × σ −1 × σ ) := σ 1 · (σ 2 × (σ 3 × (· · · × (σ −1 × σ ) · · · ))). One finds that Q ± 1 reduce to the original XXX Hamiltonian (1) at δ = 0. As can be seen from these examples, the charges and densities are defined so that there is no term proportional to the identity.
We compute the higher charges by the recursion relations where the degree of the new charge should not exceed the bound (7). The boost operator B is defined by [3] One can check that the boost operation is consistent with (10), (11). We wrote computer programs to compute the charges Q ± n up to n = 6. More details on the recursion relation can be found in Appendix A. For Q + 3 we obtain a simple expression (A10) for its density.
We will measure Q + n and the difference between Q + n (δ) and Q − n (δ). Since Q + n (δ = 0) = Q − n (δ = 0), 4 we normalize the difference as which are also the conserved charges that are polynomials in δ with integer coefficients.

III. QUANTUM CIRCUITS FOR CONSERVED-CHARGE MEASUREMENTS
We are now ready to discuss our quantum circuits. Our quantum circuit is divided into three parts: initialization, time evolution and measurement: The first step is to initialize the state before the dynamical time evolution. For our discussions of conserved charges, we can choose any state as the initial state |ψ 0 : the conserved charges will be time-independent under the noiseless time evolution regardless of the choice of the initial state. In our work we consider three types of the initial state: • The ground state |00 . . . 0 in the computational basis.
• An eigenstate |s 1 s 2 . . . s N P1P2...P N for the operator P 1 ⊗ P 2 ⊗ . . . ⊗ P N (P j ∈ {X, Y, Z}) with eigenvalues (−1) s1+s2+···+s N , with s j ∈ {0, 1}. 5 The third type is the most general, and it includes the first two as special cases. The state initialization part of the circuit (16) consists entirely of one-qubit gates, namely X sj j followed by H j (if P j = X j ), S j H j (if P j = Y j ), or none (if P j = Z j ). 6 See FIG. 1(a).

B. Time evolution
The second step is to consider the time evolution part of the circuit (16). The Trotterized time evolution is given by U(δ) in (3), which can be expressed as 5 For 1 ≤ j ≤ N − 1, the j-th qubit here is also the j-th qubit in our codes. The N -th qubit here is the 0-th qubit in the codes. 6 In the computational basis, where a rectangle represents an R-matrix. We split the R-matrix (2) into a phase and two unitary operators: We will ignore the overall phase e −i α 2 because it does not affect measurements. 7 The two unitary operators in (18) can be realized in terms of elementary gates as follows: The rotation gate is defined by R Z (α) = e −iαZ/2 , in agreement with RZGate in Qiskit. The time evolution part of the circuit (16) consists of d repetitions of (19) and (20). We note that pairs of CNOT gates cancel when (19) and (20) are combined. 8

C. Measurement protocol via Pauli words
Let us next discuss the measurement of conserved charges. Suppose that we have a quantum state described by a density matrix ρ. We then wish to evaluate the expectation value Q = tr(ρQ) of the charge Q, as well as its statistical uncertainty.
We review below a measurement protocol via "Pauli words". See, for example, [30]. Recall that the charge density (8), (9) is a sum over the product of Pauli matrices, P, P , . . . . Our idea is to measure another Pauli operator (Pauli word) W which contains P, P as a subset. The measurement of W tells us the expectation values of P, P , . . . all at once. To compute the charge density, we need to choose an appropriate set of Pauli words and repeatedly run the circuit for all the Pauli words in the set.
According to the recursion relation (12) a general conserved charge Q can be expanded as The problem thus reduces to the evaluation of the expectation values P = tr(ρP ) for all the Pauli operators P that appear in the expansion (21) of a given conserved charge Q. Let us consider a word W = L 1 . . . L N whose letters L j (j = 1, . . . , N ) are one of X, Y , and Z (but not I). We use such a "Pauli word" to specify a measurement basis (W -basis) which diagonalizes the product operator L 1 ⊗ . . . ⊗ L N [30]. For each qubit, the letter Z corresponds to the computational basis {|0 , |1 }, X to {|+ = H|0 , |− = H|1 }, and Y to {SH|0 , SH|1 }, with the eigenvalues {+1, −1} respectively. We say that a Pauli operator P = P 1 ⊗ . . . ⊗ P N ∈ {I, X, Y, Z} ⊗N is contained in the Pauli word W (and write or W ⊃ P ) if P j = L j for all j such that P j = I. 9 For given Q, we need to find a set S of Pauli words W such that all the Pauli operators P that appear in (21) are contained in some W ∈ S. Given such S, measurements along the bases specified by all words in S are sufficient to evaluate Q . Note that we allow each P to be contained in more than one word. To generate an appropriate S, we make use of a code (data_acquisition_shadow.py) provided by H.-Y. Huang [31]. 10 The standard computational basis measurement can be converted to the W -basis measurement by including the corresponding one-qubit gates. These gates are the inverse of those considered in Section III A. FIG. 1(b). Suppose that we execute the circuit n W times in the basis specified by W ∈ S. In each execution, the measurement outcome b j ∈ {0, 1} for the j-th qubit means that the eigenvalue (−1) bj of the single-qubit Pauli operator L j ∈ 7 Usually the overall scalar factor of the R-matrix is neglected in quantum circuits. Thus, our simulation corresponds to any of the integrable spin chains that have the same R-matrix as the Heisenberg XXX spin chain up to a scalar factor. For example, the bi-local deformation modifies the scalar factor, which redefines the conserved charges [29]. 8 It has recently been found that three CNOT gates are sufficient to construct this R-matrix as a quantum gate [23]. 9 For example, both X ⊗ I ⊗ Z ⊗ I = X 1 Z 3 and X ⊗ X ⊗ I ⊗ I = X 1 X 2 are contained in XXZY . 10 The code was developed to implement classical shadow [32] and derandomization [33] schemes. Note that it is not guaranteed that this code generates the most efficient Pauli word set S for any observables.
{X, Y, Z} is observed. This defines a bit string b = b 1 . . . b N ∈ {0, 1} N . We perform n W measurements in the W -basis. The result for the j-th qubit in the i-th measurement is denoted by b (W,i) j . As the estimator of Q we take where 11 Here n P := W ⊃P n W and w j (P ) = P j for P = P 1 ⊗ . . . ⊗ P N . By grouping together P, P , . . . in one word W , the measurement results for P, P , . . . are no longer statistically independent. In Appendix B, we obtain an unbiased estimator of the variance in Q est by taking this effect into account. The unbiased estimator is where n P,P := W ⊃P,P n W and We use its square-root s Q to estimate the statistical uncertainty.

IV. SIMULATION RESULTS
In this section we present the results of our simulations on a real device as well as on classical simulators. In the ideal simulation, the conserved charges Q + n and Q dif n remain constant under time evolution. On real NISQ devices, however, the charges do not stay constant due to noise. To analyze such effects, we computed the time evolution of conserved charges and performed quantum state tomography using classical simulators.
We used Qiskit -IBM's open-source software development kit for quantum computers [34] -to create circuits and execute jobs on the real quantum device ibm kawasaki of IBM. We also used Qiskit to perform classically-emulated quantum computation. Throughout this section, the error bar represents the statistical uncertainty s Q discussed in Section III C, and the value of α in (18)  In Appendix C, we provide further technical details of the simulations, and examine the aforementioned systematic errors.

B. Classically-emulated quantum computation with large noise rates
In this subsection we discuss the classically emulated quantum computation with quantum errors. We consider two kinds of noise: 1) depolarizing errors and 2) combined amplitude-and-phase damping errors. The errors are with the following initial states: |0000 Y ZXX for 4 sites, |000000 Y ZXY ZX for 6 sites, |00000000 Y ZXY ZXY X for 8 sites, |0000000000 Y ZXY ZXY ZXX for 10 sites, and |000000000000 Y ZXY ZXY ZXY ZX for 12 sites. In both (a) and (b), error bars representing statistical uncertainties are small and hidden by markers, and the dashed curves indicate the fitted functions c1e −γd + c2 with constants c1, c2, and γ. The dotted horizontal lines indicate the (rescaled) values of the charges in the noiseless case. Each data set is rescaled by a factor C, chosen conveniently for better visibility.
probabilistically inserted after each gate operation. Recall that our circuit is mostly composed of the time evolution part, which is a product of R-matrices between adjacent qubits as in (17). The R-matrix circuit has two kinds of one-qubit gates H and R Z (±α), together with the two-qubit gates CX. Thus, the errors quickly and uniformly spread over the entire circuit during the time evolution.
The state preparation and measurement parts also have one-qubit gates (X, H, and S † ). The conserved charges Q + n and Q dif n defined in (15) are measured.

Simulations with depolarizing errors
Let us first consider the noise model with depolarizing errors. A single-qubit depolarizing channel with a rate 0 ≤ p ≤ 1 acts on the density matrix ρ as 12 with Kraus operators In the simulation, after each one-qubit gate, we apply Φ depo with p = 0.0013 on the qubit. After each two-qubit (CNOT) gate, we apply the two-qubit channel Φ depo ⊗ Φ depo with p = 0.013. In FIG. 3(a) we show the results of classically emulated quantum simulations with depolarizing errors. We see that the expectation value Q of the conserved charge

Simulations with amplitude-and-phase damping errors
We define the combined amplitude-and-phase damping error, with rates λ a and λ p , on a single qubit as 13 where † denotes hermitian conjugation, and the Kraus operators are given as In our noise model we choose the values λ a = λ p = 0.018. After each two-qubit gate, we apply Φ damp ⊗ Φ damp . For simplicity, we let one-qubit gates be noise-free. In FIG. 3(b) we show the results of classically emulated quantum simulations with amplitude and phase damping errors. The expectation value Q of a conserved charge Q ∈ {Q + 1 , Q + 2 } is approximated by c 1 e −γt + c 2 for some constants c 1 and c 2 . The constant c 2 is non-zero for some charges unlike the case of depolarizing errors, and it does not depend on the initial state. These conclusions also hold for the results for Q ∈ {Q dif 1 , Q dif 2 , Q + 3 }, which we omit in FIG. 3(b) to avoid clutter. We conclude that the final state is not the completely mixed state.
We will strengthen our understanding by quantum state tomography in Section IV C and the spectrum analysis of the quantum channels in Section V A.

C. Quantum state tomography
We have observed that the expectation values of all conserved charges at late times approach a constant value independent of the choice of the initial state. This indicates that the density matrix of the underlying systems gets closer to each other. We perform another noisy simulation for the quantum state tomography (QST), which is the procedure of reconstructing the density matrix from measurements. In QST, we perform repeated measurements in some basis, and reconstruct the density matrix by solving an optimization problem. We employed qiskit-experiments framework to perform this procedure. 14 In particular, we chose the Pauli measurement basis and applied the fitter fitters.cvxpy gaussian lstsq to solve the optimization problem. We performed S = 40, 000 shots for each element of the Pauli measurement basis. 15 Since QST requires large computational resources, we performed a simulation on the system with 6 sites. Recently, another efficient method for the full quantum state tomography has been proposed in [36]. See Appendix D for the details of the quantum state tomography algorithm. FIG. 4(a) shows the self-fidelity for three choices of the initial state, namely the fidelity between the ideal density matrix constructed from the gate set at d = 0, and the density matrix after Trotter time evolution.
where the state fidelity F is defined by F (ρ, σ) = (tr √ ρσ √ ρ) 2 . There are small bumps around d ∼ 7 in the figure, which seem to be related to our choice of α = 0.3. FIG. 4(b) shows the fidelity between a pair of the density matrices at step d, evolved from different initial states, The fidelity for all pairs becomes close to 1 at a late time, and we expect that the final fidelity approaches 1 as we increase the number of shots. This fact suggests that all states approach a fixed point of the noisy time-evolution operator. Clearly, when the density matrices are identical, the expectation values of any observables must also be identical.

D. Classically emulated quantum simulations with small noise rates
According to our simulation results on a real device in Section IV A, the conserved charges rapidly decay in the first few Trotter steps. We expect that this initial decay will be slower in less noisy devices in the FTQC era, and deviates from the initial value by small corrections. In this case we expect that we can quantify the initial change of the charge Q by a linear function for a small number of Trotter steps:   . 3(a). Namely, after each one-qubit gate we apply Φ depo with p = 0.00013 on the qubit, and after each two-qubit gate we apply Φ depo with p = 0.0013 on each of the two qubits. We fit the data by dotted straight lines.
(In practice it is useful to choose the initial state and the charge so that Q d=0 = 0, in which case it is easier to extract the linear behavior (32) from numerical results.) In FIG. 5 we present classical simulation results for small noise rates. We can conclude that the values of β are in general different for different choices of charges and for different initial states. We also find that the values of β for Q diff 1 and Q + 2 are larger than those for Q + 1 . While data is limited and a more thorough analysis is necessary, one is tempted to conjecture that one obtains larger values of β by considering Q ± n for larger values of n, and also by replacing Q ± n by Q diff n .

V. ANALYSIS OF THE SIMULATION RESULTS
We now discuss the implications of the simulation results. Since the available computational resources are not limitless, some of our results also include significant statistical uncertainties. Still, our simulation results have indicated some clear patterns already.
In our simulations we have observed two characteristic features at different Trotter steps: • "Late-time behavior": at large Trotter steps the value of Q ± n approaches a constant value.
• "Early-time behavior": at small Trotter steps we have a rapid decay of Q ± n .

A. Exponential decays at late times
In our simulations with large noise (both on simulators and on a real device) we find that the expectation value of a conserved charge Q of the integrable model decays exponentially and approaches constant values: where γ > 0 and c 1 , c 2 ∈ R. Moreover, the asymptotic value c 2 is independent of the choice of the initial state ρ 0 . Such a pattern is clear especially in our noisy simulations on classical emulators shown in FIG. 3.
We can explain this behavior as follows. In our quantum simulation of the integrable circuit, we prepare an N -qubit density matrix ρ 0 ∼ |ψ 0 ψ 0 |, and repeatedly apply a quantum channel Φ (i.e. a completely positive trace-preserving (CPTP) map) corresponding to a single Trotter step, so that the expectation value after d Trotter steps where Φ d is the d-th power of Φ. From the theory of quantum channels (see e.g. [37,38]), it is known that all the eigenvalues {λ n } of Φ must have absolute values less than or equal to 1, and that Φ has λ 0 = 1 as an eigenvalue (i.e. Φ has a fixed point). Denoting the fixed-point density matrix by ρ * , we have the asymptotic behavior (33) with If we assume that there is no degeneracy in the eigenvalue λ 0 (as one expects in generic noisy situations), 16 then ρ * is determined by the quantum channel only and is independent of the choice of the initial state.
The asymptotic density matrix ρ * depends heavily on the properties of the noise. For the time evolution Φ with depolarizing errors discussed in IV B 1, it is easy to see that Φ has the completely-mixed state ρ * = (1/2 N )I ⊗N as a fixed point. 17 This gives c 2 = 0 for traceless conserved charges Q, which is consistent with the behaviors of Q in FIG. 3. The situation is different for the time evolution with the amplitude-and-phase damping errors discussed in IV B 2, 18 where the asymptotic values of the charges do not necessarily vanish. Although statistical uncertainties are larger, our simulation results from the real quantum device are consistent with the hypothesis that the asymptotic values of Q + 1 and Q dif 1 both vanish. In all our simulations, the fact that the asymptotic values do not depend on the choice of the initial state indicates that the quantum channel has a unique fixed point.
Whether the fixed point of Φ is genuinely unique is an important question. To clarify this point, we calculated the spectrum of the quantum channel Φ numerically (in a classical computer) on the system with 4 sites. The details of this computation are explained in Appendix E. FIG.6 illustrates the eigenvalue distribution of the quantum channel Φ of noiseless and noisy (depolarizing error with the error rate p = 0.018) time evolution. 19 The eigenvalues of the unitary channel are distributed on the unit circle. In contrast, all the eigenvalues of the noisy channel are inside the unit circle, except a single eigenvalue with the value 1. As we increase the Trotter step in the noisy time evolution, the eigenvalues except for 1 move more deeply toward the center.
Let us make more quantitative analysis on the value of the decay rate γ in (33). One expects that the decay under the time evolution is governed by the eigenvalue λ = λ 1 with the largest among |λ| < 1, and that the late-time decay rate γ (inverse of the relaxation time) of the charge should equal − ln |λ 1 |. 20 While this relation should hold in idealized situations (at very late times and for λ 1 isolated from other eigenvalues), we find that it does not strictly hold for our simulation results. Indeed, the values of γ obtained by fitting the data in FIG. 3(a) (36) are comparable with but somewhat larger than − ln |λ 1 | ∼ 0.20, and depend on the type of the charge. This difference appears because we obtained γ by fitting the data for all steps d, and also because as shown in FIG. 6 we have many eigenvalues λ whose absolute values are smaller but close to |λ 1 | -this means that we will have a sum of exponentials with comparable sizes.

B. Towards benchmarking protocols from early-time decays
As we discussed in Section IV D, we expect that future real-device simulations with small noise rates will exhibit the linear behavior (32) for the initial change of the charge Q. The parameter β in (32) quantifies how "noisy" the quantum device is. It is therefore natural to propose to use β as a benchmark for quantum devices in the future.
While there are already many benchmarking proposals in the literature, there are some motivations for proposing a new benchmark. First, since we are discussing spin chains this benchmark could be more suitable for evaluating the performance of quantum many-body simulations. Second, integrable models exhibit a characteristic complexity growth different from chaotic or free models [43]. In this sense, our benchmark is conceptually distinct from the randomized benchmarks related to chaotic time evolution [5][6][7][8][9]. Third, the higher conserved charges {Q ± n } can be regarded as a large number of benchmark indicators, which provides refined data about the precision of real devices.
Note that exactly-known values of some particular quantities (e.g. energy as an eigenvalue of the Hamiltonian) have often been used in calibrating quantum devices, and there are also papers on integrable-model benchmarks in the NISQ era [10,44,45]. However, the essence of integrability is not only about the solvability of the spectrum of the Hamiltonian, but also the existence of as many conserved charges as the degrees of freedom. It is therefore natural to take advantage of all such conserved charges in the benchmarking.
Let us make some preliminary comments on the benchmarking protocols. First, we should keep in mind that the parameter β depends on many parameters/data: the choice of the initial state |ψ 0 , the size δ of the Trotter step, the choice of the charge Q (i.e. Q = Q ± n or Q diff n for some positive integer n), the number of qubits (width w), and the number of Trotterized time evolution (depth d). While δ is a parameter independent from others, for the purpose of simulating the real-time dynamics of the spin chain it is natural to fix the product d · δ to a constant value.
We can then propose the following procedure: 21 • Fix the choice of the w-qubit initial state |ψ 0 , as well as the threshold value β (say 1%).
• Run the depth-d quantum circuit that approximately prepares the state U(δ) d |ψ 0 .
• Compute Q by (many copies of) the state obtained above, and estimate the values of β in (32).
• We declare pass/fail depending on whether or not β is smaller/larger than β . We obtain a multi-dimensional plot by collecting the results for different values of w, d, n and for different initial conditions.
Compared with the "volumetric benchmarks" [46] with width w and depth d, here we have the third direction given by the order n of the conserved charge, and the fourth direction given by the choice of the initial state. The charges Q n are more non-local as we increase n and could thus in principle be better at capturing more non-localized errors than the usual Hamiltonian. It remains to be seen if this new direction is of any use for benchmarking purposes.

VI. CONCLUSION AND DISCUSSION
In this work we performed the digital quantum simulation of the integrably-Trotterized spin-1/2 Heisenberg XXX spin chain with periodic boundary conditions on a real quantum computer and on classical simulators. We designed and implemented the protocols to simulate the discrete time evolution and measurements. We computed the expectation values of the conserved charges together with their statistical uncertainties. While integrability guarantees that the conserved charges remain constant under the time evolution, we found that quantum noise forces the charges to decay. We analyzed noise models and performed quantum state tomography on classical simulators, and confirmed that the decay of the conserved charges is a typical behavior of noisy systems. Further, we demonstrated that the late-time behavior of the time evolution can be understood in terms of the spectral properties of the quantum channels. We also proposed that the decay rate β for the early-time behavior (32) in weakly noisy simulations can be used to benchmark quantum devices and algorithms in the future. This benchmark setup is directly relevant to the simulation of quantum field theories.
We made modest theoretical progress regarding the integrable Trotterization of the XXX spin chain. By rewriting the boost relations that recursively generate the conserved charges, we provided an efficient method to derive explicit expressions of the conserved charges using a classical computer program. Based on these results, we obtained a simple expression (A10) for the density of Q + 3 , as well as explicit expressions for higher charges up to Q + 6 . It would be interesting to find a closed-form expression for all Q + n , by generalizing the known results for the original XXX spin chain [47].
While we focused on the spin-1/2 Heisenberg XXX spin chain in this paper, it is straightforward to generalize our method to other integrable models. In Section IV D we found that the effects of noise are different for different charges. This may be due to the wider interaction range for the density of Q + n with higher n. Similarly, it is natural to ask if the effects of noise are different for different integrable models. Turning this around, we may hope that a systematic analysis of noises in quantum devices could lead to new classifications of quantum integrable models.
There are several points where our simulation methods may be improved. While we made sufficiently many measurements (shots) for the simulation on ibm kawasaki, we should be able to make our protocols and job management more efficient, thereby scaling up the qubit size N and Trotter step d. This will allow us to compute charges higher than Q + 1 and Q dif 1 that we computed on a real device. As a related issue, it is desirable to have a more quantitative understanding of the number of necessary measurement bases (Pauli words) needed to compute higher charges. To alleviate the systematic errors noted in Section IV A and discussed in Appendix C, it is desirable to have available a hardware-level technique such as the Floquet calibration used in [21]. We expect that error mitigation techniques can alleviate the effects of noise (see [48] for a recent review). It would be interesting to see if the error-mitigated computations can make the charges better conserved. This will pave the way toward studying quantum quench via digital quantum simulation because, in the continuum limit δ → 0, the time evolution with an initial state such as (the zero-momentum version of) the Néel state can be regarded as quantum quench.  In this appendix, we study in more detail the recursion relations for the conserved charges Q ± n discussed in II B.
To use the relation (A1), we assume that the densities take the form where the indices µ 1 , . . . , µ 2n+1 are summed over the values 0, 1, 2, 3, and σ 0 a denotes the identity I on site a. Since the density is a polynomial in δ, we define the coefficients in front of δ m as These coefficients can be read off by where the trace is taken over 2n + 1 qubits.

Appendix B: Formulas for statistical uncertainties
Following definitions, it is straightforward to check that the expectation value E[ Q est ] (with respect to the probability distribution defined by many copies of the quantum state ρ) of the estimator Q est given in (22) indeed equals Q ≡ tr(ρQ).
In practice, we also need to estimate the statistical uncertainty in the estimator Q est . The estimator is computed from the outcomes of many measurements in various measurement bases. If we repeat this set of measurements many times, the value of the estimator itself fluctuates according to a probability distribution whose variance is, by definition where (B2) The summand in (B1) vanishes unless (W, i) = (W , i ). By expanding the product we obtain c Q,P c Q,P n P n P (E [Π P,W,i Π P ,W,i ] − P P ) = P,P n P,P n P n P c Q,P c Q,P ( P P − P P ) , where n P,P := W ⊃P,P n W . If the quantum state ρ is known explicitly, the expression (B3) can be easily computed. But we should also be able to estimate the variance from measurement results. To find a good estimator, let us define S P,P := 1 n P,P W ⊃P,P where P P est := 1 n P,P W ⊃P,P Π P,W,i .
Writing Π for Π P,W,i and Π for Π P ,W,i for brevity, we have where we have used the relation to obtain the last line. Using the relation once again and evaluating the expectation values, we obtain E[S P,P ] = 1 − 1 n P,P ( P P − P P ) .
This shows that (24) is an unbiased estimator of the variance (B3): with tr(P Q) = 2 N δ P Q . We should set a P = 1 for P = I ⊗N to solve the unit-trace condition.
Suppose that the measurement result of the observable P is m P , which may contain various statistical and systematic errors. We consider the following weighted log-likelihood function, This function L log has the global minimum at ρ * = 1 2 N P m P P , which is called the linear inversion solution. This ρ * is a Hermitian matrix, but not always positive semi-definite for generic {m P }. We have to solve a convex optimization problem numerically to find a consistent density matrix.
If the weights are trivial, w P = 1 (∀P ), the log-likelihood function becomes 24 One can obtain the global minimum of this function in the space of density matrices by the PSD (positive semi-definite) rescaling [51].
The measurements for a mixed state can be described by a set of orthogonal projections which is called a projection-valued measure (PVM). For example, the PVM in the computational basis is given by the binary representation of i = 0, 1, . . . , 2 N − 1 as Each element Π i is the Kronecker product of the i-th eigenvector of Z ⊗N and its dual vector. The measurements of {Π i } gives 2 N integers, 24 Use the identities m P = tr(ρ * P ) and ρ = 1 2 N P tr(ρP )P .
where S is the total number of shots. The symbol means that we estimate the LHS (a number) from the RHS (an integer-valued random variable). Since a general density matrix (D1) has 4 N −1 real parameters, we need more PVM's for quantum state tomography. We use the so-called Pauli measurement basis, which consists of 3 N PVM's, where we assumed that the number of shots is the same for each measurement observable P (k) . From the measurement results, we define the weighted log-likelihood function and look for its (local) minimum in the space of density matrices. Note that this POVM basis {Π (k) i } does not satisfy the orthogonality condition tr(P Q) = 2 N δ P Q . Even if all the weights are trial w k,i = 1 , it is not guaranteed that the PSD rescaling gives the global minimum of the log-likelihood function in the space of density matrices.
It is straightforward to construct the quantum circuits for the Pauli measurement basis (D7), because (X, Y, Z) for each qubit can be measured as in FIG. 1(b). However, the log-likelihood function (D9) contains 6 N parameters whereas the density matrix contains (4 N −1) real parameters as in (D1). The Pauli measurement basis is over-complete for quantum state tomography, and the diagonal elements of ρ are usually more often measured than off-diagonal elements.

Appendix E: Details of spectral analysis
For the spectral analysis in Section V A, it is useful to invoke the Choi-Jamio lkowski isomorphism (channel-state duality) and represent the quantum channel as an operator acting on the states |ρ representing the density matrix ρ [52,53].
The density matrix ρ is a positive operator H → H acting on the Hilbert space H, which can be represented in an orthonormal basis {|α } of the Hilbert space as ρ = α,β ρ αβ |α ⊗ β|. (E1) The vectorized state |ρ of ρ can then be defined as For our purposes, recall that the circuit for the noiseless single-step time evolution in Section 3.1 is given by 25 Such a P (k) is called a Pauli word in Section III C.
We replace each unitary gate U with an operator U = U ⊗ U * , e.g.
Let us include the effect of noises represented by a quantum channel Φ err acting on a single qubit. In terms of the Kraus operators {D k }, this can be represented as, e.g.
The noisy time evolution Φ is defined by inserting the quantum channel Φ err into the noiseless time evolution (E4): The noisy time evolution at time d is obtained by repeated action of the operator Φ and the expectation value Q d of the charge Q is given by We studied the eigensystem of two types of Φ, one with the depolarizing error and the other with the damping errors. Both errors are introduced in section IV B, and the corresponding Kraus operators are defined in eqs. 27,29, respectively. We confirmed that the analytic results calculated by the above algorithm agree with those from the quantum simulations using the noise model on the Qiskit for small values of N .