Quantum simulation of scattering in the quantum Ising model

TL;DR

This paper explores quantum simulation of the real-time dynamics in the 1D quantum Ising model, demonstrating accurate particle evolution predictions with small interactions and analyzing implementation errors on near-term quantum computers.

Contribution

It introduces a method to simulate quantum Ising model dynamics using Trotterization on quantum computers, with analytical solutions for small interactions and error analysis for near-term devices.

Findings

Discrete Bessel functions accurately model particle evolution for small J.

Boundary conditions significantly affect long-time evolution.

Error analysis for Trotter steps and quantum hardware noise.

Abstract

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.

Tables

Table 1. Table 1 : Initial State of the system

Sim Type	OBC	PBC
Free Prop	$|10000000⟩$	$|10000000⟩$
Scatt	$|10000001⟩$	$|10001000⟩$

Table 2. Table 2 : Optimistic error rates for current trapped ions

pauli ch.	1 qubit error	2 qubit error
$p_x,p_y$	0.00002	0.0001
$p_z$	0.00006	0.0003
meas. ch.	1 qubit error	2 qubit error
$p_{\text{measure}}$	0.05	-

Table 3. Table 3 : Optimistic error rates for near term superconducting qubits

pauli ch.	1 qubit error	2 qubit error
$p_x,p_y$	0.00033	0.0033
$p_z$	0.00033	0.0033
meas. ch.	1 qubit error	2 qubit error
$p_{\text{measure}}$	0.05	-

Table 4. Table 4 : χ ~ 2 superscript ~ 𝜒 2 \tilde{\chi}^{2} at different times over all sites for free propagation

metric	obc	pbc
${\tilde{\chi }}^2\left(0\right)$	1.5	0.91
${\tilde{\chi }}^2\left(80\right)$	4.41	0.5
${\tilde{\chi }}^2\left(160\right)$	2.43	0.45
${\tilde{\chi }}^2\left(240\right)$	4.4	0.2
${\tilde{\chi }}^2\left(320\right)$	0.88	0.13
${\tilde{\chi }}^2\left(400\right)$	3.2	0.28
${\tilde{\chi }}^2\left(480\right)$	4.1	0.33

Table 5. Table 5 : χ ~ 2 ​ ( t ) superscript ~ 𝜒 2 𝑡 \tilde{\chi}^{2}(t) for scattering and mixed field simulations for 8 sites; the number of degrees of freedom is 8

metric	obc scattering	pbc scattering
${\tilde{\chi }}^2\left(0.0\right)$	2.1	1.7
${\tilde{\chi }}^2\left(40.0\right)$	7.4	0.71
${\tilde{\chi }}^2\left(80.0\right)$	4.3	0.525
${\tilde{\chi }}^2\left(120.0\right)$	1.4	0.19
${\tilde{\chi }}^2\left(160.0\right)$	2.9	1.1
${\tilde{\chi }}^2\left(200.0\right)$	9.1	0.38
${\tilde{\chi }}^2\left(240.0\right)$	12.9	0.6

Equations

$${\hat{n}}_l=(1-{\hat{\sigma }}_l^z)/2.$$

$${\hat{\sigma }}_3^x|1011⟩=|1001⟩.$$

$$H_{obc}=-J\sum _{i=1}^{N_s-1}{\hat{\sigma }}_i^x{\hat{\sigma }}_{i+1}^x-h_T\sum _{i=1}^{N_s}{\hat{\sigma }}_i^z.$$

$$H_{pbc}=H_{obc}-J{\hat{\sigma }}_1^x{\hat{\sigma }}_{N_s}^x,$$

$$H_{abc}=H_{obc}+J{\hat{\sigma }}_1^x{\hat{\sigma }}_{N_s}^x.$$

$$U(t)=e^{-itH}.$$

$$⟨{\psi }_j(t)|{\hat{n}}_l|{\psi }_j(t)⟩\simeq |J_{l-j}^{(N_s)}(2Jt){|}^2,$$

$$J_n^{(N_s)}(x)=\frac{(-i{)}^n}{N_s}\sum _{m=0}^{N_s-1}{e}^{i((\frac{2\pi mn}{N_s}+x\text{cos}(\frac{2\pi m}{N_s})))}$$

$$⟨i,j(t)|{\hat{n}}_l|i,j(t)⟩\simeq |J_{l-i}^{(N_s)}(2Jt){|}^2+|J_{l-j}^{(N_s)}(2Jt){|}^2,$$

$$J_n^{(N_s)}(x)=J_n(x)+\sum _{\ell \ne 0}(i{)}^{N_s\ell }{J}_{n+N_s\ell }(x),$$

$$U(\delta t)\simeq e^{i\delta tJ\sum {\hat{\sigma }}_i^x{\hat{\sigma }}_{i+1}^x}{e}^{i\delta t{h}_T\sum {\hat{\sigma }}_i^z}+\mathcal{O}(\delta t^2).$$

$$U(t;\delta t)\simeq (e^{i\delta tJ\sum {\hat{\sigma }}_i^x{\hat{\sigma }}_{i+1}^x}{e}^{i\delta t{h}_T\sum {\hat{\sigma }}_i^z}{)}^{t/\delta t}+\mathcal{O}((\delta t)t).$$

$$U^{\text{ST}(2)}(t)=(e^{i\delta t\frac{h_T}{2}\sum {\hat{\sigma }}_i^z}{e}^{i\delta tJ\sum {\hat{\sigma }}_i^x{\hat{\sigma }}_{i+1}^x}{e}^{i\delta t\frac{h_T}{2}\sum {\hat{\sigma }}_i^z}{)}^{t/\delta t}$$

$$\begin{array}{rl}{H}_1& =-h_T\sum _{i=1}^N{\sigma }_i^z\\ H_2& =-J\sum _{i=1,3,\mathrm{5...}}^N{\sigma }_i^x{\sigma }_{i+1}^x\\ H_3& =-J\sum _{i=2,4,\mathrm{6...}}^N{\sigma }_i^x{\sigma }_{i+1}^x.\end{array}$$

$$d_l=7\ast N_t,$$

$$N_{\text{gates}}^{\text{pbc}}=2\ast (N_s\ast 2)\ast N_t$$

$$\begin{array}{rl}{N}_{\text{gates}}^{\text{obc}}& =(N_s)\ast 2-1)\ast N_t\\ & +((N_s-1)\ast 2)\ast N_t.\end{array}$$

$$\mathcal{F}(t;\delta t)=|⟨\psi (0)|(U^{\text{exact}}(t){)}^{†}{U}^{ST}(t;\delta t)|\psi (0)⟩|.$$

$$\mathcal{E}(\hat{\rho };p_x,p_y,p_z)=(1-p)\hat{\rho }+p_x{\hat{\sigma }}^x\hat{\rho }{\hat{\sigma }}^x+p_y{\hat{\sigma }}^y\hat{\rho }{\hat{\sigma }}^y+p_z{\hat{\sigma }}^z\hat{\rho }{\hat{\sigma }}^z.$$

$$⟨Z^{\text{exact}}⟩=\frac{⟨Z^{\text{noisy}}⟩}{p_{\text{measure}}}$$

$$⟨O⟩(0)=(⟨O(ϵ)⟩{)}^{\frac{r}{r-1}}(⟨O(ϵ)⟩{)}^{\frac{1}{1-r}}.$$

$$⟨\mathcal{O}(ϵ\ast r;t)⟩=A{B}^r+C.$$

$$\tilde{\chi}^{2}(t)=\frac{1}{N}\sum_{j=1}^{N}\Big{(}\frac{\langle n_{j}^{\text{exact}}(t)\rangle-\langle n_{j}^{\text{sim.}}(t)\rangle}{\delta\langle n_{j}^{\text{sim}}\rangle}\Big{)}^{2}.$$

$$⟨\mathcal{O}(ϵ\ast r;t)⟩=A{B}^r+C$$

Full text

I Introduction

There has been a fast growing interest for quantum computation in the context of high energy and nuclear physics Jordan et al. (2012); Martinez et al. (2016); Hamed Moosavian and Jordan (2018); Klco et al. (2018); Lamm and Lawrence (2018); Dumitrescu et al. (2018); Macridin et al. (2018); Raychowdhury and Stryker (2018); Stryker (2018); Yeter-Aydeniz et al. (2018); Hackett et al. (2018); Klco and Savage (2018); Roggero and Carlson (2018); Muschik et al. (2017); Kokail et al. (2018); Lu et al. (2018); Somma (2016). 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 111We currently are discussing the regime of 4 to 8 qubits.

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 Lamm and Lawrence (2018) or the Schwinger model Martinez et al. (2016); Klco et al. (2018).

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. Zhang et al. (2018). 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. Jordan et al. (2012). 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.

II The Quantum Ising Model

II.1 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 Fradkin and Susskind (1978); Kogut (1979). In this example, the classical model is the usual two-dimensional Ising model solved by Onsager Onsager (1944) and Kaufman Kaufman (1949). 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 $\hat{\sigma}^{x}$, governed by $J$, while the transverse magnetic field interactions (on-site coupling) governed by $h_{T}$ use $\hat{\sigma}^{z}$, where $\hat{\sigma}^{z}$ and $\hat{\sigma}^{x}$ are the traditional Pauli matrices. In the particle basis, we define the “particle number” at each $l$ site as

[TABLE]

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 example for $N_{s}=4$, the action of a sample operator on a sample state can be illustrated as

[TABLE]

We will see that the Hamiltonian only connects states for which the total particle number ($\hat{n}=\sum_{i}\hat{n}_{i}$) is the same modulo 2.

We define the Hamiltonian corresponding to open boundary conditions, hereinafter referred to as (OBC), as,

[TABLE]

The Hamiltonian corresponding to periodic boundary conditions (PBC) is defined as,

[TABLE]

while it is also interesting to consider a Hamiltonian with antiperiodic boundary conditions (ABC),

[TABLE]

For any of these Hamiltonians we define an operator to carry out the exact time evolution in units where $\hbar$ =1 as,

[TABLE]

II.2 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 $\hat{\sigma}^{z}$ at each site. This defines a unitary transformation that flips the sign of the operator $\hat{\sigma}^{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 number $\hat{n}$ defined above is conserved modulo 2.

When $N_{s}$ is even, it is also possible to invent a two-step transformation which changes the sign of the entire Hamiltonian. Similar equivalences appear for classical gauge theories Li and Meurice (2005). We first apply a $\hat{\sigma}^{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 $\hat{\sigma}^{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.

III Approximate evolution for $J<<h_{T}$

III.1 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\rangle$. 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 Kogut (1979). The model can also be solved exactly by performing a Wigner-Jordan transformation Kaufman (1949). However, at finite volume, boundary conditions should be treated carefully. To be more explicit, the term $a_{Ns}^{\dagger}a_{1}$ needs to be supplemented with a product of $\hat{\sigma}^{z}_{l}$ in order to reproduce the original spin Hamiltonian which requires a separate discussion for the even and odd sectors Kaufman (1949).

III.2 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. Katsura (1962); Nambu (1950) and the zero field case is examined in Refs. Lieb et al. (1961); Perk and Capel (1977); Perk (2017), time dependent $zz$ spin correlations were studied in Refs. Niemeijer (1967); Katsura et al. (1970) and $xx$ spin correlations were studied in Ref. Perk and Au-Yang (2009). If periodic boundary conditions are imposed, as in Eq. (4), Fourier modes diagonalize the perturbation. This lifts the degeneracy by a term proportional to $2J\text{cos}(2\pi m/N_{s})$. The perturbation also contains operators that connect to the 3-particle states; this leads to energy shifts $\mathcal{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 $|\psi\rangle$ and calculate $\langle\psi(t)|\hat{n}_{l}|\psi(t)\rangle$, where the calculations in the quantum mechanical approximation are relatively easy. For instance, for $|\psi(0)\rangle=|j\rangle$, we obtain

[TABLE]

where the “discrete” Bessel functions are defined as,

[TABLE]

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. Katsura et al. (1970). The approximation is accurate when $t$ is less than $\mathcal{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.

III.3 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,

[TABLE]

where $|ij\rangle=|0...0,1_{i},0...0,1_{j},0...\rangle$. 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 $\langle n_{i}(t)\rangle$ (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.

III.4 Finite volume corrections

It is possible to calculate the difference between the finite volume discrete Bessel functions $J_{n}^{(N_{s})}(x)$ and the usual infinite volume expressions $J_{n}(x)$. Using the Poisson formula, one finds that

[TABLE]

where the sum over $\ell$ runs over strictly positive and negative integers.The difference between $J_{n}^{(N_{s})}(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_{0}^{(8)}(x)$ and $J_{0}(x)$ becomes visible near $x\sim 4$. At that point the difference is almost saturated by the $\ell=\pm 1$ terms $J_{8}(x)+J_{-8}(x)$. The $\ell=\pm 2$ terms become important near $x\sim 12$.

IV Formulation of Real Time Evolution

For systems with a small number of spins, the Hilbert space of the model defined by Eq. (3) and Eq. (4) is small and this evolution operator is tractable for exact implementation on a classical computer; however, the Hilbert space of the model scales like $2^{N}$ for $N$ spins. For large $N$ ($N\gg 20$) it would only be possible to implement this on a quantum computer as the computational resources instead scale linearly.

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:

[TABLE]

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:

[TABLE]

While the estimated worse case error $\mathcal{O}(t~{}\delta t)$ is true in general, this bound over estimates the ST truncation error, which should be approximately $\mathcal{O}(Jh_{T}t\delta t)$ because the $\hat{\sigma}^{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,

[TABLE]

which can be found using the methodology proposed in Hatano and Suzuki (2005). 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 $\mathcal{O}(J^{3}t(\delta 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: $H=H_{1}+H_{2}+H_{3}$, where,

[TABLE]

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 Lloyd (1996) 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}$),

[TABLE]

and a total number of gate operations:

[TABLE]

[TABLE]

where $N_{s}$ is the number of sites and $N_{t}$ is the number of trotter steps.

V Results of Real Time Evolution

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 1.

V.1 Sources of error

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,

[TABLE]

For both OBC and PBC it appears that a trotter step of $\delta 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_{\text{1 qubit}}=1.0*10^{-4}$ and $p_{\text{2 qubit}}=5.0*10^{-4}$ Harty et al. (2014); Ballance et al. (2016)) and for slight improvements of current typical digital quantum computers IBM (2019); Rigetti (2019)

The Pauli error channel for single-qubit gates is defined in terms of the density matrix $\hat{\rho}$:

[TABLE]

The values $p_{x}$, $p_{y}$, and $p_{z}$ correspond to the probabilities of an $\sigma^{x}$, $\sigma^{y}$, or $\sigma^{z}$ error respectively occurring and $p=p_{x}+p_{y}+p_{z}$. The error channel for two qubit gates is given by $\mathcal{E}^{(2)}=\mathcal{E}\bigotimes\mathcal{E}$. The values of $p_{x},~{}p_{y},~{}p_{z}$ for the one and two-qubit gates are given in Tables 2 and 3. 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_{\text{measure}}$ given in Tables 2 and 3. Ref. Kandala et al. (2017) 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:

[TABLE]

V.2 Methods for dealing with error

Approximating the time evolution operator using the ST method introduces an error $\mathcal{O}(\delta t~{}t)$. It was suggested in Endo et al. (2018, 2018); Li and Benjamin (2017) 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. Endo et al. (2018) suggested several methods, one of which is an exponential extrapolation,

[TABLE]

In Eq. (21), $\epsilon$ 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 Endo et al. (2018, 2018); Li and Benjamin (2017) involves changing the extrapolation equation proposed in Eq. (21) to an exponential ansatz of the form:

[TABLE]

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. Endo et al. (2018). They argue that the algorithmic error rate, $\epsilon_{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.

V.3 Simulation Results

The initial state of the system for each case we looked at is given in Table 1. 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 $\chi^{2}$ value ($\tilde{\chi}^{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,

[TABLE]

These values over time are given in Table 4 and Table 5. The systematic errors from the ST approximation become significant at time scales on the order of $Jt\approx 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.

V.4 Results of simulated superconducting qubit quantum computer

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 $\sim$ $0.01$ for one qubit gates and $\sim 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 operators into our circuit. Similarities can be seen between Fig. 11 and Fig. 4 in Ref. Klco et al. (2018). The only difference in the procedure that we carried out is the extrapolation method that we used. Ref. Klco et al. (2018) used a quadratic ansatz while we used an exponential ansatz of the form:

[TABLE]

to extrapolate the noiseless observable. We used priors of $A=0.0\pm 0.5$, $B=0.0\pm 1.0$, and $C=0.5\pm 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.

V.5 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\delta 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.

Acknowledgements.

This work was supported in part by the U.S. Department of Energy (DOE) under Award Number DE-SC0019139. We thank the members of this grant for stimulating discussions. In particular we thanks Stephen Jordan, Nathalie Kclo, and Martin Savage for their input. We thank Jacques Perk for the suggestion of references regarding the XY-model. JUY was supported by the US Department of Energy (DOE), Office of Science, Office of High Energy Physics, under Award Numbers DE-SC0009998.