Simulating quantum many-body dynamics on a current digital quantum computer

TL;DR

This paper demonstrates the use of current IBM quantum computers to qualitatively simulate quantum many-body dynamics, revealing localization and interaction effects despite hardware limitations.

Contribution

It is the first to benchmark and showcase qualitative quantum many-body simulations on existing digital quantum computers.

Findings

Current quantum hardware cannot yet accurately simulate large-scale dynamics quantitatively.

Qualitative features of localization and interaction effects are observable.

Simulations are limited to small system sizes comparable to exact diagonalization.

Abstract

Universal quantum computers are potentially an ideal setting for simulating many-body quantum dynamics that is out of reach for classical digital computers. We use state-of-the-art IBM quantum computers to study paradigmatic examples of condensed matter physics -- we simulate the effects of disorder and interactions on quantum particle transport, as well as correlation and entanglement spreading. Our benchmark results show that the quality of the current machines is below what is necessary for quantitatively accurate continuous time dynamics of observables and reachable system sizes are small comparable to exact diagonalization. Despite this, we are successfully able to demonstrate clear qualitative behaviour associated with localization physics and many-body interaction effects.

Equations

$$⟨\psi (t)|{\hat{O}}_j|\psi (t)⟩,⟨\psi (t)|{\hat{O}}_j{\hat{O}}_k|\psi (t)⟩,$$

$$|\psi (t)⟩=e^{-i\hat{H}t}|\psi (0)⟩,$$

$$\hat{H}=-J\sum _{j=1}^{N-1}\left({\hat{\sigma }}_j^x{\hat{\sigma }}_{j+1}^x+{\hat{\sigma }}_j^y{\hat{\sigma }}_{j+1}^y\right)+U\sum _{j=1}^{N-1}{\hat{\sigma }}_j^z{\hat{\sigma }}_{j+1}^z+\sum _{j=1}^N{h}_j{\hat{\sigma }}_j^z,$$

$$M_j(t)=⟨\psi (t)|{\hat{\sigma }}_j^z|\psi (t)⟩.$$

$$N_{\text{half}}(t)=\sum _{j=1}^{N/2}⟨\psi (t)|\frac{{\hat{\sigma }}^z+1}{2}|\psi (t)⟩,$$

$$C_{jk}(t)=⟨\psi (t)|{\hat{\sigma }}_j^z{\hat{\sigma }}_k^z|\psi (t)⟩-⟨\psi (t)|{\hat{\sigma }}_j^z|\psi (t)⟩⟨\psi (t)|{\hat{\sigma }}_k^z|\psi (t)⟩.$$

$$F_Q(t)=\sum _{jk}{s}_j{s}_k⟨\psi (t)|{\hat{\sigma }}_j^z{\hat{\sigma }}_k^z|\psi (t)⟩-{\left(\sum _j{s}_j⟨\psi (t)|{\hat{\sigma }}_j^z|\psi (t)⟩\right)}^2,$$

$$S_{\text{vN}}\approx \frac{5}{32}{F}_Q,$$

$$M_{\text{GHZ}}=|⟨\hat{X}\hat{Y}\hat{Y}⟩+⟨\hat{Y}\hat{X}\hat{Y}⟩+⟨\hat{Y}\hat{Y}\hat{X}⟩-⟨\hat{X}\hat{X}\hat{X}⟩|\le 2,$$

$$e^{-i\hat{H}t}\approx \hat{U}(\Delta t{)}^M\hat{U}(\delta t),$$

$$e^{-i\hat{H}t}=e^{-i\hat{A}t}{e}^{-i\hat{B}t}+\mathcal{O}(t^2),$$

$$\hat{H}=-J\sum _{j=1}^{N-1}\left({\hat{\sigma }}_j^x{\hat{\sigma }}_{j+1}^x+{\hat{\sigma }}_j^y{\hat{\sigma }}_{j+1}^y\right)+U\sum _{j=1}^{N-1}{\hat{\sigma }}_j^z{\hat{\sigma }}_{j+1}^z+\sum _{j=1}^N{h}_j{\hat{\sigma }}_j^z,$$

$$e^{-i\hat{H}t}={\left(e^{-i\hat{H}\Delta t}\right)}^M,$$

$${\hat{A}}_j=e^{-i{h}_j{\hat{\sigma }}_j^z\Delta t},{\hat{B}}_j=e^{-i(U{\hat{\sigma }}_j^z{\hat{\sigma }}_{j+1}^z-J({\hat{\sigma }}_j^x{\hat{\sigma }}_{j+1}^x+{\hat{\sigma }}_j^y{\hat{\sigma }}_{j+1}^y))\Delta t}.$$

$$e^{-i\hat{H}\Delta t}=\left(\prod _j{\hat{A}}_j\right)\left(\prod _{j\text{ even}}{\hat{B}}_j\right)\left(\prod _{j\text{ odd}}{\hat{B}}_j\right)+\mathcal{O}((\Delta t{)}^2),$$

$$e^{-i\hat{H}t}=e^{-i\frac{t}{2}\hat{A}}{e}^{-it\hat{B}}{e}^{-i\frac{t}{2}\hat{A}}+\mathcal{O}(t^3).$$

$$\begin{array}{c}{\hat{A}}_j=e^{-i{h}_j{\hat{\sigma }}_j^z\frac{\Delta t}{2}},{\hat{B}}_j=e^{-i(U{\hat{\sigma }}_j^z{\hat{\sigma }}_{j+1}^z-J({\hat{\sigma }}_j^x{\hat{\sigma }}_{j+1}^x+{\hat{\sigma }}_j^y{\hat{\sigma }}_{j+1}^y))\frac{\Delta t}{2}},\\ {\hat{C}}_j=e^{-i(U{\hat{\sigma }}_j^z{\hat{\sigma }}_{j+1}^z-J({\hat{\sigma }}_j^x{\hat{\sigma }}_{j+1}^x+{\hat{\sigma }}_j^y{\hat{\sigma }}_{j+1}^y))\Delta t},\end{array}$$

$$e^{-i\hat{H}\Delta t}=\left(\prod _j{\hat{A}}_j\right)\left(\prod _{j\text{ even}}{\hat{B}}_j\right)\left(\prod _{j\text{ odd}}{\hat{C}}_j\right)\left(\prod _{j\text{ even}}{\hat{B}}_j\right)\left(\prod _j{\hat{A}}_j\right)$$

$$+\mathcal{O}((\Delta t{)}^3),$$

$$⟨\psi (t)|{\hat{\sigma }}_j^z|\psi (t)⟩=\sum _{\lambda :j=\uparrow }|{\alpha }_{\lambda }{|}^2-\sum _{{\lambda }^{\prime }:j=\downarrow }|{\alpha }_{{\lambda }^{\prime }}{|}^2,$$

Full text

Introduction

Quantum computers are general purpose devices that leverage quantum mechanical behaviour to outperform their classical counterparts by reducing the computational time and/or the required physical resources Nielsen2010. Excitement about quantum computation was initially fuelled by the prime-factorization algorithm developed by Shor Shor1995, which is most popularly associated with the ability to attack currently used cyber security protocols. More importantly, it provided a paradigmatic example of dramatic exponential improvement in computational speed when compared with classical algorithms. It has since been realised that the potential power of quantum computers could have far reaching applications, from quantum chemistry and the associated benefits for medicine and drug discovery Kandala2017, to quantum machine learning and artificial intelligence Biamonte2017. Even before reaching these lofty goals, there may also be practical uses for noisy intermediate-scale quantum (NISQ) devices Preskill2018.

These quantum devices can be implemented in a large number of ways, for example, using ultracold trapped ions Cirac1995; Martinez2016; Nam2019; Wright2019; Bruzewicz2019; Figgatt2018, cavity quantum electrodynamics (QED) Houck2012; Barends2014; Blais2007; Zhu2017, photonic circuits Aspuru2012; Sun2018; Tambasco2018, silicon quantum dots Kane1998; West2019; Yang2019, and theoretically even by braiding, as yet unobserved, exotic collective excitations called non-abelian anyons Nayak2008; Sarma2006; Lahtinen2017. One of the most promising approaches is using superconducting circuits Nakamura1999; Wendin2017; Krantz2019, where recent advances have resulted in devices consisting of up to 72 qubits, pushing us ever closer to realising so-called quantum supremacy Boixo2018. The apparent proximity of current devices to this milestone makes it timely to review the current capabilities and limitations of quantum computers.

Richard Feynman’s original idea was to simulate quantum many-body dynamics – a notoriously hard problem for a classical computer – by using another quantum system Feynman1982. Over the last couple of decades this approach of using a purpose built quantum simulator has been extremely successful in accessing physics beyond the reach of numerics on a classical computer. Most notable are cold atom experiments with optical lattices Bloch2008; Bakr2009; Schreiber2015; Choi2016; Bordia2017; Mitra2017; Cooper2018 where the natural evolution of the atoms corresponds to a high accuracy to that of a local Hamiltonian of choice. This has, for example, allowed us to study Hubbard model physics Bakr2009; Mitra2017 and many-body localized systems Choi2016; Schreiber2015; Bordia2017 in two dimensions. More recently, there have also been advances in trapped ion quantum simulators, which have the benefit of being able to implement long range interactions Cirac1995; Martinez2016; Garttner2017, and have been used to study the Schwinger-mechanism of pair production/annihilation in 1D lattice QED Martinez2016 and Floquet time crystals Zhang2017.

Universal quantum computers are also increasingly looking like a feasible setting for simulating quantum dynamics. One of the biggest advantages of using a quantum computer for this purpose is the flexibility it offers. A single quantum device could in principle perform simulations that currently require several different experiments, using disparate methods. Furthermore, it should be possible to access new physics not currently accessible, most notably, the simulation of lattice gauge theories with dynamical gauge fields. These are ubiquitous in the theoretical study of strongly-correlated quantum matter but require multi-body couplings, which have so far proven difficult to achieve in experiment Wiese2013; Zohar2016.

A particularly exciting opportunity has been provided by IBM in the form of an online quantum computing network called IBM Q. This consists of a set of small quantum computers of 5 and 16 qubits that are availably freely to the public, two 20 qubit machines accessible by IBM Q partners IBMQ, and the qiskit python API qiskit for programming the devices. The publicly available resources have already resulted in a spread of results, such as calculating the ground state of simple molecules OMalley2016; Kandala2017, creating and measuring highly entangled many qubit states Wang2018; Choo2018, implementing quantum algorithms Bravo-Prieto2019; Doronin2019; Amico2019, and simulating non-equilibrium dynamics in the transverse-field Ising Zhukov2018; Cervera-Lierta2018, Heisenberg Lamm2018 and Schwinger Klco2018 models, as just a few examples. Given the infancy of quantum computing efforts, all these results are understandably small scale and of limited accuracy.

IBM is not alone in their efforts to make quantum computing more mainstream, with Microsoft introducing the $Q$-sharp programming language q-sharp, Google developing the Circq python library circq, and Rigetti providing their own Quantum Cloud Service and Forest SDK built on Python pyquil. Rigetti and IonQ also provide selective public access to hardware – based on superconducting qubits and trapped ions, respectively. All of these resources are allowing a lot of hands on experience with quantum computers from researchers and the general public all around the world. Furthermore, it highlights that there are currently several parallel efforts of research and development from industry focussed on quantum computing, quantum programming and cloud-based services that are flexible and prepared for future hardware.

Given the current stage of development, and the immense expectation from the public and physics communities, it is timely to critically assess and benchmark the state-of-the-art. In this article we consider the far-from-equilibrium dynamics of global quantum quenches simulated on an IBM 20 qubit quantum computer. The models that we consider are of central importance to condensed matter physics and display a wide range of phenomenologies. By measuring a range of physical correlators, we can assess the capabilities and limitations of current quantum computers for simulating quantum dynamics.

Results

**Setup

**In this article we study global quantum quenches Essler2016; Vasseur2016, that is we calculate local observables and correlators of the form

[TABLE]

where $\hat{O}_{j}$ are local operators and the time-dependent states are

[TABLE]

and where $|\psi(0)\rangle$ is the initial state, which differs globally from an eigenstate of $\hat{H}$. In other words, we can consider $|\psi(0)\rangle$ to be the ground state of a (time-independent) preparation Hamiltonian, $\hat{H}_{0}=\hat{H}-h\hat{V}$, where $\hat{V}$ is a global perturbation, and the perturbation strength $h$ is instantaneously quenched from zero at time $t=0$.

We consider one dimensional spin-1/2 chains, consisting of $N$ spins, initially prepared in either a domain wall configuration, $|\cdots\!\downarrow\downarrow\downarrow\uparrow\uparrow\uparrow\!\cdots\rangle$, or Néel state, $|\cdots\!\uparrow\downarrow\uparrow\downarrow\uparrow\!\cdots\rangle$. Quenches from these initial states are particularly easy to prepare due to the local tensor-product structure of the initial state and have been studied in cold atom experiments Schreiber2015. The time evolution after the quantum quench will be governed by a Hamiltonian of the form

[TABLE]

with $J>0$, and $\hat{\sigma}^{\alpha}_{j}$ are the Pauli matrices with eigenvalues $\pm 1$. This model can also be written in terms of hard-core boson, or spinless fermion Hamiltonian via the Jordan-Wigner transformation.

We will consider four distinct cases for the Hamiltonian parameters:

• (i)

The XX chain, defined by $U=0$ and $h_{j}=0$. This is a uniform, non-interacting model.

• (ii)

Disordered XX chain for $U=0$ and $h_{j}$ uniformly randomly sampled from the interval $[-h,h]$. This is the prototypical lattice model for Anderson localization Anderson1958.

• (iii)

XXZ spin chain, defined by $U\neq 0$ and $h_{j}=0$. This is a particular case of the Heisenberg model for quantum magnetism.

• (iv)

XXZ chain with linear potential, defined by $U>0$ and $h_{j}=hj$, i.e., a linearly increasing potential. This model was recently studied in the context of many-body localization without disorder Schulz2018; VanNieuwenburg, where Wannier-Stark localization Wannier1962 is induced by the linear potential.

This family of Hamiltonians covers a large range of physics in condensed matter, from the integrable limit Essler2016, to many-body quantum magnetism Vasiliev2018, to that of many-body localization Nandkishore2015; Abanin2017. These Hamiltonians are perfect testbeds for digital quantum simulation since they can be directly simulated on a quantum computer – the spins-1/2 of the former correspond directly to the qubits of the latter, and the local connectivity of the qubits is suited for the local form of the Hamiltonian.

In our simulations we achieve the time evolution using a Trotter decomposition of the unitary time evolution operator $\hat{U}(t)=e^{-i\hat{H}t}$. Our simulation proceeds by the following steps. First, we create the initial state, which is done by applying Pauli-X operations to the default initial state of the IBM machine, $|\cdots\uparrow\uparrow\uparrow\cdots\rangle$. Second, we discretize time and the evolution becomes a product of discrete evolution operations $\hat{U}(t)=\hat{U}(\Delta t)\cdots\hat{U}(\Delta t)$. Note that in our plots we include additional data points by varying $\delta t$ in the final discrete evolution step. Third, we trotterize the discrete operators $U(\Delta t)$ into a product of one- and two-qubit unitaries. Finally we decompose these unitaries into the CNOT and single qubit rotation gates that can be directly applied on the IBM devices. Please see Fig. 1 for a schematic of this procedure, and see Methods for more details.

Once we have constructed the state $|\psi(t)\rangle$, we then perform measurements, which allows us to construct the quantities of interest in Eq. (1). We will only consider correlators of $\hat{\sigma}^{z}_{j}$ operators, where we only need to make measurements in a single ($z$-basis) for the spins. In the following we will use 8192 measurements per data point, which means that the statistical error for these local correlators is $\sim 0.01$, which is too small to be included in our figures.

We consider three different types of dynamical quantities:

• •

The local magnetization

[TABLE]

In the case of a domain wall initial state, we will also compute

[TABLE]

which is equal to zero at $t=0$ and grows as the domain wall spreads.

• •

The connected equal-time correlator

[TABLE]

Note, the connected form of this correlator measures the quantum correlators between distant spins but is not sensitive to the classical correlations in our initial states.

• •

The quantum Fisher information (QFI). We will consider a particular case of the QFI for a pure state, namely,

[TABLE]

where $s_{j}=+1$ for left half of the sites $j$ and $s_{j}=-1$ for the right half at $t=0$. More generally, the QFI (for a pure state) is defined as the variance, $4\left(\langle\hat{O}^{2}\rangle-\langle\hat{O}\rangle^{2}\right)$, where $\hat{O}$ is a sum of local operators, which each have a spectrum of unit width. In our case we have $\hat{O}=\frac{1}{2}\sum_{j}s_{j}\hat{\sigma}^{z}_{j}$. The QFI is an entanglement witness Hyllus2012; Toth2012; Hauke2016, and for our chosen definition in Eq. (7) is also closely related to the von Neumann entanglement entropy Klich2009; Song2011; Song2012.

**The IBM Quantum Computers

**The quantum computer that we use is the latest 20-qubit IBM device, codenamed ibmq_poughkeepsie. It consists of a two-dimensional array of qubits that have local connectivity. We can perform arbitrary single qubit rotations and controlled-NOT (CNOT) gates between connected qubits, see Methods. For the data presented in this paper the average readout errors, CNOT errors, and T2 (dephasing) times were approximately $4\%,2\%$ and $90\mu$s, respectively. An important point to note is that the IBM machines are recalibrated on an approximately daily basis, which means the data can vary across days. Crucially, we find that the our results are qualitatively reproducible, and we compare data obtained across three consecutive days in the Methods.

To benchmark the accuracy of the simulation, we compare the data with a numerical implementation of the Trotter evolution, as well as continuous-time exact diagonalization (ED). The errors in our results are strongly influenced by the number of CNOT gates in the corresponding quantum circuit. One of the reasons for this is that the fidelities of these two qubit gates are an order of magnitude worse than the single qubit gates. The CNOTs also take a longer real-world time to implement than the single qubit gates. The increased implementation time of the circuits increases the potential for errors due to energy relaxation and dephasing, parametrized by the T1 and T2 times, as well as other environmental effects and cross-talk. On the IBM device, we use $N=6,8,10$ of the qubits as our system with 4 or 5 symmetric Trotter steps. These qubits are chosen as the connected subset with the lowest average CNOT errors such that the single qubit measurement errors and T2 decoherence times do not exceed a certain threshold. Please see Methods for more details about the quantum device and details of the algorithm used to select the qubits. All data presented below is available in Ref. IBMdata.

**Local Magnetization

**First, we consider results for the uniform XX spin chain Hamiltonian with $U=0$ and $h_{j}=0$ (case (i)), quenched from a domain wall configuration, shown in Fig. 2. Figure 2(a-c) show a comparison of the magnetization of the fourth (middle) and sixth (end) spins of the chain as computed by exact diagonalization with continuous time evolution, a numerical implementation of the Trotter decomposition and the corresponding data from the IBM machine. The data from the machine is further split into the raw data (orange triangles) and constrained data (red squares). The constrained data only considers those measurement outcomes that have the same total magnetization as the initial state – which is a conserved quantity – that is, we restrict to the physical Hilbert space of the Hamiltonian we are simulating. We discuss this rudimentary error mitigation method further in the context of quantifying the accuracy of the quantum computer (see Fig. 6(b)), and it will be used in all subsequent figures. See Methods for more details.

The data in Fig. 2(a-b) show that while all curves show reasonable agreement at short times – for instance, we have a delayed decay of the magnetization in Fig. 2(b) – the accuracy of the IBM data becomes bad very quickly. For times $Jt>1$ the magnetization as measured on the machine approaches zero, which is the expected average value if the system thermalizes or if we were to randomly sample states. The agreement between the numerical results obtained from exact diagonalization (ED) and trotterization shows that the inaccuracy of the results is due to the machine and not our approximation of the evolution. This rapid decay in the accuracy with number of Trotter steps was also observed in Ref. Zhukov2018 for the transverse field Ising model on the 5 and 16 qubit IBM machines Leeb2018.

In Fig. 2(c) we consider the evolution over a longer time window, up to $Jt=10$, rather than $Jt=2.5$. We see that the data very quickly approaches [math] and remains there, indicating that the system has equilibrated. From this figure we can see that due to the quality of the machine we want to take $\Delta t$ in our Trotter steps as large as possible so that we can reach longer times without thermalizing. However, for the Trotter evolution to be accurate we want $\Delta t$ as small as possible. The $\Delta t$ that we use are chosen (without much optimization) to maintain a reasonable accuracy in both the IBM data and the numerical trotterization.

Next, we consider the site dependence of the magnetization shown in Figs. 2(d-f), for $N=6,8,10$ sites and compared with the numerically computed Trotter evolution shown in Fig. 2(g). In these figures we see a clear qualitative agreement between the experimental and numerical results at short times, particularly by the presence of a linear light-cone causality structure for the spreading of the domain wall. This qualitative agreement, however, also worsens at longer times, and as we increase the system size. In particular, in Fig. 2(d) we can see a marked decrease in accuracy every three time steps. The origin can be explained as follows. Each block of three time steps is computed for the same number of time steps with $\Delta t$ varying in the final Trotter steps (as explained earlier). It is when we add an additional time step for the next block – and thus increase the number of gates in the quantum circuit – that we see a drop in the accuracy. This behaviour is also seen in Fig. 6(b) where we show the number of measurements that are in the physical Hilbert space. There is a clear decrease in the percentage after the introduction of each new Trotter step.

**Short-Time Many-Body Physics

**While the results in Fig. 2 may at first seem discouraging for quantitative large-scale dynamical simulations, we will show in this section that we may still observe qualitative behaviour associated with non-trivial quantum phenomena. Here, we consider $N_{\text{half}}(t)$, defined in Eq. (5), after quenching from the domain wall initial state for three different cases.

First, we consider the disordered XX chain (case (ii)), which is known to exhibit Anderson Localization in 1D for all values of the disorder strength, $h$ Kramer1993. As a consequence of increasing the disorder strength, the extent of the spreading of the domain wall, and consequently the growth of $N_{\text{half}}$ is reduced (indicated by a black arrow). We are able to reproduce this behaviour qualitatively for short times on the IBM machine as shown in Fig. 3(a). The corresponding numerical Trotter results are shown in the inset, which shows that while the accuracy of the results is quite low, the qualitative behaviour is still captured. Once again, we see that the data is biased towards the scrambled value as the number of Trotter steps is increased, which in the present case is $1.5$.

Next, we consider the uniform XXZ chain (case (iii)) with $U>0$ and $h_{j}=0$. In contrast to the previous cases, this Hamiltonian is interacting and thus describes true many-body physics. At short times, the spreading of the domain wall is also hindered by the energy cost of the additional interactions between neighbouring spins. Once again, we can see this behaviour qualitatively in the experimental results, shown in Fig. 3(b). Note that while the short time results are similar to the previous case, it is due to a different physical mechanism and is a many-body effect. The long-time behaviour would be starkly different from that of the Anderson localized case, which is, however, beyond the current capabilities of the IBM quantum computers.

In the third case we combine features of the previous two models and include both on-site potential energies and interactions. We will consider the XXZ spin chain with a linear potential (case (iv)), and both $U>0$ and $h>0$. If we compare with $U=0$, that is, a linear potential alone, then the eigenstates will be localized Schulz2018; Wannier1962; VanNieuwenburg, and thus the spreading of the domain wall will be limited. If we have $U>0$, then the two energy costs can compensate. For instance, consider flipping the middle two spins, then there will be an increase in energy due to the potential but a decrease in the interaction energy. Therefore, the presence of interactions makes it easier for the domain wall to spread resulting in an increase of $N_{\text{half}}$ as a function of $U$ (see black arrow). This simple argument of energetics is confirmed by the numerical results in the inset of Fig. 3(c), and is qualitatively reproduced in the experimental data shown in the main figure. Note, however, that this trend is less pronounced than in the previous two cases. This is in part due to the small scale of the changes in the exact results, as well as the bias towards the thermal value of 1.5 at long times.

In this data, the addition of disorder and interactions lead to similar qualitative behaviour on the time scales that we have considered. However, at longer times there is a clear difference between the two cases. In the former we have localization behaviour leading to the long-time persistence of the initial imbalance, whereas interactions generically lead to ergodic and thermalizing behaviour, resulting in the loss of this information in local observables. With the current devices we are unfortunately unable to distinguish these two different regimes.

**Spreading of Correlations – Light Cone

**In the previous section, all the data correspond to some combination of local average magnetizations. Here, we look at the correlations between pairs of spatially separated spins, and how these correlations spread. Lieb and Robinson showed that for a local Hamiltonian, correlations can spread at most linearly and display a light-cone causality structure Lieb1972. For instance, for tight-binding spinless fermions, correlations spread with a speed proportional to their maximum group velocity, $v=4J$ Essler2016.

In Fig. 4 we compare ED, numerical trotterization and experimental results for the connected spin correlator defined in Eq. (6), for $h_{j}=U=0$. We measure correlations between the first and the $j^{\text{th}}$ spin after a quench from a charge density wave initial state. The ED results show a clear ballistic spreading of correlations, which is also qualitatively reproduced by the numerical trotterization and the IBM data. As with all previous results, the agreement between the numerical and IBM results is best at short times, but the IBM data is able to capture the point where the correlations reach the system size. Furthermore, for this free model, there are only significant correlations along the light-cone and not within it, which is also approximately captured by the experimental data. There does, however, seem to be a slightly faster light-cone velocity, which is most evident in the shift of the position of the peak on site 6. This indicates that there is an effective renormalization of the Hamiltonian parameters, particularly $J$, due to the errors in the machine.

**Quantum Fisher Information

**Finally, we consider the quantum Fisher information defined in Eq. (7) starting from both the Néel and domain wall initial states, shown in Figs. 5(a) and (b) respectively. In this section we will consider the model of (case (i)), for which the quantum Fisher information has two important relations to entanglement, which we outline below. Note that the quantum Fisher information has a more general definition for mixed states, which reduces to our definition for pure states. We do not consider the more general definition since we are simulating unitary evolution from a pure quantum state and we are not able to differentiate between the unitary and non-unitary errors occurring in our circuits.

For non-interacting models (as is the case for the Hamiltonian of (case (i))), there is a direct relationship between the bipartite von Neumann entanglement entropy and the magnetization fluctuations Klich2009; Song2011; Song2012. The former is defined by $S_{\text{vN}}=-\text{Tr}[\rho^{A}\ln\rho^{A}]$, where $\rho^{A}$ is the reduced density matrix for half of the system. The variance of the half chain magnetization is proportional to our definition of the QFI and we have the approximate relation

[TABLE]

see Ref. Song2012 for details of this cumulant expansion. In MBL systems the QFI – using instead the staggered magnetization for the operator $\hat{O}$ in Eq. (7) – also appears to mimic the bipartite entanglement entropy and grows logarithmically after a quantum quench Smith2016a.

The QFI is also a multi-partite entanglement witness Pezze2009; Toth2012; Hyllus2012, and in particular, if the state is separable then we have that $F_{Q}\leq N$, which is known as the shot noise limit in quantum metrology Helstrom1976; Holevo1982; Giovannetti2004. For a general entangled state, however, the QFI is bounded by $F_{Q}\leq N^{2}$, and a value of $F_{Q}/N\geq m$, indicates at least $m+1$-partite entanglement. Note that the QFI is sensitive to the choice of the operator $\hat{O}$, and for example, if we choose $\hat{O}\propto\sum_{j}\hat{\sigma}^{z}_{j}$, the total magnetization, then $F_{Q}=0$, since our models conserve the total magnetization.

Let us now consider the IBM results for the QFI, starting with a quench from a Néel initial state shown in Fig. 5(a). We first note that the numerical results for the QFI (black) closely follow the bipartite von Neumann entanglement entropy (blue). The results from the IBM machine are also able to reproduce this behaviour quite accurately, characterized by the linear growth with a maximum just after $Jt=1$ due finite size effects. The fact that we measure $F_{Q}/N>1$ also implies entanglement in the state on the IBM quantum computer.

In Fig. 5(b) we consider a quench from a domain wall. In this case both the numerical and experimental simulations have the same qualitative behaviour but are further off in absolute value, as compared with Fig. 5(a). The QFI computed on the IBM machine once again is able to reproduce the behaviour of the von Neummann entanglement entropy.

The entanglement entropy is generally a difficult quantity to measure. It typically requires some form of state tomography, which consists of a set of measurements in a number of different bases that grows exponentially with system size, rendering it impractical for large systems. Furthermore, even with low error rates, the resulting density matrix may be unphysical Smolin2012; Choo2018. The quantum Fisher information may provide an alternative in certain circumstances because it can be significantly easier to compute – for the definition that we consider we only need to measure in a single basis.

**Quantifying the Accuracy of the Quantum Computer

**While in the previous sections we showed that the current IBM device can qualitatively reproduce physical behaviour, it is also important to develop practical quantitative measures of their accuracy. These measures should allow us to track the evolution of the quantum computers as they are developed and improved, as well as potentially providing further insight into the various sources of error that are present.

Going beyond the reported gate errors, one of the simplest things to measure is the violation of the Mermin inequality Mermin1990,

[TABLE]

where $\hat{X},\hat{Y}$ are the $x$ and $y$ Pauli-operators, and we omit the consecutive site labels on the operators. This inequality should hold for a classical theory with local realism. If we consider the GHZ state $|\psi\rangle=\frac{1}{\sqrt{2}}(|\!\uparrow\uparrow\uparrow\rangle-|\!\downarrow\downarrow\downarrow\rangle)$, then $M_{\text{GHZ}}=4$ and this bound is maximally violated and can be used to demonstrate the "quantumness" of the machine. In Fig. 6(a) we show the values computed across 4 consecutive days for the first three qubits used for the $N=6$ simulations. This data shows that we are consistently able to violate the Mermin inequality.

As can be seen in the methods, the fidelities of individual gates do not necessarily reflect the accuracy of a simulation across many qubits. It may therefore make sense to consider more practical measures of quality that more directly relate to the simulations we are performing. We use the physical conservation laws of the evolution to improve the accuracy of our simulations by throwing away measurements of unphysical states. The number of measurements that are kept/thrown away can also be taken as a quantitative measure of the effective accuracy of the machine since for a perfect quantum computer we should find that all measured states are in the physical Hilbert. In Fig. 6(b) we show the percentage of measured states that are physical for 4 consecutive days, showing the variation in the effective quality of the device.

As an unbiased measure of the quality of our simulation we suggest to compute the identity in the form $\hat{I}=\sqrt{\hat{U}^{-1}(t)\hat{U}(t)}$. By implementing a circuit to perform the forward and backward time evolution, the probability of returning to the initial state measures the accuracy of the implementation of the unitary $\hat{U}(t)$. If the implementation is non-unitary then we should expect decay with an increasing number of Trotter steps. We should also expect decay of this quantity if there are only unitary errors since this quantity takes the form of a Loschmidt echo, which measures the sensitivity of the system to perturbations. It does not rely on any special properties of $\hat{U}(t)$, such as the presence and knowledge of conserved quantities. It therefore provides an unbiased measure of the quality of the implementation of $\hat{U}(t)$ and of the quantum device.

We show the computation of the identity performed on the IBM machine across four consecutive days in Fig. 6(c). In an ideal machine this quantity should be identically equal to 1 for all times. However, with each additional Trotter step – corresponding to the observed plateaux – the accuracy drops significantly. We also note that the behaviour observed in Figs. 6(b) and (c) are very similar, and both reflect the accuracy of the simulations in the previous sections.

Discussion

Probably the most striking feature of our results from the IBM machines is the low quantitative accuracy when compared with exact numerics. Considering the limited system sizes and time scales that we can reach; it highlights the current limitations of these quantum devices. Most importantly, this shows that whilst the number of qubits is now reaching limits beyond the capabilities of classical computers, the error rates and/or the isolation of these quantum computers is not yet sufficient for useful computations, at least for accurately simulating quantum dynamics. Due to the large array of possible error sources, pinning down the most damaging for our simulations is a difficult task, and is an import practical area for future research.

Despite this unfortunate conclusion, we are still able to access a range of qualitative physical behaviours demonstrating non-trivial simulations of quantum dynamics. We were able to compute a range of expectation values and two-point correlators, and observe behaviour associated with localization, many-body interactions, and the ballistic spreading of quantum correlations. We also observed the compensation of different energy costs due to on-site potentials and neighbouring site interactions and witnessed the generation of entanglement due to unitary evolution through the QFI.

The goal of using quantum computers for dynamical simulations is to be able to access systems intractable using classical algorithms. The limitations of the current devices that we have observed in our results demonstrate the need for improved quality of the machines and not simply adding more qubits. In this type of simulation, the number of gates, and thus the real execution time of the quantum circuits, grows linearly with the system size and with the number of Trotter steps. This means that we can estimate that we would need at least an order of magnitude improvement in a combination of the gate fidelities and/or T1/T2 times to get close to achieving this goal.

One of the biggest challenges facing the field of quantum computation is how to deal with errors. Although this is primarily an engineering issue to increase the quality, isolation and control of the devices, there is also the theoretical contribution of error correction methods Nielsen2010; PreskillChapt7; Devitt2013. A particularly promising avenue for error correction is to use surface codes Fowler2012; Bravyi2017. One big advantage of these methods is the moderately low fidelities required for them to work effectively. As the size and quality of the quantum computers increases, it is hoped that these error correction schemes will allow us to rapidly increase their scale, and with it their utility. In the meantime, there may also be more room for practical error mitigation schemes, such as the one that we have used, to get the most out of NISQ devices.

While we have considered a range of correlation functions and physical mechanisms, there is still much that can be learnt about the current quantum computers and how well they can simulate quantum dynamics for condensed matter systems. In particular, there are physical mechanisms beyond those that we have considered here. For instance, models with gauge coupling to dynamical gauge fields Martinez2016; Klco2018; Wiese2013; Zohar2016, and the physics of confinement Kormos2016; Liu2019. In these settings there is also hope that interesting physics can be extracted from the short-time dynamics and thus may be suited to the current machines. The combination of disorder and interactions, resulting in the many-body localized phase, is also currently a particularly active area of research Abanin2017; Nandkishore2015; Choi2016; Schreiber2015. The investigation of the transition between MBL and ergodic dynamics may also benefit in the future from quantum computation. It is notoriously difficult to study numerically due to the requirement of large systems and/or long-time simulations Pal2010; Luitz2014; Herviou2018.

In conclusion, digital quantum simulation is still in its infancy and we have shown that it requires an order of magnitude improvement in fidelity and coherence until it will realistically outperform classical computers, at least applied to dynamical problems of interest in condensed matter physics. However, while it is hard to predict the pace of technological progress, our results will serve as a useful benchmark for improvements in the foreseeable future; and in the long run they will provide a snapshot of capabilities at the beginning of a new quantum simulation era.

METHODS

**Implementation: Trotterized Evolution

**For our global quench protocol, we first need to prepare the initial state of our system. Since the IBM quantum computers are initialized in the state $|\!\uparrow\uparrow\!\cdots\,\rangle$ by default, both of our choices of initial states are tensor product states in the $z$-basis and thus can be created by applications of the Pauli $X$ gate. Next, we need to implement the time evolution, which proceeds by three main steps that we will briefly outline here.

First, we discretize time, that is, we split the time evolution operator, $\hat{U}(t)=e^{-i\hat{H}t}$, into a sequence of discrete operators, i.e., $\hat{U}(t)=\hat{U}(\Delta t)\cdots\hat{U}(\Delta t)$ with fixed $\Delta t$. Each application of the discrete operator $\hat{U}(\Delta t)$ is called a Trotter step. This is illustrated in Fig. 1(a) where we apply an increasing number of Trotter steps to reach later times.

In our simulations we fix $\Delta t$ to fix the accuracy of our approximation. However, since we can only implement up to $5$ Trotter steps, we add additional data points by varying $\Delta t$ in the final Trotter step. More explictly, consider Trotter steps $M$ and $M+1$, we can add additional data points by using the evolution operator

[TABLE]

where $\delta t=\Delta t/r$, where $r$ is the number of data points we want between $t=M\Delta t$ and $t=(M+1)\Delta t$. Since the accuracy of the trotter decomposition $U(\delta t)$ is better than that of $U(\Delta t)$ (see below), these extra data points will have errors intermediate between that of the $M^{\text{th}}$ and $(M+1)^{\text{th}}$ steps.

Second, we perform a Trotter decomposition of (trotterize) the unitary evolutions, that is, we approximately decompose the operator $\hat{U}(\Delta t)$ into a sequence of unitaries that act on at most two neighbouring qubits. In the following we use either a basic or symmetric Trotter decomposition, shown in Fig. 1(b) and (c), respectively. For $m$ Trotter steps of length $\Delta t$, the error of the symmetric decomposition is of order $\mathcal{O}(m(\Delta t)^{3})$, compared with $\mathcal{O}(m(\Delta t)^{2})$ for the basic decomposition. Note, however, that due to the symmetric structure, when we apply several Trotter steps we can combine several layers of gates. This means that we only need an extra two layers of gates compared with the basic decomposition regardless of the number of Trotter steps.

Third, we must efficiently decompose these two qubit operators into the gates that can be directly implemented on the quantum device, which are the CNOT gate and arbitrary single qubit unitaries. An efficient decomposition is found in Ref. Vatan2004a, which we summarise below. The result is that if $U\neq 0$ then $\hat{B}$ and $\hat{C}$ (defined in the figure caption) can be implemented using three CNOTs and with $U=0$ this can be reduced to only two CNOTs.

**Trotter Decomposition

**In this paper we use a Trotter decomposition (commonly known as a trotterization) of the unitary time evolution operator $\hat{U}(t)=e^{-i\hat{H}t}$. That is, we want to approximate these operators by a sequence of more easily implemented operators, namely those that act on at most two qubits.

As a starting point, consider a Hamiltonian of the form $\hat{H}=\hat{A}+\hat{B}$, where $[\hat{A},\hat{B}]\neq 0$. Then, since these operators do not commute in general, we have that

[TABLE]

which can be naturally extended to Hamiltonians that are sums of more than two terms. To use this fact for trotterized evolution we can use the two steps outlined in the main text. We will now go through the details of these steps in more detail, for the Hamiltonian Eq. (3):

[TABLE]

which we rewrite here for convenience.

First, we discretize time and split the evolution operator $\hat{U}(t)$ into a product of discrete evolution operators $\hat{U}(\Delta t)$, i.e.,

[TABLE]

where $\Delta t=t/M$ and $M$ is the number of "Trotter steps". Since the Hamiltonian commutes with itself, Eq. (13) is exact. To perform time evolution, we will typically fix $\Delta t$ and increase the number of Trotter steps $M$ to reach later times.

The second step is to approximate each of these discrete time operators in a similar manner to Eq. (11). For notational simplicity, let us first define the operators

[TABLE]

Using these operators, we can make the approximation

[TABLE]

which corresponds to the schematic quantum circuit show in Fig. 1(b) in the main text, and which we will refer to as the basic trotterization. If we wish to evolve to time $t=M\Delta t$, then we find that the error is $\mathcal{O}(\Delta t)$, which is controlled by the size of the Trotter step $\Delta t$. We can therefore improve the accuracy of the approximation by decreasing $\Delta t$, however, this must be balanced against the cost of needing more Trotter steps, as explained in the main text.

To improve the accuracy of our simulations we can use better approximations to the discrete evolution operators by way of higher-order Trotter decompositions Hatano2005. The leading error term in Eq. (11) is due to the non-zero commutator $[\hat{A},\hat{B}]$. By compensating for this error, we can increase the order of the leading error term.

The only higher-order decomposition that we will consider is the symmetrized Trotter step. Let us again start with the simple case of $\hat{H}=\hat{A}+\hat{B}$. The symmetric decomposition would then be

[TABLE]

The error in this decomposition is of higher-order due to the symmetry which ensures that $U(-t)=U(t)^{\dagger}=U^{-1}(t)$. This means that the even-order error terms vanish and the leading error is $\sim(\Delta t)^{3}$. See Ref. Hatano2005 for more details and for an iterative method for constructing higher-order decompositions.

For the Hamiltonian Eq. (3) in question, let us again make some definitions to simplify notation

[TABLE]

then the symmetric Trotter decomposition is

[TABLE]

which is shown schematically in Fig. 1(c) of the main text.

**Measurement

**When making a measurement we will find the system in one of the many-body states $|\lambda\rangle$ in this basis, e.g., $|\!\uparrow\downarrow\downarrow\cdots\rangle$ or $|\!\downarrow\uparrow\downarrow\cdots\rangle$. By performing multiple runs and measurements, we can approximate the probability of measuring each of the many-body states, that is, we can extract $|\alpha_{\lambda}|^{2}$, where $\alpha_{\lambda}$ is the probability amplitude for the state $|\lambda\rangle$. These probabilities can then be used to construct the observables. To be more concrete, consider the expectation value of the operator $\hat{\sigma}^{z}_{j}$ on site $j$. This is computed as follows,

[TABLE]

where the sums are over all states with the $j^{\text{th}}$ spin up or down respectively, and $|\alpha_{\lambda}|^{2}$ is the proportion of measurements for which we found the state $\lambda$. In all the following experiments we will use 8192 measurements per data point, which means that the statistical error for these local correlators is $\sim 0.01$, which is too small to be included in our figures.

**Error Mitigation: Physical Hilbert Space

**As we noted in the previous section, due to the presence of conserved quantities in the Hamiltonian evolution, the Hilbert space of states splits into those that are physically allowed by the evolution and those that are not. In particular, the models we consider have conserved net magnetization $S_{z}=\sum_{j}\hat{\sigma}_{j}^{z}$. This fact turns out to be advantageous, and allows us to perform rudimentary error mitigation – a term that we use to distinguish it from scalable error correction methods, since in this case we deal only with the lowest order errors. The idea is to simply disregard any measurements for states outside of the physical Hilbert space. Let use present a simplified argument for why this might be a good thing to do, where we will first assume that only bit-flips can occur.

A single bit-flip in the course of the evolution will take us outside of the Hilbert space, and let us denote the probability of this happening as $\Delta$. However, the lowest order of errors within the physical Hilbert space is $\Delta^{2}$, i.e., we need at least two bit-flip errors to get back to the same total magnetization. Hence, by discarding counts outside of the physical Hilbert space we reduce the leading order error to $\Delta^{2}$. If the probability, $\Delta$, is sufficiently small, then we can rely on these perturbative arguments. However, if the error rate is large enough, then multiple bit-flip errors can become significant and the error mitigation will be ineffective.

In Fig. 6(b) we show the percentage of measurements that are within the physical Hilbert space, that is, the percentage of measurements that are kept. For the first data points, we are simply measuring the initial state, which has an average measurement fidelity per qubit of $\sim 95\%$ leading to a value of approximately $(95\%)^{6}\approx 70\%$. At long times the number of retained states stabilises and approaches a value that corresponds to the percentage of states in the total Hilbert space that are physical – in other words, the percentage probability that a randomly selected state is in the physical subspace. Once this number of discarded measurements is reached, the errors are very large and the error mitigation is no longer effective.

Bit-flips are not the only type of errors that could occur. As an example, there are also phase errors, which do not necessarily change the net magnetization. However, we note that the constrained data does typically have improved accuracy, and so we use the restriction to the physical Hilbert space for all our subsequent data.

**Quantum Circuits

**Here, we will go over some of the details necessary to implement the trotterized evolution operators on the IBM devices. We will only cover those elements of direct relevance to this paper and refer the reader to Ref. Nielsen2010 for an introduction to quantum circuits and quantum computation. We will first introduce the quantum gates that can be implemented on the IBM quantum computers, and then decompose the two qubit unitary operations that appear in our trotterized evolution operator in terms of the elementary one and two qubit gates.

A quantum circuit consists of an array of quantum channels – which represent the physical qubits – and a series of quantum gates that are applied to them. These quantum gates are unitary operators that can be applied to one or more of the qubits. The IBM quantum devices can implement the CNOT gate along with an arbitrary single qubit gate, parametrized by three phases. The cobination of these gates forms a universal set that is, any $N$ qubit gate can be implemented using a combination of these gates, and in principle an arbitrary quantum computation can be performed.

There is a collection of single particle gates that are useful for writing quantum circuits. We write down a list of the most frequently used gates and how they are implemented on the IBM machines. Consider the computational basis to be the tensor product of single qubit states in the $z$-basis, i.e., $\{|\!\uparrow\rangle,|\!\downarrow\rangle\}$. All matrix forms of the gates are given in this basis and all measurements are made in this $z$-basis. It is important to note that gate multiplication reads left to right, whereas matrix multiplication is right to left, i.e.