Quantum simulation of scattering in the quantum Ising model

We discuss real time evolution for the quantum Ising model in one spatial dimension with $N_s$ sites. In the limit where the nearest neighbor interactions $J$ in the spatial directions are small, there is a simple physical picture where qubit states can be interpreted as approximate particle occupations. Using exact diagonalization, for initial states with one or two particles, we show that for small $J$, discrete Bessel functions provide very accurate expressions for the evolution of the occupancies corresponding to initial states with one and two particles. Boundary conditions play an important role when the evolution time is long enough. We discuss a Trotter procedure to implement the evolution on existing quantum computers and discuss the error associated with the Trotter step size. We discuss the effects of gate and measurement errors on the evolution of one and two particle states using 4 and 8 qubits circuits approximately corresponding to existing or near term quantum computers.


I. INTRODUCTION
There has been a fast growing interest for quantum computation in the context of high energy and nuclear physics [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. One important motivation is to calculate real-time evolution of states in large Hilbert spaces which cannot be handled with standard sampling methods. The long term goals include jet physics and early cosmology. However, in the near term it is important to demonstrate that it is possible to make progress towards these major goals using quantum computers or quantum simulation experiments with a limited number of qubits [18].
Numerical lattice gauge theory started in the late 70's by studying Z 2 (Ising) gauge theories on 3 4 lattices and has steadily developed as a reliable tool that today allows different collaborations to compare numerical estimates for hadronic processes with errors of a few percent. It thus seems natural to start the study of real time evolution using the quantum Ising model in 1+1 dimension [5] or the Schwinger model [2,4].
In the following we propose to consider simple cases of time evolution for the quantum Ising model with a number of sites of the same order as the number of qubits in devices existing or expected to exist in the near future. The model has a second order phase transition which allows the use of finite size scaling (FSS) to extract interesting information using systems with a small number of lattice sites. This strategy is explained in Ref. [19]. The main goal is to provide reliable benchmark calculations in situations that will allow comparison among different platforms.
The paper is organized as follows. In Sec. II, we present the model and in Sec. III we present some perturbative results. In Sec. IV the Suzuki-Trotter formulation of the time evolution operator for the quantum Ising model is derived. In Sec. V we discuss how artificial noise is introduced into our simulations as well as the results of our simulations for both free propagation and "scattering" of particles. We follow the methodology inspired by the one laid out in Ref. [1]. We first prepare highly localized "wavepackets" which can be considered a two-particle state. Then we let the packets spread and "interact" using the corresponding Hamiltonian.

A. Hamiltonian and boundary conditions
The one-dimensional quantum Ising model is the standard example of a quantum field theory with continuous time that is obtained from a classical lattice model with one extra dimension corresponding to the Euclidean time [20,21]. In this example, the classical model is the usual two-dimensional Ising model solved by Onsager [22] and Kaufman [23]. The Hilbert space of the quantum model is a tensor product of qubits and the connection to quantum computing is immediate (see Eq. (2) for an illustration) The connection to the classical model makes the choice of a basis where the nearest neighbor interactions in the spatial direction are diagonal very natural. We call this choice the "spin basis". For reasons that will become clear soon, we use a representation where the other term, often referred to as the transverse magnetic field term, is diagonal. We call this choice the "particle basis." The two representations are connected by a Hadamard unitary transformation.
In the particle basis, the nearest neighbor interactions (particle hopping) use adjacent pairs ofσ x , governed by J, while the transverse magnetic field interactions (onsite coupling) governed by h T useσ z , whereσ z andσ x are the traditional Pauli matrices. In the particle basis, we define the "particle number" at each l site aŝ We will use these quantum numbers to specify the Hilbert space which is a direct product of two-dimensional qubit spaces at each of the N s spatial sites. Just to give an arXiv:1901.05944v2 [hep-lat] 4 Feb 2019 example for N s = 4, the action of a sample operator on a sample state can be illustrated aŝ We will see that the Hamiltonian only connects states for which the total particle number (n = in i ) is the same modulo 2.
We define the Hamiltonian corresponding to open boundary conditions, hereinafter referred to as (OBC), as, The Hamiltonian corresponding to periodic boundary conditions (PBC) is defined as, while it is also interesting to consider a Hamiltonian with antiperiodic boundary conditions (ABC), For any of these Hamiltonians we define an operator to carry out the exact time evolution in units where =1 as, B. Symmetries The model has a Z 2 global symmetry corresponding to flipping all the spins in the spin basis. In the particle basis, this corresponds to multiplying the states byσ z at each site. This defines a unitary transformation that flips the sign of the operatorσ x at each site. As such operators come in nearest neighbor pairs, the transformation leaves the Hamiltonian invariant regardless of boundary conditions. This is equivalent to saying that the particle numbern defined above is conserved modulo 2.
When N s is even, it is also possible to invent a twostep transformation which changes the sign of the entire Hamiltonian. Similar equivalences appear for classical gauge theories [24]. We first apply aσ x transformation at each site. This changes the sign of the on-site term and leaves the hopping term unchanged. In a second step, we apply aσ z on every other site. The full transformation flips the sign of both terms of the Hamiltonian. Consequently, all the states appear in pairs with opposite signs. This property appears clearly in Fig. 1 for N s = 4. In this case, the 16 states split into approximately degenerate groups of 1 (0 particle), 4 (1 particle), 6 (2 particles), 4 (3 particles or 1 hole), and 1 (4 particles) when J << h T . In the next section, we discuss the splitting using degenerate perturbation theory. Perturbation theory allows for very accurate calculations of real-time evolution when J is sufficiently small. Here h T is the transverse magnetic field and J is the strength of the nearest-neighbor interaction. As J/h T increases the degenerate energy levels split.

III. APPROXIMATE EVOLUTION FOR J << hT
A. Approximate particle description In the limit where J = 0, we obtain a very simple picture for the quantum Ising model. The energy is then the sum of the on-site energies. We have a unique ground state where all sites have an energy −h T and so E (0) = −N s h T . We now have degenerate "one-particle" states where one on-site state with energy +h T can be placed at N s locations. If the h T energy is located at the site j, we call this state |j . These states have an energy −(N s − 2)h T . Similarly we have N s !/(n!(N s − n)!) totally antisymmetrized states with n "particles" and an energy (−N s + 2n)h T . The effect of the nearest neighbor interactions can be included perturbatively [21]. The model can also be solved exactly by performing a Wigner-Jordan transformation [23]. However, at finite volume, boundary conditions should be treated carefully. To be more explicit, the term a † N s a 1 needs to be supplemented with a product ofσ z l in order to reproduce the original spin Hamiltonian which requires a separate discussion for the even and odd sectors [23].

B. One particle
At order J in the one-particle sector we have a particle hopping that stays in the one-particle sector. It is worth noting that the particle conservation picture makes this model in the small J limit equivalent to the XY model which has been studied thoroughly. In particular the eigenstates and energies of this model are discussed in Refs. [25,26] and the zero field case is examined in Refs. [27][28][29], time dependent zz spin correlations were studied in Refs. [30,31] and xx spin correlations were studied in Ref. [32]. If periodic boundary conditions are imposed, as in Eq. (4), Fourier modes diagonalize the perturbation. This lifts the degeneracy by a term proportional to 2Jcos(2πm/N s ). The perturbation also contains operators that connect to the 3-particle states; this leads to energy shifts O(J 2 /h T ). If we neglect these second order effects, we have a simple approximate quantum mechanical behavior.
We can then prepare the system in an initial state |ψ and calculate ψ(t)|n l |ψ(t) , where the calculations in the quantum mechanical approximation are relatively easy. For instance, for |ψ(0) = |j , we obtain where the "discrete" Bessel functions are defined as, which corresponds to the usual definition in the limit of large N s . In fact these "Bessel" functions appear in the XY model for the zz correlations as shown in Ref. [31]. The approximation is accurate when t is less than O(h T /(J 2 )) (see Fig. 2). The implication of this is that pair creation in this model is driven by the hopping parameter J, rather than the size of the model. The peaks in the occupation values shown in Fig. 2 suggest that for J = 0.02 perturbation theory would be accurate for a system up to 32 sites because approximately 4 resurgences happen before noticeable discrepancies begin to appear.

C. Two particles
The results for one-particle states can be generalized to two-particle states provided that ABC are used. Using lowest order degenerate perturbation theory with ABC we find that the occupation number is, where |ij = |0...0, 1 i , 0...0, 1 j , 0... . Agreement is excellent for long time scales, with only a small discrepancy for (c) in Fig. 2. It is interesting that after a time, such that Jt is of order 1, the three types of boundary conditions start to give very different values of n i (t) (see Fig. 3). For open boundary conditions, one can compare the situation with that of an ideal gas where the forces exerted on the particles are due to the walls and generate the pressure.

D. Finite volume corrections
It is possible to calculate the difference between the finite volume discrete Bessel functions J (Ns) n (x) and the usual infinite volume expressions J n (x). Using the Poisson formula, one finds that where the sum over runs over strictly positive and negative integers.The difference between J (Ns) n (x) and J n (x) is small for small argument. This is illustrated for N s = 8 and n = 0 in Fig. 4. One sees that the difference between J  Since a quantum computer cannot exactly implement the operator given in Eq. (6), we need to use the Suzuki-Trotter (ST) approximation to evaluate the time evolution. We use the first order approximation in order to limit the gate depth of the system: However we have to apply this operator multiple times to evolve the system to some final time t. This iterative process leads to a new expression for the time evolution operator: (12) While the estimated worse case error O(t δt) is true in general, this bound over estimates the ST truncation error, which should be approximately O(Jh T tδt) because theσ z terms only add phase-shifts to the state-vectors of the Hilbert space, and does not affect any measurement of the the basis state. The next order correction to the ST formula is, which can be found using the methodology proposed in [33]. The second order ST approximation essentially becomes the first order approximation; at this point we can justify the error beginning at the second order ST approximation being O(J 3 t(δt) 2 ).
It is relatively straightforward to implement the simplest ST approximation as a quantum circuit (see Fig. 5). The Hamiltonian is split as follows: The H 1 term can easily be implemented as a single moment in a quantum circuit. While H 2 and H 3 commute, they have to be implemented as separate moments in the quantum circuit because each of these contain terms which entangle two qubits. This corresponds to the idea [34] of splitting the Hamiltonian into pieces that can be implemented easily separately and when combined correspond to the original Hamiltonian with a Trotter error. On the other hand, the on-site terms can be executed in a single moment in a quantum circuit as these are single qubit operators. The implementation of this circuit for an arbitrary number of qubits has a gate depth (d l ), and a total number of gate operations: where N s is the number of sites and N t is the number of trotter steps. We examined two different cases of the 1-D Ising model with 8 sites (with both OBC and PBC): the time evolution of a single particle initial state and scattering of two particles. For all cases we examined the system with the nearest neighbor coupling J = 0.02 and on-site coupling h T = 1.0. We define the initial states for the system in Table I  Because the ST approximation is an iterative process we want to know how many operations can be carried out before imperfections in the approximations and the noisy gates of a quantum computer will produce substantial issues with our simulations. A first step is to measure the fidelity of the Trotter operator with the exact evolution operator over the time scales of interest in these processes (see Fig. 6A and Fig. 6B). The fidelity of the ST operator is, For both OBC and PBC it appears that a trotter step of δt = 10.0 for J = 0.02 is satisfactory to describe the time evolution for small time scales.
It is important to have an understanding of how frequent quantum gate errors will be when applying the ST operator. We expect the number of gate errors to increase the more Trotter steps we apply. In Fig. 7, we show the average number of gate errors (Pauli Channel) at a given evolution time for a state of the art trapped ion system (p 1 qubit = 1.0 * 10 −4 and p 2 qubit = 5.0 * 10 −4 [35,36]) and for slight improvements of current typical digital quantum computers [37,38] The Pauli error channel for single-qubit gates is defined in terms of the density matrixρ: E(ρ; p x , p y , p z ) = (1−p)ρ+p xσ xρσx +p yσ yρσy +p zσ zρσz . (19) The values p x , p y , and p z correspond to the probabilities of an σ x , σ y , or σ z error respectively occurring and p = p x + p y + p z . The error channel for two qubit gates is given by E (2) = E E. The values of p x , p y , p z for the one and two-qubit gates are given in Tables II and  III. We also introduce measurement errors into our simulations. These are caused by misidentifying the state that the qubit is in (i.e. reading a 1 as a 0 or a 0 as a 1). We implemented this by changing the readout value with a chance p measure given in Tables II and III. Ref. [39] identifies a procedure to address the readout error in the supplementary material; we simplify their result to using the following rescaling of the measured readout:   Approximating the time evolution operator using the ST method introduces an error O(δt t). It was suggested in [40][41][42] that by using simulations at multiple different Trotter steps and noise levels it is possible to systematically reduce the uncertainty in the measured quantities.
In order to minimize the noise error, Ref. [41] suggested several methods, one of which is an exponential extrapolation, In Eq. (21), is the noise rate for the system which is dependent upon the probability of gate errors occurring as in Eq. (19) and r is a scale factor such that r > 1. Due to the computational overhead of carrying out this methodology of error mitigation we only demonstrate a modification of it in our results for free propagation on a four site lattice. Our modification to the error mitigation scheme proposed in [40][41][42] involves changing the extrapolation equation proposed in Eq. (21) to an exponential ansatz of the form: This form should retain the same general behaviors of the ansatz proposed in Eq. (21).
The methods for reducing algorithmic errors are more computationally intensive and discussed in Ref. [40]. They argue that the algorithmic error rate, N , scales as 1/N , where N is the number of trotter steps. Due to the increased computational demand of this error mitigation method we currently have not implemented it.

C. Simulation Results
The initial state of the system for each case we looked at is given in Table I. We evolved the system using exact diagonalization of the Hamiltonian as well as simulated using a quantum virtual machine implemented through the python library QISKIT published by IBM, which implemented the noise corresponding to the Pauli channel for the system. The results from the QISKIT simulations were checked for consistency by implementing the noise channel using matrix operations instead of the QISKIT library. The results of the simulation are show in Fig.  8A through Fig. 9B. We use a reduced χ 2 value (χ 2 ) as a metric for how well the simulated quantum computer predicts the exact results, we do this by comparing the difference at different sites, These values over time are given in Table IV and Table  V. The systematic errors from the ST approximation become significant at time scales on the order of Jt ≈ 6. These results suggest that it is possible to simulate a real time scattering event and measure the matrix elements of a correlation function, for a simple field theory.   In addition it would be interesting to see if it is possible to extract results using current superconducting qubit quantum computers. In this case we used p x = p y = p z = 0.0005 for the one qubit gates and p x = p y = p z = 0.004 for the two qubit gates; this corresponds to a gate error of ∼ 0.01 for one qubit gates and ∼ 0.04 for two qubit gates. In this case we need to address the issues of noisy quantum gates because of the number of two qubit gates we have in our system. To do this we simulated the system at 4 different noise levels by introducing noisy identity QISKIT simulation for current trapped ions, red diamonds: numpy simulation for current trapped ion, blue triangle: QISKIT simulation for near future superconducting qubit quantum computers, cyan right arrow: numpy simulation for near future superconducting qubit quantum computers, gray bars: exact diagonalization operators into our circuit. Similarities can be seen between Fig. 11 and Fig. 4 in Ref. [4]. The only difference in the procedure that we carried out is the extrapolation method that we used. Ref. [4] used a quadratic ansatz while we used an exponential ansatz of the form: O( * r; t) = AB r + C (24) to extrapolate the noiseless observable. We used priors of A = 0.0 ± 0.5, B = 0.0 ± 1.0, and C = 0.5 ± 0.5. For proof of concept we worked at 4 sites with J = 0.02 and h T = 1.0 and we took 8000 measurements for each data point at each noise level. The results of these simulations are shown in Figs. 10 and 11. The errors found at later times in Fig. 12 are likely due to the signal and the noise level being so close together that it is difficult for a fit to yield an accurate noiseless extrapolation.

E. Continuation to larger J
We have examined tentatively regions where J = 0.2 and h T = 1.0. The most noticeable effect is the particles hop between sites far more quickly. While this case is still far away from the continuum limit, i.e. J = h T , in this regime pair creation and annihilation becomes more frequent. This implies that the particle conservation picture breaks down quicker and second order effects in degenerate perturbation theory become significant. In addition we find that the particle occupation at given sites is not as regular as in the case where J << h T .
The most significant change in working with a larger value of J is the Trotter time steps must be shrunk because the Trotter truncation error still scales in a similar manner as the regime J = 0.02. We expect some issues arising from noisy simulations will be similar to those encountered in the Jδt = 0.1 shown in Fig. 12, with noiseless extrapolation. Specifically when the lowest noise simulation observable is close to, or crosses, the observed value for inflated noise simulations which no longer have a discernible signal, the noiseless extrapolation method produces substantially larger uncertainties for the observable.

VI. CONCLUSION
We have demonstrated through simulations on an emulated quantum computer that it is possible to use current trapped ion systems to simulate the real-time evolution of the quantum Ising model with both 4 and 8 sites, and in the near future it will be possible to be simulated on quantum computers using superconducting qubits. Currently, the density matrix renormalization group and tensor networks are the only methods we have of examining real-time scattering; however, in the near future quantum computers will be able to join this group of tools so that we can examine these systems in real time. We have derived a simple perturbative expression that can be used to check the consistency of the results done on a system of trapped ions or in the near future on superconducting qubits. These perturbative expressions can be used for much larger systems and can be easily handled analytically and numerically.
Much work remains to be done in order to study the real-time evolution of interacting particles close to the continuum limit. We plan to examine related theories, such as the O(3) non-linear sigma model where the triplet and singlet states could be implemented with a pair of qubits, or slight modifications to the Ising model such as changes in the transverse field, because these models allow us to examine a richer volume of observables such as phase shifts, scattering cross sections, and bound states.