Efficient tensor network simulation of IBM’s largest quantum processors

We show how quantum-inspired 2d tensor networks can be used to efficiently and accurately simulate the largest quantum processors from IBM, namely Eagle (127 qubits), Osprey (433 qubits) and Condor (1121 qubits). We simulate the dynamics of a complex quantum many-body system— specifically, the kicked Ising experiment considered recently by IBM in Nature 618, p. 500–505 (2023)—using graph-based Projected Entangled Pair States (gPEPS), which was proposed by some of us in PRB 99, 195105 (2019). Our results show that simple tensor updates are already sufficient to achieve very large unprecedented accuracy with remarkably low computational resources for this model. Apart from simulating the original experiment for 127 qubits, we also extend our results to 433 and 1121 qubits, and for evolution times around 8 times longer, thus setting a benchmark for the newest IBM quantum machines. We also report accurate simulations for infinitely-many qubits. Our results show that gPEPS are a natural tool to efficiently simulate quantum computers with an underlying lattice-based qubit connectivity, such as all quantum processors based on superconducting qubits.


I. INTRODUCTION
We are currently witnessing an unprecedented technology race to develop practical large-scale quantum computers.While several hardware architectures have been developed, the largest available quantum processors are those built with superconducting qubit technology [1].In this setting, IBM's quantum roadmap is particularly promising, with the delivery of increasingly-larger quantum processors every year: Eagle with 127 qubits in 2021, Osprey with 433 qubits in 2022, and Condor with 1121 qubits and expected by the end of 2023.These are presently among the most powerful quantum machines worldwide.Furthermore, a large effort is being devoted to mitigate errors in the processors, to become able to run longer quantum circuits and therefore increase quantum volume.Such error mitigation was pushed to an unprecedented level in a recent paper [2], where the IBM team simulated the dynamics of a kicked quantum Ising model on a 127-qubit 2d lattice that matched the connectivity topology of Eagle's quantum computer.These results are a great step forward towards practical quantum computation in superconducting quantum processors.However, and unlike originally thought, they are still far from any sort of quantum advantage: as pointed out by several authors, the experiment can be simulated efficiently by purely classical means [3][4][5][6][7][8][9], and specially by methods using quantum-inspired tensor networks [10,11].
In this paper we go one step further, and show how 2d tensor networks based on Projected Entangled Pair States (PEPS) [12][13][14] can be used to simulate IBM's largest quantum processors: Eagle (127 qubits), Osprey (433 qubits) and Condor (1121 qubits).We show this by simulating the kicked Ising experiment mentioned above, with unprecedented accuracy and not just for 127 qubits, as in the original proposal, but for the larger quantum processors and longer evolution times, setting new bench-marks for those machines.We use graph-based Projected Entangled Pair States (gPEPS) [15], a type of 2d tensor network algorithm which provides great flexibility in adapting to new lattices, both of finite and infinite size.We conclude that gPEPS is a natural tool to efficiently and accurately simulate slightly-entangled quantum computations on quantum computers that have an underlying lattice-based qubit connectivity.

II. MODEL
We implement a simulation of IBM kicked quantum Ising model.Specifically, we consider the dynamics generated by the spin-1/2 Hamiltonian with Z i , X i being the Z and X Pauli matrices at site i, coupling J, transverse magnetic field h, and where the the sum of interactions is over nearest neighbors ⟨i, j⟩ on a lattice matching the topology of IBM's quantum processors.A first-order trotterization of the time evolution leads to the unitary operator with θ h a parameter controlling the relative strength of the field with respect to the spin-spin interaction.Starting from an initial state with all spins in the |0⟩ state (i.e., all "up"), we simulate the dynamics by applying the unitary operator U (θ h ) multiple times, therefore generating the state arXiv:2309.15642v3 [quant-ph] 2 Apr 2024 after n applications of the operator, m being the number of spins in the lattice.

III. METHOD
Here we simulate the dynamics of the above model on finite heavy-hexagon lattices with open boundary conditions and 127, 433 and 1121 vertices, respectively matching the connectivity of qubits in Eagle, Osprey and Condor, see Fig. 1.For this, we adapt the gPEPS method [15]-initially proposed for infinite systems-to finite-size lattices.We also use gPEPS to study the heavy-hexagon lattice in the thermodynamic limit with a unit cell of 10 sites.The gPEPS algorithm is a quantum-inspired tensor network method that allows to easily simulate systems on generic lattices with desired dimensionality.As such, it is a natural evolution of the original iPEPS algorithm [13,14] to simulate two-dimensional quantum lattice systems and also proposed years ago by some of us.gPEPS makes use of the simple update of tensors [16][17][18] and a mean-field approximation for expectation values.These approximations are accurate for slightlyentangled 2d quantum lattice systems.Though they can be systematically improved, e.g., by using full and fastfull [19] updates and corner transfer matrices [14], we have observed that the simplest of our approximations is already capable of simulating the system at hand with large accuracy.In our approach, the bond dimension of the PEPS tensor network is χ, and is also the truncation parameter in our simulations: larger the χ, the larger the allowed entanglement per bond.For comparison, we have also studied the effect of re-gauging the PEPS using belief propagation (BP) after each trotter step, as proposed in Ref. [3].

A. Benchmarking
First, we simulate the 127-site heavy-hexagon lattice from Fig. 1(a).Using gPEPS we perform the unitary evolution (U (θ h )) n up to n = 5 trotter steps, followed by computing expectation values using a mean-field approximation.We reproduce Fig. 3 from Ref. [2] in our Fig. 2, where we compare the outcome of our simulations with those obtained from experimental calculations performed on the IBM Eagle quantum processor.Additionally, we benchmark our findings against other tensor network methods.Comparison of our average magnetization values with the available light cone-based exact solution [2] shows exceptional precision (∼ 10 −15 of absolute error), with each data point taking on average 2 seconds to run on a standard desktop PC (Windows 11, Intel i7-11700 @2.50GHz, 16 GB RAM).Our results not only surpass the outcomes of IBM's quantum simulations, but they also outperform some of the best state-of-the-art tensor networks methods in both precision and speed.
Additionally, to study the effect of Belief Propagation (BP) gauging we have independently simulated the unitary evolution of 5 trotter steps, where we do BP gauging after each trotter step.We find that BP does not improve accuracy [20], as can be seen in Fig. 2(a), even though the average computational time per point increased to 9.2 seconds.
We also computed the expectation value of "higherweight" observables (Appendix.B) reported in the IBM experiment, as shown in Figs.2(b,c).Here, we have also included the tensor network results from Ref. [3] for comparison.In these plots we provide the expectation values for the Weight-10 and Weight-17 operators, acting respectively on 10 and 17 lattice sites, across a range of θ h values.The plots show that we obtain better precision than the quantum processor with a small bond dimension χ = 32, requiring an average compute time of 10 seconds per computed data point.As expected, we see an increase in accuracy by ramping up the bond dimension to χ = 64 and χ = 128 in Fig. 2(c).
However, in the range θ h ∈ (π/8, π/4) in Fig. 2(b), the Weight-10 observable expectation value with χ = 32 is found to be more accurate than the higher bond dimension results, the explanation can be found in the Appendix.A.
Next, we have studied the case in which the unitary evolution spans more than 5 trotter steps.This involves simulating the state corresponding to the extended-depth quantum circuit, as shown in Figure 4 of Ref. [2].We computed the Weight-17 observable after 6 trotter steps and compare it with the result obtained by the 127qubit quantum processor, and our results can be found on  Fig. 3(a).Comparison against the exact result [2] shows that gPEPS simulations with bond dimension χ = 64 already outperform the Eagle quantum processor in accuracy.As expected, we have observed further enhancements in accuracy for larger bond dimension χ = 128.
To further test the algorithm we performed simulations involving longer unitary evolutions, with 20 trotter steps, and across a range of θ h values.We computed the expectation value of the Weight-1 (single-site) operator, as shown in Fig. 3(b).Notably, we achieved numerically ex-act results for the Clifford points θ h = 0 and θ h = π/2 when using bond dimension χ = 64.While an exact solution for longer unitary evolution remains elusive, we were able to compare with tensor network results in the infinite bond dimension limit [3], obtained from finiteentanglement scaling.While our χ = 64 bond dimension accurately captures points for θ h ≲ 3π/16, noticeable deviations become apparent beyond this regime.As shown in the figure, we could improve the accuracy significantly by increasing the bond dimension to χ = 128, 256 and  512.Finally, to show the reliability of our method we have studied the finite-entanglement scaling of ⟨Z 62 ⟩(θ h ) with the inverse bond dimension (1/χ), as shown in Fig. 3(c) for θ h = 1.0 and θ h = 0.7.We have plotted the values for 9 different bond dimensions and the fitting of extrapolation is done for the 5 highest bond dimensions.While for θ h = 0.7 we start to see a tendency to saturation at large bond dimension, at θ h = 1.0 we see no clear evidence yet of bond dimension saturation, which indicates that the quantum state has large entanglement.

B. Large systems
Next, and thanks to the computational efficiency of gPEPS, we have successfully simulated larger IBM quantum systems involving 433 (Osprey) and 1121 (Condor) qubits, corresponding to the heavy-hexagon lattices in Figs.1(b,c).While the original IBM experiment in Ref. [2] was implemented only for 127 qubit system, we understand our results as a benchmark for future experiments on these larger quantum processors.In addition, we have also implemented gPEPS for the system with infinitelymany qubits, by assuming translation invariance in the PEPS tensor network with a unit cell of 10 sites.For all these sizes, we have simulated the unitary evolution of 5 trotter steps and computed a number of observables.Our results can be found in Fig. 4. First, we computed the average magnetization in the Z-direction for χ = 32, as shown in Fig. 4(a).We can see that there are minimal differences for all sizes (127, 433, 1121 and infinite), indicating that the bulk of the system is already quite close to the thermodynamic limit already for 127 qubits.For this plot, the average simulation time of one data point for sizes 127, 433, 1121 and infinite are respectively 2, 8.3, 50, and 0.17 seconds.Next, in order to test possible boundary effects, we compare in Fig. 4(b) the results for a Weight-10 observable near the open boundary of the heavy-hexagon lattice, for the three finite-size systems, again for χ = 32.The composition of the observable and its calculation is discussed in the Appendix.B.As we can observe, there is no appreciable difference in the result, signalling again that even the smallest lattice is already close to the thermodynamic limit.Last but not least, we have also computed the expectation value of a Weight-17 operator deep inside the bulk of the system, for all lattices (including the infinite one for χ = 64), and the results are in Fig. 4(c).Minimal difference among the results of different system sizes show that 127 qubit system is already very close to the thermodynamic limit.

C. Longe time evolution
All the above results motivate us to test the limit of our simulation method.One should expect that longtime evolutions may create a large amount of entanglement that is hard to be captured by the gPEPS technique.This could be captured as a loss of convergence with the bond dimension χ in our simulations, setting then a benchmark: a quantum computer claiming quantum advantage in simulating this model should (at the very least) be able to compute time evolutions longer than those for which gPEPS loses convergence.Therefore, to understand the limit of our method we computed the results in Figs.(5) and (6), respectively for the time evolution of the magnetization of a site deep in the bulk and the average magnetization over all sites, for the three considered sizes (127, 433 and 1121 qubits).The results are for a large number of Trotter steps (between 37 and 39), and for the largest bond dimension that we could achieve for each size.χ = 270 is the maximum common bond dimension that we could simulate for all three system sizes.Bond dimension scaling is shown in Figs.7.As we can see in the plots, for sizes we find the same result: the study of consecutive bond dimensions shows that the observables converge even for a large number of Trotter steps.This is indeed surprising and seems to indicate that the gPEPS technique is particularly suited to capture the entanglement structure of the heavy hexagon lattice.The reason behind this may be that this lattice, after all, can be quite well approximated by a tree-like structure with no loops (Appendix.C).And the properties of such loop-free structures can be captured with large precision by the simple tensor update that we use in gPEPS.

V. CONCLUSIONS
In this paper we have simulated IBM's kicked Ising experiment [2] on heavy-hexagon lattices with 127, 433, 1121, and infinitely-many qubits, using a quantum-inspired tensor network technique tailored to higherdimensional systems.Our method uses the so-called gPEPS algorithm, which is remarkably efficient and accurate.Our method not only reproduces the results of the original experiment for 127 qubits but also settles new benchmarks for large quantum circuits on IBM's Eagle, Osprey, and Condor quantum processors.We conclude that gPEPS is a natural tool to efficiently and accurately simulate slightly entangled quantum computations on quantum computers with an underlying lattice-based qubit connectivity, be it in 2 or higher dimensions, going much beyond the capabilities of other tensor network structures.In particular, it is an ideal tool to classically simulate quantum computers based on superconducting qubits.A relevant question triggered by our results is whether quantum processors based on artificial qubits (e.g., superconducting qubits, quantum dots, etc.), and with an underlying lattice-based connectivity, can reach a sufficiently-low noise level so that they cannot be simulated classically by some tailored tensor network algorithm.It would also be interesting to assess gPEPS in the simulation of other types of quantum hardware with allto-all qubit connectivity, such as trapped ions and neutral atoms.shorter-range correlations, owing to the absence of loops.In the Figs.8, we present a comparison between a patch of 2D heavy-hexagon lattice and a 2D square lattice, illustrating the environment of the blue edge (pair of spins).It is evident that for a 2D square lattice, the 'n-th'neighbor environment contains loops for n > 1. Conversely, in the case of a heavy-hexagonal lattice, the 'nth'-neighbor environment exhibits loops for n > 5. Thus, for models with short-range correlations, such as the Ising transverse field away from the critical point, the heavyhexagon lattice behaves akin to a tree structure locally.Consequently, local tensor updates employed in gPEPS, as well as local measurements, excel in capturing the realtime dynamics of the kicked Ising experiment.

FIG. 1 :
FIG. 1: (Color online) Different heavy-hexagon lattices, corresponding to the topology of qubit connectivity of three IBM quantum processors: (a) Eagle, with 127 qubits; (b) Osprey, with 433 qubits; (c) Condor, with 1121 qubits.Every dot in the lattices corresponds to a superconducting qubit, and every link corresponds to a qubit-qubit coupling.

10 Z
FIG.3:(Color online) Comparison of the gPEPS simulation with higher number of trotter steps with the Eagle quantum processor and various other tensor network methods.(a) Weight-17 observable computed after 6 trotter steps with respect to the state |ψ(θ h , 6)⟩.The bottom plot shows the absolute difference between our simulation and the available exact result.(b) Weight-1 expectation value computed after 20 trotter steps with respect to the state |ψ(θ h , 20)⟩.Because of the absence of exact result for this simulation, we have computed the absolute difference between our simulation and the BP-approximation tensor network state approach with χ → ∞, presented in the bottom subplot.(c) Finite-entanglement scaling of Weight-1 observable expectation value ⟨ψ(θ h , 20)|Z 62 |ψ(θ h , 20)⟩ with respect to the inverse of bond dimension (1/χ) for two distinct θ h values.Labeling of qubits is done sequentially, from left to right and top to bottom, starting with 0.

10 −6 10 0
FIG.4: (Color online) Results of simulating various IBM quantum chips with higher number of qubits using gPEPS: Eagle processor with 127 qubits, Osprey with 433 qubits, Condor with 1121 qubits, and the heavy-hexagon lattice in thermodynamic limit.(a) Average magnetization, (b) Weight-10 observable near the open edge, (c) Weight-17 observable deep inside the bulk.The structure of the Weight-10 and Weight-17 observables is discussed in (Appendix.B).
FIG. 5: (Color online) Long time evolution of the magnetization for a site in the bulk at θ h = 1.0 and the three different sizes: (a) 127 qubits, up to χ = 560 and 39 Trotter steps, (b) 433 qubits, up to χ = 370 and 38 Trotter steps, (c) 1121 qubits, up to χ = 270 and 37 Trotter steps.Lower panel shows relative errors with respect to the maximum achievable bond dimension.

OOO
FIG. 6: (Color online) Long time evolution of the average magnetization of the whole system at θ h = 1.0 and the three different sizes: (a) 127 qubits, up to χ = 560 and 39 Trotter steps, (b) 433 qubits, up to χ = 370 and 38 Trotter steps, (c) 1121 qubits, up to χ = 270 and 37 Trotter steps.The lower panel shows relative errors with respect to the maximum achievable bond dimension.

FIG. 7 :
FIG.7:(Color online) Bond dimension scaling of the average magnetization for different system sizes.
Color online) Comparing the gPEPS method in simulating the kicked transverse field Ising model against the 127-qubit IBM Eagle quantum processor and various other tensor network methods.The operator expectation values shown in (a) Average Magnetization, (b) Weight-10 observable, and (c) Weight-17 observable, are computed with respect to the state |ψ(θ h , 5)⟩.Each bottom plot shows the absolute difference between the light-cone based exact results and the results obtained through simulations (gPEPS and Eagle processor).Labelling of qubits is done sequentially, from left to right and top to bottom, starting with 0.