Scattering in the Ising Model Using Quantum Lanczos Algorithm

Time evolution and scattering simulation in phenomenological models are of great interest for testing and validating the potential for near-term quantum computers to simulate quantum field theories. Here, we simulate one-particle propagation and two-particle scattering in the one-dimensional transverse Ising model for 3 and 4 spatial sites with periodic boundary conditions on a quantum computer. We use the quantum Lanczos algorithm to obtain all energy levels and corresponding eigenstates of the system. We simplify the quantum computation by taking advantage of the symmetries of the system. These results enable us to compute one- and two-particle transition amplitudes, particle numbers for spatial sites, and the transverse magnetization as functions of time. The quantum circuits were executed on IBM 5-qubit superconducting hardware. The experimental results with readout error mitigation are in very good agreement with the values obtained using exact diagonalization.


I. INTRODUCTION
The Ising model is a quintessential spin system within which one can simulate and study many-body interactions. The model allows for simulating spin-spin physics and the calculation of properties such as magnetization and spin-frustration. For instance, the Ising model is of critical importance in the study of high T c superconductors since it allows one to study the electrical transport properties near a quantum critical point (Ising-nematic) which helps one understand strong electronic interactions in these systems [1]. To obtain information about the non-equilibrium dynamics of isolated many-body systems, the time evolution of the transverse magnetization as well as the entanglement entropy of the evolved states have been evaluated to study the domain wall melting in the ferromagnetic phase of transverse Ising chains [2]. The Ising model also serves as a useful arena for the study of more complex quantum field theories on a lattice. For example, scattering in a spin system on a lattice holds many parallels with scattering between particles in high energy physics experiments [3][4][5][6][7].
In a different perspective, the Ising model itself is used as a generic quantum computer model for the adiabatic quantum computers and quantum annealers [8]. Therefore, the dynamics and the correlations of the quantum entanglement between large number of spins are of interest related to the operation of quantum annealers [9].
Computing scattering amplitudes, transition rates, and other physical quantities involving quantum fields are hard tasks for classical computers. Quantum computers promise exponential speedup, however the approach with quantum simulators often revolves around the computation of real-time evolution based on Trotterization which is of limited utility on NISQ (noisy intermediatescale quantum [10]) hardware [11]. Previous studies have simulated real-time dynamics of interactions [3,[5][6][7] and evolution of disordered Hamiltonians [12] with this method. In this type of simulation, the number of gates grows linearly with the system size and the number of Trotter steps. Therefore, the noise in the system grows as the system size grows.
As an alternative perspective for NISQ devices, the variational quantum simulation of real time, imaginary time, and generalized time evolution of quantum systems have also been studied [13,14]. The introduction of variational quantum simulation for studying real-time dynamics of quantum systems and comparison to Trotterization method was first done in [15] where it was claimed to provide an improvement over the Trotterization method. Specifically in Ref. [14], the authors introduced a variational quantum simulation of open system dynamics and numerically tested the algorithm with a 6qubit 2D transverse-field Ising model under dissipation. A comprehensive review of these variational algorithms can be found in [16].
Others have simulated the Ising model both variationally [17] and via direct diagonalization within the quantum circuit [18].
Here, we use the Quantum Lanczos (QLanczos) algorithm [19] to calculate transition probabilities and scat-2 tering amplitudes in the one-dimensional transverse Ising model with periodic boundary conditions. We use the quantum imaginary-time evolution algorithm (QITE) to provide a basis for the Hilbert (Krylov) space employed in the QLanczos algorithm. We tune the QITE step size, and thus the total noise in the circuit, by using a hybrid quantum-classical approach to the algorithm. Using this technique we also compute occupation numbers and the transverse magnetization.
Extending the results of this study for the use of Ising model as a generic model would be an interesting research but we will leave this as a future study.
The hybrid quantum-classical version of the imaginarytime evolution was first proposed in [20] where the nonunitary imaginary-time evolution operator was approximated by a parameterized Ansatz state, and the parameters to obtain the ground state were found using a variational method. The QITE algorithm proposed in [19] has certain advantages, because it does not require costly optimization or ancilla qubits. When it comes to its implementation on NISQ devices, it has disadvantages over the method of [20] because of increasing circuit depth at each QITE step which raises the impact of noise from short coherence time, cross-talk between qubits, etc. Recent efforts have sought to economize the circuit depth in the QITE algorithm [21][22][23] to reduce the impact of these noise sources. In [21], we employed a method that simplified the quantum circuit needed for the unitary updates of the QITE step, thereby reducing the gate depth and noise. Here, we follow a slightly different approach for quantum circuit simplification.
For N s = 3, 4 spatial sites in the Ising spin chain with periodic boundary conditions, we used the QLanczos algorithm to compute the eigenvalues and eigenstates of the system so that transition probability, occupation number, and transverse magnetization could be calculated. We computed energy expectation values as functions of imaginary time on the IBM Q 5-qubit Yorktown device. These expectation values were obtained using QITE, and were subsequently fed to the QLanczos algorithm. We benchmarked these results against exact calculations, and obtained good agreement when error mitigation was employed.
Our discussion is organized as follows. In Section II, we introduce the model and the physical quantities to be computed. We discuss the Hilbert space and the simplifications afforded by symmetry. In Section III, we discuss the QITE and QLanczos algorithms, and the details of our hybrid classical-quantum implementation. In Section IV, we discuss the implementation of our quantum algorithm including error mitigation. In Section V, we discuss our results. Finally, in Section VI, we summarize our conclusions.

II. PRELIMINARIES
In this Section, we introduce the Ising model we used in our work and define the physical quantities we computed. We also discuss details of the Hilbert space and the simplifications one can take advantage of due to symmetry.

A. The model
The Ising model Hamiltonian with periodic boundary conditions (PBC) can be written as where X i , Y i , Z i are the Pauli matrices at the ith site, i = 0, 1, . . . , N s − 1, N s is the number of spatial sites, J is the nearest-neighbor coupling strength, and h T is the transverse magnetic field. We impose periodic boundary conditions by identifying X Ns = X 0 . At each site, we place a qubit on which the Pauli matrices act, and define the occupation number of the ith site by n i = I−Zi 2 with corresponding eigenstates |n i , where n i = 0, 1 (|0 (|1 ) denotes an unoccupied (occupied) site). A vector in the computational basis |x (x = 0, 1, . . . , 2 Ns − 1) is specified by the sites which are occupied corresponding to the digits of x equal to 1 (e.g., for N s = 4, the state |0000 has no particles, whereas |0101 consists of two particles at sites 1 and 3).

B. Unitary Time Evolution
To study the time evolution of the system, we prepare it in the initial state |initial , evolve it for time t with the evolution operator U(t) = e −iHt , and then measure it, thus projecting it onto a state |final . This process leads to the quantum computation of the transition probability In particular, in this work we study single-particle propagation and two-particle scattering. In both cases, we prepare the system in the computational basis state |initial = |x in . For single-particle propagation, x in contains a single digit equal to 1, whereas for two-particle scattering, it contains two digits equal to 1. At the end of the quantum computation, the measurement projects the system onto a different computational basis state |final = |x fin . Being in the computational basis, both initial and final states are easy to construct. However, the unitary U(t) is difficult to implement. We use the QLanczos algorithm to accomplish this, which is based on the quantum imaginary-time evolution (QITE) algorithm [19].
To calculate the transition probabilities (2), we employ a hybrid quantum-classical algorithm to solve the eigenvalue problem of the Hamiltonian (1), The unitary evolution operator is expressed in terms of the eigenvalues and eigenstates of the Hamiltonian (1) as Let t be the unitary transformation from the eigenstates of H to the computational basis. Its matrix elements are All components of the eigenstates |ψ I are real, therefore, t Ix ∈ R. This will simplify the computation of the components of the eigenstates. Scattering data can be expressed in terms of transition amplitudes between an initial and a final state, both members of the computational basis, |x in and |x fin , respectively. A transition amplitude over time t, can be calculated classically using the matrix t (eq. (5).
We obtain It should be noted that, while this calculation leads to more accurate results for NISQ devices, as we will demonstrate, for a large number of qubits, it may be more efficient to use other approaches, such as Trotterization on the evolution unitary U(t).
The time evolution of the occupation number for the ith site (i = 1, . . . , N s ) can be calculated using the expression (4) of the evolution operator. We obtain the average in the state |x at time t, where y i is the ith digit in the binary expansion of y. We deduce the transverse magnetization as One can also simulate the thermal evolution of the system [18] by computing the ensemble average of any operator O at finite temperature, T , where β = 1 k B T , k B is the Boltzmann constant, and Z = I e −βE I is the partition function. The phase transition can also be studied by using the probability of the system being in the ferromagnetic state, P FM , as an order parameter, as studied in [5] using a trapped ion quantum computer. We leave these calculations to a future study.
C. Symmetry of the system Next, we discuss the symmetry of the system and explain how it can be utilized to reduce the number of steps in quantum computations.
A conserved quantity of the system is parity, is the total occupation number. Indeed, it is easy to check that parity commutes with the Hamiltonian (1), Therefore all eigenstates of the Hamiltonian have definite parity, starting with the ground state that has even parity ((−) F = +1).
The Hamiltonian (1) is also symmetric under permutations of the sites, P : i → (i + 1)modN s , and reflection around, say, i = 0, R : i → (−i)modN s . If |ψ I is an eigenstate of the Hamiltonian (eq. (3)), then P|ψ I and R|ψ I are also eigenstates of H belonging to the same eigenvalue E I . If the energy level E I is non-degenerate, then the corresponding eigenstate must be invariant under permutation and reflection of the sites. Moreover, since R 2 = I, each energy level consists of states which are either even or odd under reflection of the spatial sites.
Let us first consider the case N s = 3. The ground state must be parity and reflection even. Since the ground state is non-degenerate, it must also be invariant under permutation of the sites. It follows that it has to be of the form There is also an excited state of this form but with different coefficients.
Other excited states are obtained by flipping all three qubits in the above expression, Next, consider an excited state which is odd under parity and reflection. These properties are incompatible with symmetry under permutation of sites, indicating that the energy level is degenerate. It is a double degeneracy with the space spanned by {|ψ 1 , P|ψ 1 } (P 2 |ψ 1 is a linear combination of the other two states, since P 3 = I, and so P 2 = −P − I). We may choose so that P|ψ 1 = 1 √ 2 (|100 − |001 ). Thus, we were able to determine the states of an excited level solely from symmetry considerations.
By the same token, there is another degenerate energy level which is obtained by flipping all three qubits, with states and easily checked to be parity and reflection even, as well as invariant under permutation. There is an excited state of the same form as the ground state and orthogonal to it.
Another excited state is of the form which is parity odd, reflection even, and invariant under permutation. There is a degenerate energy level spanned by the states which are parity and reflection odd. Another excited state is of the form which is parity odd, reflection even, invariant under permutation and orthogonal to the excited state (18). Another degenerate energy level is spanned by the states all of even parity. Another set of parity and reflection odd, degenerate higher energy level states are To access the states of the remaining energy levels, it is advantageous to flip the sign and use −H as the Hamiltonian and start by computing its ground state which corresponds to the highest energy level of H. The same symmetry considerations apply to the Hamiltonian with flipped sign, −H, and one obtains expressions for the higher-level states of H that are similar to the lower-level states obtained above.

III. ALGORITHMS
As mentioned earlier, to calculate the energy levels and corresponding eigenstates of our system we will use a hybrid quantum-classical method based on the QLanczos algorithm which uses the QITE algorithm first proposed in [19]. Therefore, in this Section we will give a brief overview of these quantum algorithms.

A. Quantum Imaginary Time Evolution (QITE)
We start by discussing the QITE algorithm whose classical counterpart was introduced in order to simulate the dynamics of many-body systems. It is advantageous to separate the Hamiltonian into local, but non-commuting, components, H = m h m . The number of these local terms in the Hamiltonian scales polynomially with the number of particles in the many-body system. Since we are only dealing with a small number of qubits, there is no need to split the Hamiltonian in our case.
QITE relies on evolution in imaginary time. To implement it, we need to set t → −iβ in eq. (4) and define the imaginary-time evolution operator U = e −βH which is no longer unitary. Starting with the state |Ψ 0 , the evolved state is found in n steps each evolving the system in imaginary time ∆τ , where n = β ∆τ , with c n being a normalization constant (c −2 n = Ψ 0 |U 2 |Ψ 0 ). In the zero-temperature limit (β → ∞), this state converges to the ground state of the system.
The QITE algorithm simulates this non-unitary imaginary-time evolution by approximate unitary updates. Thus, the sth step of the imaginary-time evolution, with s = 1, 2, . . . , n and c 0 = 1, can be approximated as where A[s] can be written in terms of Pauli operators (σ ∈ {X, Y, Z}) involving N s qubits as Once the a[s] coefficients are calculated, these unitary updates can be implemented on a quantum computer. These coefficients can be calculated up to order O(∆τ 2 ) by solving a linear system of equations (S + S T ) · a = b, where with I = {i 1 , . . . , i Ns }, and the expectation values evaluated at the state computed in the previous step, |Ψ s−1 .
These expectation values involve strings of Pauli matrices and can be evaluated with quantum algorithms recursively. By solving this linear system of equations classically, we obtain the minimum distance between |Ψ s and the unitary update (25) to lowest order in ∆τ [19]. A solution of the linear system of equations can also be found with a quantum algorithm, but we will not do this here as our focus is implementation on NISQ hardware. This kind of quantum algorithm would require implementation of unitary operations with a circuit depth that NISQ hardware could not handle. In our previous work [21], we found out that these unitary updates for the systems we considered were in the form of a unitary coupled cluster (UCC) Ansatz. This is also the case for the current Ising spin chain model.
The initial state |Ψ 0 determines which eigenstate of the system the QITE algorithm will converge to. It will converge to the ground state as long as |Ψ 0 has a finite overlap with it. For convergence to an excited state, |Ψ 0 must be orthogonal to the ground state. As we discussed in Section II C above, utilizing the symmetry of the system helps us make an educated choice of initial state. In our Ising model, we can exploit the parity and reflection symmetries to choose an initial state for QITE that will be orthogonal to low-level states and therefore converge to the desired energy level. This minimizes the number of required calculations.
The vector b has 3 Ns elements and S is a 3 Ns × 3 Ns matrix, therefore we need to perform 3 Ns (3 Ns + 1) measurements in order to calculate all elements in b and S. Since the Hamiltonian is real, so are these matrix elements. Therefore, in the calculation of b, only the elements that have an odd number of Y Pauli matrices will contribute while the rest will vanish. Similarly, the S + S T matrix elements which have an even number of Y Pauli matrices will not contribute. Additionally, the S + S T matrix is symmetric and its diagonal elements are all the same. Using this information, we can reduce the number of measurements significantly.
As explained in more detail in the next section, although all of the above steps can be performed on quantum hardware, in view of limited resources, we first computed the coefficients a[s] in the unitary updates using quantum simulation. We then implemented the unitary updates with a quantum circuit that produced |Ψ s from |Ψ 0 on quantum hardware, aided by the initialize function in the IBM Qiskit library. The partial use of quantum simulation limited the error produced by quantum hardware. If all steps are implemented on NISQ hardware, then the error we are reporting here will be larger and depend on the NISQ device used.

B. Quantum Lanczos (QLanczos) Algorithm
Next, we apply the QLanczos algorithm which uses the measurement outcomes of the QITE algorithm in order to obtain all the eigenstates of the system , including excited states. The classical Lanczos algorithm uses the Krylov space K spanned by a set of vectors {|Φ , H|Φ , H 2 |Φ , . . . }. In its quantum version (QLanc- The number of required QLanczos states |Φ l in the Krylov space is determined by the number of eigenstates of the system that have non-zero overlap with the initial state, |Ψ 0 . Our numerical calculations showed that having a smaller number of QLanczos states in the Krylov space than the number of eigenstates with non-zero overlap will result in convergence if we sample from states at a high number of QITE steps (large s).
After filling the Krylov space with QLanczos states obtained from QITE, we form the overlap (T ) and Hamiltonian (H) matrices whose elements can be calculated in terms of the energy expectation values obtained from quantum hardware (QITE), respectively, as and where r = l+l 2 , and l, l are even. The normalization constants can be calculated recursively in terms of expectation values using For experimental computation of the normalization constants we expanded Φ r |e −2∆τ H |Φ r ≈ 1 − 2∆τ Φ r |H|Φ r + O(∆τ 2 ) while keeping in mind that c 0 = 1. The expectation value Φ r |H|Φ r was calculated on a quantum computer experimentally using the states generated by the QITE algorithm. The details as to how we obtained the QITE states can be found in Section IV. For more accurate results, one can also measure the expectation values for higher powers of the Hamiltonian. The experimentally calculated normalization constants were then used to obtain the matrix elements (29) and (30). Thus, all matrix elements of T and H were computed with a quantum circuit as expectation values evaluated in the states generated by the QITE algorithm. We then solved the generalized eigenvalue equation classically and found approximations to the eigenvalues and corresponding eigenstates of the system Hamiltonian, which depended on the choice of initial state |Ψ 0 .  (1), where c −1 l |Φ l . Given the state |Ψ[E] , one can recover the approximation to the corresponding energy level using This expression for E is redundant, because we have already derived E from Eq. (32). However, due to noise the results for the energy levels deduced from (32) are numerically unstable. Thus, to obtain E, after obtaining the eigenvector x (E) from (32) classically, we engineered |Ψ[E] (Eq. (33)) by building a quantum circuit that we implemented on quantum hardware and calculated the energy expectation value (34) experimentally by performing measurements. The quantum circuit was built using Quantum Programming Studio [24] and the hardware noise from CNOT gates was reduced by Richardson extrapolation [15] in which the noise is increased purposefully by introducing double CNOT gates corresponding to the each CNOT gate in the quantum circuit and then the extrapolation of the energy expectation value was calculated to obtain the noiseless energy expectation value.
Even though this process improves the numerical stability and accuracy of the experimental results it adds to the total run time of the classical computation.
Although the noise introduced by quantum hardware increases as the system size grows, making it hard to avoid numerical instabilities, one can improve the numerical stability of the eigenvalues of the generalized eigenvalue equation (32) by applying error mitigation techniques such as Richardson extrapolation at each QITE step, or by increasing the order in the series expansion used in the calculation of the normalization constants in (31), or using a different quantum circuit simplification algorithm than Qiskit's initialize function resulting in a shorter quantum circuit with fewer CNOT gates. Work in this direction is in progress.

IV. QUANTUM PROGRAM
To calculate the time evolution of various physical quantities, we need the eigenvalues and eigenstates of the system. In our previous work, we demonstrated the practical calculation of the energy spectrum of manybody chemical and nuclear systems by implementing the QITE/QLanczos algorithm on NISQ devices [21]. Here, we extend our work to the calculation of energy levels and corresponding eigenstates of the Ising model Hamiltonian (1).
Since the QLanczos algorithm makes use of output from the QITE algorithm, we start with the calculation of energy expectation values of imaginary-time evolution with different initial states informed by symmetry considerations of the system. Using the QITE algorithm outlined above, we calculate the unitary updates (eqs. (25) and (26)) at every imaginary-time step using a small value of the imaginary-time parameter ∆τ and the Hamiltonian (1). Starting with the state |Ψ 0 , after s unitary updates, we obtain the state which we implement with a quantum circuit. We simplified these circuits following the methods discussed in [25], as implemented with the initialize function in the IBM Q Qiskit library. Examples of 3-and 4-qubit quan-tum circuits for the states (35) are depicted in Fig. 1 in terms of single-qubit rotation gates R y (θ) and two-qubit CNOT gates. At every imaginary-time step, the angles change, as they depend on the state |Ψ s , but the depth of the circuit remains the same. Therefore, in terms of economizing the number of gates and operations in the quantum circuit, our results are similar to those in our earlier work [21]. It should be noted that, depending on the topology of the quantum hardware, interactions between physical qubits matching those in the quantum circuit implementing (35) may not be readily available, necessitating the addition of SWAP gates to the circuits in Fig. 1.

A. Error Mitigation
Running the quantum circuits on NISQ devices brings errors of various sources such as noise from the implementation of the circuit gates and noise due to the measurement readout errors. To mitigate these errors in the measurements error mitigation strategies are employed. In this work, we only use a readout error mitigation technique in calculation of the energy expectation values at each QITE step. One can use further error mitigation strategies such as Richardson extrapolation as we did in ref. [21] or reduced density matrix purification ( [26]) to improve the results obtained using the QITE algorithm.
In this paper, we use local readout error mitigation strategy that we used in our previous work [21] in which the corrected expectation values of the Pauli terms is calculated using where p(x) is the probability of each qubit outcome and it takes 2 N values. Here, we only consider the expectation values for Z terms since we do the measurements in Z basis. The terms with X and Y Pauli operators are rotated to be measured in Z basis. We define the symmetric and anti-symmetric combinations of the probability of i-th qubit flipping from 0 to 1 (p i (0|1)) or from 1 to 0 (p i (1|0)) as with p(1|0) = # of states expected in |1 measured in |0 # of shots (38) or vice versa for p(0|1).

A. QITE and QLanczos Results
The experiments for the QITE algorithm were run on 5-qubit IBM Q Yorktown hardware. The number of shots for the each experiment was 8192 and each experiment was run N runs = 3 times to calculate the statistical error in the measurements. The reason for choosing this quantum computer out of other IBM Q's cloud accessible devices is its periodic topology as seen in Fig. 2. Using a quantum computer with periodic topology reduces the number of required SWAP gates for our periodic Ising spin chain Hamiltonian which reduces the number of required CNOT gates. This is important because CNOT gates are the dominant source of the error in a quantum circuit. For comparison, Honeywell's ion trap quantum computer offers connectivity between all physical qubits. Therefore, the error in this type of quantum system might be smaller since it does not require the addition of SWAP gates for the type of interaction Hamiltonian considered here. The basis gates which can be directly implemented on IBM Q Yorktown quantum computer are single-qubit gates U and the two-qubit CNOT gate, where is a general three-parameter single-qubit gate. In Fig.  1, we used the single-qubit rotation gate R y (θ) which can be expressed in terms of the basis gates as R y (θ) = U (θ, 0, 0). As mentioned in the Section III, simulating each QITE step requires significant number of measurements on hardware. Even using the aforementioned properties of b and S matrices there needs to be a large number of measurement done to apply the QITE algorithm on hardware. For example, for N s = 3 we were able to reduce the number of measurements from 756 to 187 at every QITE step. Due to limitations in cloud access to the quantum hardware (such as long queue and connection interruptions) we simulated the quantum circuits for the states |Ψ s and implemented them on quantum hardware to obtain the energy expectation values for various values of imaginary time. With full implementation on a NISQ device, additional errors will occur. To estimate these additional errors, we considered a generic case and fully implemented it on simulated quantum hardware. We obtained energy expectation values for various values of imaginary time for three sites, N s = 3, using the initial state |Ψ 0 = |100 , and the Ising model with parameters J = 0.6 and h T = 1. We implemented the QITE algorithm and obtained the operator A[s] from measurements on the noisy simulator of the same backend. We used N shots = 8192 and the calibration parameters from 04/24/2020. In Fig. 3 we compare the convergence of the energy expectation values to the first excited state energy in three different cases, (a) from exact calculation of the state |Ψ s as well as energy expectation values, (b) from a noisy simulation of the state |Ψ s and exact energy expectation values, and (c) from a noisy simulation of both the state |Ψ s and readout error mitigated (notated as ROEM in the rest of the paper) energy expectation values. The energy expectation values obtained using methods (a) and (b) are very close to each other, showing that the main source of additional error is due to measurements. It follows that the use of simulated states does not introduce significant errors. However, the ROEM energy expectation values obtained from measurements on quantum hardware differ from results from noiseless simulations. In what follows, we use noiseless simulated states, implement their quantum circuits on quantum hardware, and perform measurements to obtain energy expectation values. The results of measurements on these quantum circuits produced by the QITE algorithm as a function of imaginary time for different initial states are depicted in Figs. 5, 6, and 7 for N s = 3 and N s = 4 spatial sites of our Ising model with parameters h T = 1 and J = 0.6. In the N s = 3 case, the initial states |Ψ 0 are chosen as |100 , |010 , 1 √ 3 (|110 + |101 + |011 ), and |111 shown in Fig.   5 (a)-(d), respectively.
The QITE algorithm converges to the minimum of the symmetry group that the initial state belongs to. Therefore, it might be challenging to access higher-value energy levels using the QITE and QLanczos algorithms. To facilitate the algorithm's convergence to higher levels, we reversed the sign of the Hamiltonian (1) so that high energy levels turn into low levels whereas the corresponding eigenstates remain the same. We applied this strategy to calculate some of the high energy levels and corresponding eigenstates of our system, e.g., for the 4th and 5th excited states in the 3-qubit (N s = 3), and the 15th excited state in the 4-qubit (N s = 4) case. The results of this strategy for the QITE algorithm, including energy expectation values, can be seen in Figs. 7(e) and 6. In these examples, since we are looking for the minimum of the reverse Hamiltonian −H, we chose the initial states |Ψ 0 to be reflection and parity symmetric, namely |110 , |011 , |101 , and |0000 , respectively.
As mentioned in Section II C, some of the eigenstates are completely constrained by the symmetry of the system, therefore calculating them is redundant. For example, in the N s = 4 case for parameters J = 0.6 and h T = 1 the zero eigenvalue is degenerate and the corresponding exact eigenstates are given by eq. (21). Similarly, the eigenstates correponding to the degenerate energy level −2 are given analytically by eq. (19). We took advantage of the exact expressions for these eigenstates in our calculations. Although, we used the exact eigenstates obtained using symmetry constraints, we measured the energy expectation values for the eigenstates demonstrated as the first state in (21) and states in (19) on hardware (IBM Q Yorktown) using the quantum circuits seen in Fig. 4 (a), (b), and (c), respectively. These circuits were run on hardware N runs = 3 times with each run having N shots = 8192 on 08/12/2020 with qubit layout [q 0 , q 1 , q 2 , q 3 ] = [1, 0, 3, 2] and the ROEM average energy values obtained are 0.037 ± 0.006, −2.06 ± 0.02, and −2.01±0.01 (where the ± error is the standard deviation of the mean) compared to the exact eigenvalues of 0 and -2, respectively. For the same coupling and magnetization parameters the states expressed in (22) correspond to eigenvalue 2 and it is degenerate. Although these states correspond to 3 occupied sites and since we study single particle propagation and two-particle scattering only we did not need them in our calculations we obtained an experimental ROEM mean value of 2.05 ± 0.03 for the first state in (22) using the circuit in Fig. 4(d). The experiments were run on IBM Q Yorktown, N runs = 3 times with each run having N shots = 8192 on 08/13/2020 with qubit layout [q 0 , q 1 , q 2 , q 3 ] = [0, 1, 2, 3].
The circuit for | 0 i = |1110i |1011i is as seen below.
The circuit for | 0 i = |0010i |1000i is as seen below.
The circuit for | 0 i = |0101i |1010i is as seen below.
|0i  In our current study, we used two-dimensional Krylov spaces. Although, depending on the choice of initial state, convergence might take longer for a low-dimensional Krylov space, adding more dimensions causes numerical instabilities and does not guarantee convergence to eigenstates of the Hamiltonian (1). In- terestingly, we were able to observe convergence to the eigenvalues of the system by using a three-dimensional Krylov space together with our uncertainty criterion to exclude spurious states (∆E ≤ δ). However, we did not obtain three distinct energy eigenstates. In general, results were numerically more accurate in two-dimensional Krylov spaces for the QLanczos algorithm. Adding more dimensions decreased the number of cases where offdiagonal T matrix elements were < 1 resulting in spuri-ous eigenstates. Our numerical results indicate that using two-dimensional Krylov space is the optimal choice for the implementation of the QLanczos algorithm on noisy quantum devices of this particular system with the parameters used in this study. Further application of the error mitigation strategies, such as Richardson extrapolation (an example of application of Richardson extrapolation to QITE algorithm can be seen in [21]), might improve the numerical stability and can provide faster convergence in higher-dimensional Krylov spaces. As mentioned in Section III, we decide on the convergence to the eigenstates and eigenvalues of the system and discard spurious states by using the uncertainty criterion, ∆E ≤ δ. Two examples involving the ground and excited states that are specific to a given initial state, |Ψ 0 , are shown in Fig. 9. Specifically, the uncertainty ∆E is shown for N s = 4 and various values of (l, m), where l, m are even integers and label the basis states of the Krylov space, which is spanned by {|Φ l , |Φ m }. Results of the 3 different runs on IBM Q Yorktown hardware are shown. We keep increasing l and m, which correspond to QLanczos states with higher QITE steps, until ∆E < 1, and we choose the eigenvalues and eigenstates that give the minimum uncertainty.
After the application of this process, as explained earlier, we ran each quantum circuit corresponding to the each eigenvector obtained from our hybrid quantumclassical QLanczos algorithm on quantum hardware, specifically on IBM Q Vigo, Casablanca, Manhattan devices depending on their availabilities. This gave us the experimental energy eigenvalues for N s = 3 and N s = 4 with parameters J = 0.6 and h T = 1 PBC Ising model. As a result, using either the exact states obtained from symmetry or using our QLanczos algorithm with a Krylov space of size 2 we obtained experimental energy eigenvalues as seen in Fig. 8

B. Time Evolution Results
Finally, we obtained the coefficients t Ix in (5) from the measurements of each component of the eigenvector. These measurements give the absolute value of each coefficient, |t Ix |. Since they are all real, in order to determine them, we need to find the sign. This requires additional measurements with an ancilla qubit, but they introduce no errors because of the binary nature of the sign. We used these coefficients in Eqs. (7), (8), and (9) to calculate the transition amplitudes, occupation number, and average magnetization as functions of time.
We summarized our method to calculate the transition amplitudes, occupation number, and transverse magnetization using QLanczos algorithm in the pseudocode below in Fig. 10. If the algorithm fails to find eigenvalues and corresponding eigenvectors of the Hamiltonian, then the initial parameters should be changed. If the choice of initial state is informed by symmetry considerations and the uncertainty keeps decreasing at each step, one can start the algorithm with a larger s max . In Fig. 10 we denoted the measurements needed to be performed on quantum hardware as [Q] and the classical computations were denoted by [C]. As discussed earlier, due to constraints in quantum hardware access, we calculated A[s] in Step 3 classically and used the Qiskit initialize function to find the quantum circuit in Step 4 of the pseudocode in order to reduce the depth of the circuit.
Here, we present our experimental data obtained from data on the IBM Q Yorktown, Vigo, Casablanca and Manhattan hardware for transition probability amplitudes , occupation number at each spatial site, and average transverse magnetization for number of spatial sites N s = 3 and N s = 4.
We chose the parameters of the system Hamiltonian in (1) to be h T = 1 and J = 0.6. In ref. [3] it was found that errors arising from quantum hardware become worse as the coupling J increases. In this section we present results which show small hardware errors even as one moves away from the weak coupling regime.
In Figs. 11 (a)-(d), we show numerical values of the transition amplitudes calculated from given exact |initial and |final states, and compare them with values obtained from experimental data produced by the QLanczos quantum algorithm that calculates energy eigenvalues and corresponding eigenstates. Figs. 11 (a) and (c) show the one-particle propagation probability, and Figs. 11 (b) and (d) show the probability of two-particle scattering.
Similarly, in Fig. 12 we show a comparison between the numerical value of occupation numbers at various spatial sites calculated from a given exact |initial state and the one calculated experimentally from the energy eigenvalues and corresponding eigenstates using the QLanczos algorithm.
It should be noted that when the particles are initially at sites 0, 2 (i.e., |initial = |1010 ) the time evolution of the occupation number at even (odd) sites is the same, i.e., n 0 (t) = n 2 (t) ( n 1 (t) = n 3 (t) ).
Finally, in Fig. 13 (a)-(d), we present a comparison between the numerical value of the average transverse magnetization calculated from a given exact |initial state and the experimental average transverse magnetization obtained from energy eigenvalues and corresponding eigenstates using the QLanczos algorithm.
Figs. 11, 12, and 13 demonstrate that for number of sites N s = 3, the exact results and experimental data are in very good agreement. For a larger system (N s = 4), the exact and experimental data are still in good agreement. In the latter case, the quantum circuit used to calculate the energy expectation values includes more single-qubit rotation and CNOT gates, which result in more error in the measurements. This can be seen by comparing Fig. 7 with Fig. 5 in Section III.

VI. CONCLUSION
In this work, we discussed a hybrid quantum-classical method to calculate physical properties of the Ising spin chain model as a function of time, such as transition amplitudes, occupation numbers at various sites, and transverse magnetization, using the QLanczos algorithm as a tool. We took advantage of the symmetry of the system to simplify the quantum computation of the eigenvalues and eigenstates of the Hamiltonian of the system which were then used for the computation of various physical quantities of interest. We ran experiments for the  b)). The parameters are set to J = 0.6 and hT = 1. For Ns = 4, the 2nd, 3rd, and 6th energy levels were obtained from the exact eigenstates (Eqs. (19) and (21)). The remaining eigenvalues were obtained from the eigenstate quantum circuits ran on IBM Q Yorktown, Vigo, Casablanca and Manhattan devices. The Richardson extrapolation error mitigation strategy was applied to reduce effect of the quantum hardware noise. The experiments were run Nruns = 3 times and the error bars represent one standard deviation.  Fig. 7(a)) and (b) |0000 (cf. with Fig. 7(f )). Experimental results from Nruns = 3 runs on IBM Q Yorktown hardware.
QITE algorithm on IBM Q Yorktown hardware eigenvector quantum circuits on IBM Q Vigo, Casablanca and Manhattan devices for N s = 3 and N s = 4 spatial sites.
Our results show good agreement with the exact values of the physical quantities of interest. It should be pointed out that although the use of the initialize function in the IBM Qiskit library gives energy expectation value calculations at each QITE step which are very close to the exact value in the noiseless simulator case, our results show how different the noisy simulator and the hardware data can be from each other as well as exact calculations. Our data constitute the first demonstration of quantum imaginary-time evolution in a 4-qubit system on NISQ hardware, and can be useful for benchmarking purposes.
Notably, the use of the symmetry of the system in simplifying the QITE and QLanczos algorithms reduces the number of steps in the quantum calculations, which leads to a significant reduction in error due to NISQ hardware. The QITE and QLanczos algorithms converge to the minima determined by the symmetry subgroup of the chosen Kubra December 9, 2020  initial state. Further, higher excited states can be obtained by reversing the sign of the Hamiltonian as needed. These two features enabled us to find energy levels that otherwise were difficult to compute due to the numerical difficulty associated with increasing the number of vectors in the Krylov space.