Efficient Tensor Network Simulation of IBM’s Eagle Kicked Ising Experiment

We report an accurate and efficient classical simulation of a kicked Ising quantum system on the heavy hexagon lattice. A simulation of this system was recently performed on a 127-qubit quantum processor using noise-mitigation techniques to enhance accuracy [Y. Kim et al., Nature, 618, 500–5 (2023)]. Here we show that, by adopting a tensor network approach that reflects the geometry of the lattice and is approximately contracted using belief propagation, we can perform a classical simulation that is significantly more accurate and precise than the results obtained from the quantum processor and many other classical methods. We quantify the treelike correlations of the wave function in order to explain the accuracy of our belief propagation-based approach. We also show how our method allows us to perform simulations of the system to long times in the thermodynamic limit, corresponding to a quantum computer with an infinite number of qubits. Our tensor network approach has broader applications for simulating the dynamics of quantum systems with treelike correlations.

Popular Summary

We perform a classical simulation of spin dynamics on a heavy hexagon lattice—a setup recently simulated on a quantum processor and used as evidence for “utility” in a pre-fault-tolerant era. Our method directly maps the geometry of the system to a tensor network and makes the assumption that the system is free of looplike correlations.

We demonstrate how our results are significantly more accurate than those achieved by the quantum processor and other classical methods. This success can be directly attributed to the surprising loop-free nature of the correlations in the system and allows us to perform dynamical simulations in the thermodynamic limit and to long times.

Our tensor network approach has broader applications for simulating the dynamics of quantum systems with treelike correlations and benchmarking quantum processor designs.

Figure 1

Left: structure of the Eagle quantum processor, which consists of a $6\times 3$ heavy hexagon lattice with two additional qubits ($113$ and $13$) added to the bottom left and right corners of the lattice. Right: tensor network structure used for our simulations of heavy hex lattices, with the network structure directly reflecting the lattice. On-site tensors ${\mathrm{\Gamma }}_v$ are colored in blue and possess physical, uncontracted, indices of dimension $2$ (represented by their dangling legs) and virtual indices of dimension $\chi$ (represented by the edges of the network), which are shared with neighboring tensors. Positive, diagonal bond tensors ${\mathrm{\Lambda }}_e$ live on the edges $e$ between the site tensors and are colored in gray.

Figure 2

Dynamics of the kicked Ising model on small heavy hex lattices of varying size. (a)–(c) Results from the BP-evolved TNS for several bond dimensions and systems of one to three heavy hexes are compared to exact state vector solutions at ${\theta }_h=\pi /4$. The results labeled with ${\chi }^{\ast }$ are obtained by computing expectations values of the BP-evolved TNS using exact contraction, while the results labeled with $\chi$ are obtained by computing expectation values of the BP-evolved TNS using BP contraction. In both cases, the states are evolved by applying gates using the BP approximation. Top plots show dynamics of the magnetization on the indicated (red ring) site, bottom plots show the BP error estimate [based on the spectrum of the edge environment—see Fig. 7 and Eq. (12)] for the $\chi =32$ TNS along the indicated edge $e$ versus the Trotter step. The dotted faded red line shows the relative error between the $\chi =32$ TNS magnetization approximated by BP and the exact magnetization while the faded green line shows the relative error between the $\chi =32$ TNS magnetization approximated by BP and the magnetization obtained by contracting the same TNS exactly. Insets in the bottom plots show the first $50$ singular values of the edge environment after $15$ Trotter steps. (d) BP error estimate approximated using a boundary MPS contraction scheme (see Appendix pp3 for details) of the TNS for $n\times m$ lattices of $n$ rows and $m$ columns of heavy hexes. (Top) TNS with $\chi =12$ and a boundary MPS contraction scheme with maximum MPS bond dimension $D=12$ at ${\theta }_h=\pi /4$. (Bottom) TNS with $\chi =8$ and a boundary MPS contraction scheme with maximum MPS bond dimension $D=8$ at ${\theta }_h=3\pi /8$.

Figure 3

Comparison for classically verifiable systems of our BP-approximated tensor network state approach to simulating the dynamics of the kicked transverse field Ising model on a heavy hex lattice versus the Eagle quantum processor and alternative tensor network methods. Expectation values with respect to the state $|\psi ({\theta }_h,5)⟩$ [i.e. following five Trotter steps of the dynamics of the model—see Eqs. (1) and (3)] are plotted, alongside exact results determined from light-cone simulations. (a) Average magnetization. (b) Weight-10 observable. (c) Weight-17 observable. The bottom plots show errors defined as the absolute difference between the simulation result and the exact result. For some data points the error from our TNS simulation is too small to fit on the scale of the plot and so these points are not marked. Circled, annotated points denote, for a given ${\theta }_h$, the memory required to store the state of the system at the given bond dimension $\chi$ and the walltime associated with performing the simulation and calculating the relevant observable on a Macbook M1 Pro.

Figure 4

Comparison for deeper circuits of our BP-approximated tensor network state approach to simulating the dynamics of the kicked transverse-field Ising model on the heavy hex lattice versus the Eagle quantum processor and alternative tensor network methods. Expectation values calculated following a number of Trotter steps of the dynamics of the model—see Eq. (1)—are plotted. (a) Weight-17 stabilizer after six steps of evolution. Here exact data is now available [18] and our BP data for $\chi =500$ is within ${10}^{-4}$ of the exact result for all ${\theta }_h$. (b) Weight-1 observable after $20$ steps. The shaded region shows the difference between our finite $\chi =500$ bond-dimension data and the data extrapolated to infinite bond dimension, where we believe the true answer lies. (c) Top and bottom plots show observables in (b) at ${\theta }_h=0.7$ and ${\theta }_h=1.0$, respectively, as a function of inverse bond dimension of the TNS. Red dashed lines represent a least-squares fit of the form $A+B/\chi$ taken on the data, and we take $A$ to be the predicted value of the observable in the limit $\chi \to \mathrm{\infty }$. Even in the limit $\chi \to \mathrm{\infty }$ there will generally be some deviation from the exact result due to the BP approximation that we use for evolving the state and computing expectation values (see Sec. 2). Our analysis of the errors due to BP for this system, however, suggest this deviation is likely to be very small. (d)–(f) Dynamics of $⟨Z_{62}⟩$ using the BP-approximated TNS approach versus a MPS approach (with bond dimension $2500$) with light-cone depth reduction (orange) for ${\theta }_h=0.6,0.8$, and $1.0$, respectively. Results from other methods at depth $20$ are shown as black circles (Eagle processor [10]), blue circles (truncated Pauli strings [20, 21]), purple diamonds (TNO [18]) and orange pentagons (MPO [19]). Inset shows average gate error from the MPS approach (pink circles) and absolute difference between the BP-approximated TNS and the MPS result (solid gray line).

Figure 5

Tensor network state for the infinite heavy hex lattice. The unit cell is a five-site tensor network. By adding in appropriate periodic boundary conditions and simulating the kicked Ising model with the BP-approximated TNS method we recover results for simulation of the infinite lattice with the BP-approximated TNS method. (Top right) Dynamics of $⟨Z_{62}⟩$ for ${\theta }_h=0.9$ simulated using the BP-approximated TNS method at bond dimension $\chi =400$. Crosses correspond to the dynamics of $⟨Z_3⟩$ on the periodic boundary condition (PBC) unit cell (where site $3$ is marked) with $\chi =400$ corresponding to an infinite heavy hex lattice while the dashed line corresponds to the dynamics of $⟨Z_{62}⟩$ on the $127$ qubit heavy hex lattice with $\chi =400$. (Bottom right) Dynamics of the bipartite entanglement entropy density $s$ [see Eq. (2) for the definition], calculated using the infinite BP-approximated TNS method and extrapolated to infinite bond dimension. Different curves correspond to different values of ${\theta }_h$, ranging in steps of $0.1$ from $0.1$ to $0.9$. The partition is shown by the red dotted line on the infinite lattice. The shaded area shows the difference between the extrapolated result (dotted line) and the finite $\chi =800$ result. The inset shows an example extrapolation in $1/\chi$ for ${\theta }_h=0.8$ after $40$ Trotter steps.

Figure 6

Steps to perform a “simple update” on the edge of a TNS in the Vidal gauge. The example pictured is two neighboring sites of degree $2$ (site tensor ${\mathrm{\Gamma }}_v$) and $3$ (site tensor ${\mathrm{\Gamma }}_w$), respectively. (i)–(ii) The bond tensors, site tensors, and gate are combined into a single composite tensor ${\mathrm{\Theta }}_{v,w}$. (ii)–(iii) A SVD is performed on ${\mathrm{\Theta }}_{v,w}$ and singular values of the resulting bond tensor ${\tilde{\mathrm{\Lambda }}}_{v,w}$ can be discarded in order to truncate the bond. (iii)–(iv) Resolutions of identity, using the original bond tensors, are inserted on the exposed edges. (iv)–(v) The inverse bond tensors are absorbed, resulting in the updated local tensors ${\tilde{\mathrm{\Gamma }}}_v$ and ${\tilde{\mathrm{\Gamma }}}_w$.

Figure 7

Forming the edge environment from the norm network of a tensor network state. One of the edges $e$ is split and all other indices of the network are contracted over, reducing the cut network to a single matrix where a singular value decomposition can be performed.

Figure 8

Comparison of various approaches to simulating the kicked transverse-field Ising model on the $127$-qubit Eagle processor geometry (Fig. 1). Approaches include a quantum processor (Eagle Processor [10]), MPS [10] and isoTNS [10] methods using Schrödinger evolution, a TNS (BP) approach using Schrödinger evolution and evolved with BP (this work), a Mixed TNS-TNO (BP) [21] approach using a combination of Schrödinger and Heisenberg evolution and evolved with BP, Sparse Pauli Dynamics [20, 21], a MPO [19] approach using Heisenberg evolution, a $31$-qubit full state simulation (31 Qubit Sim. [22]), and a TNO (SU) [18] approach using Heisenberg evolution and evolved with simple update, which is closely related to the BP approximation [14, 38], and contracted exactly to compute expectation values.

Figure 9

(Top) Two possible mappings from the qubits of the Eagle processor to the sites of a MPS ansatz. The left diagram shows the ordering used in Ref. [10] while the right diagram shows the ordering we use to obtain the results in Figs. 4–4, which we find leads to higher accuracy for a given MPS bond dimension. (Bottom) MPS ansatz with the site numbers corresponding to those in the above orderings.

Figure 10

Approximating the edge environment of an edge $e$ of a TNS on the heavy hex lattice of bond dimension $\chi$ using a boundary MPS-style scheme [50]. The contraction of the norm network is done by successively (top to bottom and bottom to top) approximating the contraction of pairs of rows of the heavy hex lattice as MPSs of bond dimension $D$. The resulting contraction is then reduced to the contraction of a pair of MPSs incident to the given row of the TNS.

Figure 11

Dynamics of the single-site magnetization of the kicked Ising model on heavy hex lattices of varying sizes. The TNS is evolved using the BP approximation and the magnetization is approximated with BP (solid lines) and boundary MPS (dashed lines) with the relevant parameters displayed above the plot. The middle plots show the relative error between calculating the expectation values with BP versus boundary MPS. Bottom plots show the BP error estimate based on Eq. (12) for an edge incident to the site where the magnetization is calculated.