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 (Nature volume 618, p. 500–505 (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 tree-like correlations of the wavefunction 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 tree-like correlations.

Introduction.Quantum computers are fundamentally noisy in nature and incur errors with each operation.To remedy noise and make the computer run as desired, one option could be quantum error correcting codes.In spite of significant advances in quantum technologies over the last decade, this type of error correction is currently out of reach.As such, huge efforts have been devoted to finding out whether current noisy quantum technologies could provide a practical advantage over classical computers without error correction.The emphasis here is on practicality, since expectation values of local observables in large classes of shallow quantum circuits can formally be computed in polynomial time on a classical computer [1,2] and so can the output distribution of a noisy random quantum circuit [3].The latter refutes the claim of violation of the extended Church-Turing thesis in Ref. [4], and makes it exceedingly unlikely any such violation can be achieved through noisy quantum computation without error correction.Noise can also have a similar effect as limiting entanglement, allowing noisy quantum computers to be simulated by tensor network methods [5][6][7].While such works put constraints on the power of noisy quantum devices, they also show that it may be difficult to classically simulate them in practice when the noise is sufficiently small.At the same time, error mitigation techniques have been proposed to further improve the accuracy of noisy quantum devices by classical post processing [8,9].With these techniques one can construct a better classical (biased) estimator for the desired noiseless result by paying the overhead of having to run several noisy quantum circuits.
The power of noise-mitigation techniques was recently demonstrated in quantum simulations of the dynamics of a two-dimensional (2D) transverse field Ising model on a 127-qubit heavy-hexagon lattice [10] (see Fig. 1).Expectation values of a number of observables were extracted after a few Floquet cycles using zero-noise extrapolation techniques [8].The experimental results were found to be much more accurate than several classical tensor network approaches applied to the same problem, even when the tensor network methods utilized significant computational resources.
In this work we adopt a tensor network approach that respects the qubit connectivity of the heavyhex lattice to simulate the dynamics of this kicked Ising model on systems of various sizes.Our tensor network ansatz is approximately contracted with belief propagation, which works best when correlations in the system remain "tree-like".By quantifying these correlations we explicitly show how this assumption becomes increasingly valid as the lattice size increases.This allows our method to achieve results, for the 127 qubit problem, to a much higher degree of accuracy than previously reported classical results or the error mitigated quantum computer.Even at large circuit depths where there are no exact results to benchmark against, we utilize extensive MPS calculations and boundary MPS approaches to provide substantial evidence that our results are still highly accurate and the correlations are tree-like.Finally, we conclude by showing how our method can be used to accurately simulate the thermodynamic limit of the problem to long-times, and therefore can simulate a high-depth quantum circuit involving an infinite number of qubits.Our work demonstrates the effectiveness of a belief propagation tensor network approach for solving many-body dynamics problems.We anticipate our chosen methodology will find success and serve as a benchmark when applied to problems with locally tree-like correlations and limited entanglement.
Model and Ansatz.Our focus here is on the dynamics of the Trotterized kicked transverse field Ising model given by the unitary where Z and X denote Pauli operators and ⟨v, v ′ ⟩ indicates that v and v ′ are neighbors on the corresponding lattice.The lattice we are concerned with is that of the 'heavy-hex' lattice which corresponds to a hexagonal lattice decorated with additional qubits along the edges (see Fig. 1).The dynamics of this model was recently simulated on the IBM Eagle quantum processor [10], which corresponds to a lattice of 6 × 3 heavy-hexagons plus two additional qubits.
In order to simulate this system on a classical computer, we adopt a tensor network approach that respects the qubit connectivity of the heavy-hex lattice (see Fig. 1).We fix a maximum amount of entanglement in the system by limiting the bond dimension χ of the network: each application of U (θ h ) will double the bond dimension of the network and thus it is necessary to limit the bond dimension and perform truncations.We then evolve our tensor network state (TNS) by application of the gates in U (θ h ) under the belief propagation (BP) approximation [11][12][13]; referring to the resulting TNS as a BP-approximated TNS.Unless otherwise stated, we also extract expectation values from the TNS using belief propagation.Explicit details of our BP-based method are provided in the Methods section and in Ref. [14].The BP method is fully controlled on trees but incurs a potentially small but uncontrolled approximation when there are loops in the network; the error from this approximation is, in general, smaller for larger loops (see Methods for further discussion on this).Here, we demonstrate that for a sufficiently large lattice, even at significant circuit depths, the correlations in this model remain 'tree-like' in the sense of the BP approximation giving very accurate results.Let us present these results.
Results.We start by considering lattices with a small number of heaxy-hexagons, where an exact state vector simulation is possible and our method can be directly benchmarked (see Fig. 2).Specifically, we compute the dynamics of the on-site magnetization for θ h = π/2 and lattices consisting of 1, 2, and 3 heavy-hexagons respectively.We also compute the separability of the 'edge environment' of the BP-approximated TNS along a chosen edge (see the BP error estimate section for a definition).This environment corresponds to the contraction of the TNS down to the given edge and when it is completely separable the belief propagation assumption is exact.The separability thus gives us an estimate of the error stemming from the BP approximation.
For short circuit depths n < 6 our method gives perfect agreement with the exact simulation and the edge environments are all completely separable because the light cone of circuit does not reach around the loops.For larger circuit depths there is deviation between the exact dynamics and the BPapproximated TNS dynamics.This can be explicitly characterized by a decrease in the separability of the edge environments, since the BP method approximates the environments used to perform the gate evolution as outer products of environments coming from incoming edges of a region of the system.Our small lattice simulations demonstrate something quite remarkable: as the number of heavyhexagons increases, the BP approximation at fixed χ improves significantly and the edge environment becomes highly separable even out to 20 Trotter steps.Moreover, in Fig. 2d) we use boundary MPS (see Boundary MPS for a description of the method) to show that this increased separability persists with increasing system sizes.
The accuracy of belief propagation does tend to improve with increasing system size (for instance, belief propagation on a finite ring directly gives the exact results of the thermodynamic limit, independent of system size) and may go some way towards explaining the increased accuracy of the BP simulation as we increase the number of heavy hexes.Nonetheless, the heavy-hex lattice still has finite loops in the thermodynamic limit and thus our results are indicative of something further: a level of interference between the loops in the lattice as we increase the system size.We note that Ising models have been observed to display some slow thermalization and confinement properties under a quantum quench [15,16] and we believe a similar effect is occurring here.
We now consider the 127 qubit heavy-hexagon lattice which corresponds to that of the IBM Eagle processor and we compare our method to the experimental quantum simulation results first from Fig. 3 of Ref. [10].Our results on smaller system sizes will be useful in explaining the accuray of our results here.In  The tensor network state (TNS) and gate evolution methods we use simulate the full 127-qubit system and result in highly accurate expectation values.In Fig. 3(a) we compute the expected value of the average single-site magnetization and show that we can obtain an accuracy1 of ∼ 10 −14 with a simulation that runs in less than 10 seconds on a laptop computer.Importantly, even for Figs.3(b) and Figs.3(c), where we consider higher weight observables which could be more strongly affected by loop correlations, we are still able to calculate these observables to a remarkable accuracy using our BP-approximated TNS.Specifically, we obtain values of these higher weight observables to orders of magnitude better accuracy than the quantum processor with a simulation that takes less than 4 minutes (for these observables) to run on a laptop and a state that takes up, at most, 0.3GB of memory.The remarkable accuracy we are able to achieve corroborates with our earlier analysis of the error of the BP approximation as a function of system size (see Fig. 2), where we found that the BP error decreases as we increase the system size for this model and lattice.
Turning next to larger numbers of Trotter evolution steps we show results in Fig. 4 for properties also computed by the quantum processor.For the n = 6 Trotter step simulation (Fig. 4a), where a weight-17 observable is measured, exact results are now available [17] and our BP-approximated TNS results at χ = 500 are within 10 −4 of these results for all values of θ h we plot.For n = 20 (Fig. 4b) exact data is currently unavailable and we push to sufficiently large bond dimensions to perform an extrapolation of the results from BP-approximated TNS to infinite bond dimension; demonstrating the reliability of a linear extrapolation in 1/χ for select θ h in Fig. 4c).We capture the known exact values at the Clifford points θ h = 0 and θ h = π/2.Notably, for this circuit in the region π 8 ≤ θ h ≤ 3π 8 there have been discrepancies between various classical methods and the quantum processor.In Appendix, we discuss the different classical simulation methods further and show a comparison in Fig. 8 of the various results [10,[17][18][19][20][21].
Beyond the results we presented above giving evidence for the separable nature of the edge environ-FIG.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 θ h = π/4.The results labelled with χ * are obtained by computing expectations values of the BP-evolved TNS using exact contraction, while the results labelled with χ are obtained by computing expectations 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 χ = 32 TNS along the indicated edge e versus the Trotter step.The dotted faded red line shows the relative error between the χ = 32 TNS magnetization approximated by BP and the exact magnetization while the faded green line shows the relative error between the χ = 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 the Appendix on Boundary MPS for details) of the TNS for n × m lattices of n rows and m columns of heavy hexes.Top) TNS with χ = 12 and a boundary MPS contraction scheme with maximum MPS bond dimension D = 12 at θ h = π/4.Bottom) TNS with χ = 8 and a boundary MPS contraction scheme with maximum MPS bond dimension D = 8 at θ h = 3π/8.ments in this system, we present further evidence that the BP-approximated TNS method is highly accurate throughout the whole phase diagram in Fig. 4d-f).Specifically, we compute the dynamics at every Trotter step of ⟨Z 62 ⟩ for several θ h and compare the BP-approximated TNS to our own MPS calculations as an independent check.Our MPS approach combines multiple non-trivial techniques: i) utilizing light cone depth reduction to calculate ⟨Z 62 ⟩ at every Trotter step, ii) using a higher bond dimension than that in Ref. [10], and iii) implementing an improved site ordering to lower the entanglement and gate error.We find that the difference between the BPapproximated TNS and MPS is directly correlated with the error from the MPS method and that both methods agree closely when the error in the MPS method is itself small.When the MPS error is small, one can consider MPS to be exact as it makes no uncontrolled approximations.The fact the BP-based method agrees with the MPS method when MPS exhibits very small errors suggest the BP error is also minimal.This is clearest for θ h = 0.6 in Fig. 4d) where we are able to push our MPS simulations to a bond dimension large enough for the average gate error to stay below 10 −4 .At all times for θ h = 0.6, there is clear agreement between BP-approximated TNS and MPS, yet disagreement at depth 20 versus the other methods shown [17][18][19][20][21].This agreement is only possible if the state possesses tree-like correlations and thus reinforces our earlier results on the general accuracy of the BP approximation for the dynamics of this system on this lattice.For larger θ h , in particular the θ h = 1.0 results shown in Fig. 4f, the BP results agree with the new MPS results until about step 10 where the MPS error starts to grow.We can be confident that the discrepancy there is due to the MPS method, not BP, because MPS is a controlled method that self-reports a significant error for larger steps and because the BP agrees with the result of ∼ 0 at step 20 which is predicted by a range of classical methods [17][18][19][20][21] due to the system being ergodic in this regime.
Dynamics of the infinite heavy-hex lattice -One of the powerful features of tensor network methods is their ability simulate infinite lattice models as long as they possess some form of translational invariance [22][23][24].Here we present results on the dynamics of the kicked transverse field Ising model on an infinite FIG. 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 |ψ(θ h , 5)⟩ (i.e following 5 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 θ h , the memory required to store the state of the system at the given bond dimension χ and the walltime associated with performing the simulation and calculating the relevant observable on a Macbook M1 Pro.
heavy-hexagon lattice, corresponding to a quantum computer with an infinite number of qubits.Again, we approximate the dynamics and take expectation values using the BP approximation.Given the evidence that we have presented on the accuracy of the BP approximation for the finite case of this system, especially for larger system sizes, we expect that these results are highly accurate.
For the infinite heavy-hex lattice there is a 5-site unit cell which can be tiled to produce the infinite lattice (see Fig. 5).It can be shown that simulating a periodic system on the unit cell using belief propagation corresponds to simulating the infinite lattice under the belief propagation approximation, where the sites of the periodic system constitute a unit cell of the infinite system [14].This is analogous to the standard approach to simulating infinite systems with the simple update tensor network method [25].Therefore we take the single unit cell, impose periodic boundary conditions, and present results for the belief propagation approximated dynamics.Fig. 5 illustrates this idea, showing the dynamics of the magnetization of the infinite system compared to the expected magnetization of a representative site near the center of the finite system.The extremely close agreement between the magnetization of the infinite and finite heavy-hex lattices provides strong evidence that boundary effects are minimal in the 127-qubit model and the results are already very close to those of the thermodynamic limit.
In Fig. 5 we also show the time-dependence of the bipartite entanglement entropy per edge s across a bipartition of the infinite lattice.Our infinite BPapproximated TNS method gives us an estimate of this quantity via the spectrum {λ i } of the bond tensors along the edges separating the partitions: which requires summing over the spectrum of only one of the bond tensors due to the fact that they are all identical along the cut being made.We observe that the entanglement in the system shows a sharp linear growth at short times before slowing down significantly and potentially saturating over a large time scale.This long-lived plateu is consistent with studies of the decay of the magnetisation in a smaller system [21] which suggests that the time-to-decay scales exponentially with the inverse of θ h .This long-time slow growth means that we can accurately simulate the infinite quantum processor up to large circuit Here exact data is now available [17] and our BP data for χ = 500 is within 10 −4 of the exact result for all θ h .b) Weight-1 observable after 20 steps.The shaded region shows the difference between our finite χ = 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 θ h = 0.7 and θ 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/χ taken on the data, and we take A to be the predicted value of the observable in the limit χ → ∞.Even in the limit χ → ∞ 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 the Methods section).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 ⟨Z62⟩ using the BP-approximated TNS approach versus a MPS approach (with bond dimension 2500) with light cone depth reduction (orange) for θ 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 [19,20]), purple diamonds (TNO [17]) and orange pentagons (MPO [18]).Inset shows average gate error from the MPS approach (pink circles) and absolute difference between the BP-approximated TNS and the MPS result (solid grey line).
depths for smaller values of θ h ≲ π/4.For larger values of θ h , the entanglement grows sufficiently quickly that, with our current resources, we are unable to accurately determine the entanglement entropy in the system beyond ∼ 25 Trotter steps.
State ansatz.Our ansatz for the wave function of the system is a tensor network which directly reflects the 'heavy-hex' qubit connectivity of the processor (see Fig. 1).The physical properties of a tensor network are invariant under a gauge symmetry corresponding to the insertion of any invertible matrix G e and its inverse G −1 e on any contracted, internal bond/edge e of the network.We take advantage of this symmetry to keep the tensor network in the Vidal gauge [25,26].This gauge corresponds to the choice of positive, diagonal bond tensors Λ e residing on the edges of the network and the on-site tensors Γ v of the network obeying certain isometric properties (see Eqs. ( 5) and ( 6)).These isometric properties are important for maintaining accuracy during the evolution of the network.We use |ψ(θ h , n)⟩ to denote the TNS of the system after n ≥ 1 applications of where |ψ(0)⟩ is the initial state of the system.We use the same initial state as in Ref. [10]: |ψ(0)⟩ = | ↑↑ . . .↑⟩.The single-site X rotations in U (θ h ) can be applied to |ψ(θ h , n)⟩ exactly and the two-site gates are applied approximately using the simple update FIG. 5. Tensor network state for the infinite heavy-hex lattice.The unit cell is a 5-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 ⟨Z62⟩ for θ h = 0.9 simulated using the BP-approximated TNS method at bond dimension χ = 400.Crosses correspond to the dynamics of ⟨Z3⟩ on the periodic boundary condition (PBC) unit cell (where site 3 is marked) with χ = 400 corresponding to an infinite heavy-hex lattice while the dashed line corresponds to the dynamics of ⟨Z62⟩ on the 127 qubit heavy-hex lattice with χ = 400.Bottom) 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 θ 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 χ = 800 result.The inset shows an example extrapolation in 1/χ for θ h = 0.8 after 40 Trotter steps.
[25] procedure (see the Methods section), which involves truncating the internal indices of the TNS to keep them less than or equal to a prescribed maximum bond dimension χ.Following a single Trotter step, the tensor network is regauged using belief propagation, a well established statistical inference algorithm which can be formulated for tensor networks [12] which we find improves the accuracy of the simple update procedure [14,27].Regauging before applying each gate with simple update would be most accurate and would be equivalent to performing each gate application with environments computed from the BP fixed point, however in practice we find that is too computationally expensive.Regauging after each Trotter step provides a good balance between accuracy and speed.Belief propagation allows us to rapidly determine the necessary transformation matrices by performing message passing over the net-work representing ⟨ψ(θ h , n)|ψ(θ h , n)⟩.The use of belief propagation as a method to efficiently 'regauge' a tensor network was recently formulated by some of the current authors in Ref. [14] and is closely related to other known gauging methods [27][28][29][30].We emphasize that the results presented here could have been achieved with those known gauging methods.
Measuring expectation values.We measure singlesite expectation values of the gauged state |ψ(θ h , n)⟩ using a rank-one approximation for the environments of local regions of the network (see Methods for more details).Such a method is only guaranteed to be exact in the limit of a tree network.Nonetheless, as demonstrated by our results here, the large loop structure of the heavy-hexagon lattice and dynamics of the model means that the network is locally treelike and therefore amenable to such an approximation.In Ref. [10] specific higher-weight observables were also measured to highlight that non-trivial results can be obtained in the "regime of strong entanglement" -such observables can be difficult to accurately measure for loopy tensor networks.We can, however, exploit the Clifford properties of the circuit at θ h = π 2 to transform the problem of measuring higher-weight observables after n Trotter steps into one of measuring a single-site observable after evolution by n Trotter steps with the propagator U (θ h ) and n Trotter steps with the propagator U ( π 2 ) † .
Specifically, the operator , (4) meaning this observable can be calculated, with respect to |ψ(θ h , 5)⟩, by measuring ⟨Z 13 ⟩ with respect to the state U 5 π 2 † |ψ(θ h , 5)⟩.We can reach this state straightforwardly by performing a 10 (5 + 5) Trotter step tensor network simulation.We use this 'extended time evolution' method to measure all higher-weight observables presented in Ref. [10] and note that this procedure is entirely generic: any long Pauli string can be generated out of a single Pauli by application of a Clifford circuit.Naturally this extended time evolution further increases the entanglement we must deal with as each application of U π 2 will double the bond dimension if no truncation is performed.Thus an upper bound on the bond dimension χ max needed to evolve the system by n Trotter steps and the measure a string which can be generated by n ′ applications of U π 2 to a single-site Pauli operator is 2 n+n ′ .Generically, longer strings will require larger values of n ′ .Figures 3 and 4 here, however, show that we can use values of χ ≪ χ max and still get very accurate results.
We would also like to emphasize that these higherweight observables could be measured with other tensor network methods.For instance, the planar nature of the TNS on the heavy-hex lattice means the boundary MPS [31] method can be directly applied to the norm of the state |ψ(θ h , n)⟩ to measure a desired higher-weight observable.We used the boundary MPS method in this paper to approximately compute the edge environment and also found extremely close agreement between BP and boundary MPS when measuring single-site observables (see Boundary MPS).
Conclusion.We have shown that a 127 qubit simulation of Floquet dynamics of the kicked Ising model on a heavy-hexagon lattice, recently simulated on a quantum processor in Ref. [10], can be performed accurately and with minimal computational resources with tensor networks.Our work demonstrates the importance of adopting a tensor network ansatz which reflects the spatial connectivity and entanglement structure of the system.The chosen gate evolution method for this ansatz is based on contracting the network under the belief propagation approximation.The computational scaling is O(Lχ 4 ) for a given Trotter step evolution, where L is the number of qubits and χ is the bond dimension.Like the commonly used and closely related simple update gate evolution method, this method directly assumes the lattice is locally tree-like, in other words assuming that loops have a minimal effect on the local properties of the network.We have presented evidence that this assumption becomes increasingly valid with increasing system size.This makes our gate evolution method highly accurate and reveals a striking loop-free behaviour to the dynamics of this kicked Ising model.We leveraged this understanding of the model and our method to perform accurate simulations of the long time dynamics of the model on an infinitely large heavy-hex lattice, corresponding to a quantum processor with an infinite number of qubits.Following an initial linear growth in entanglement, we find a remarkable result that the entanglement of the system appears to saturate over a large timescale.
Looking forward, we expect that the beliefpropagation based method employed here will allow efficient, accurate tensor network simulations of a range of dynamics problems and not just the one considered here.There are a number of geometriesespecially in dimensions higher than two where more neighbors typically causes more mean-field-like behavior [22]-where the effect of loops is small and we anticipate our method will find success.Moreover, while it is typically difficult to quantify the error stemming from the BP approximation, here we have performed a number of supplemental calculations (i.e.state vector calculations for small lattices, matrix product state simulations, and quantification of the BP error for increasing system sizes) to ascertain the accuracy of our results and we believe that such a methodology will be fruitful when considering other systems with a BP-based method.
Another takeaway from our results is that there can be many complementary routes toward classical simulation of quantum many-body systems, especially for those with physical structure such as separation of energy scales or locality.In addition to well-known properties such as low entanglement, shallow circuit depth, or low T-gate count (in the case of nearly Clifford circuits), our work highlights the importance of lattice topology as another key property that can be exploited.We also emphasize that tensor network methods are not limited to one-and two-dimensional systems, and higher-dimensional, non-planar systems can actually become more mean-field-like, allowing tensor network approaches to again work well [32,33].
Finally, we would like to comment that our work opens up new directions in which highly flexible and computationally inexpensive tensor network approaches can be used to benchmark new quantum processor designs and can better delineate which many-body quantum systems could become difficult for classical computing techniques.The software we used is part of the forthcoming ITensorNetworks.jlpackage [34] being developed at the Flatiron Institute (see the Computing resources and software packages section), which enables rapid testing and deployment of tensor network methods on arbitrary graphs.This software is continually being developed with the aim of tackling the type of simulation presented here.

The Vidal gauge
Our tensor network ansatz for the state of the system consists of both local tensors Γ v on the sites of the network and bond tensors Λ e which live on the edges of the network.The state is in the 'Vidal' gauge which is characterized by a set of isometric constraints on certain groupings of local tensors.Specifically, for every edge e = (v, v i ) connecting vertices v and a neighboring vertex v i , if we group the following tensors they form an isometric tensor obeying where the right-hand side denotes the identity matrix and we have used the set {v 1 , v 2 , . . .v d } to denote the d neighbors of vertex v. Having the TNS obey this property is important for maintaining accuracy when applying two-site gates and when taking expectation values.Below we detail the methods we use for performing these operations and then describe how we use belief propagation to maintain the Vidal gauge during our simulations.We emphasize that working in the Vidal gauge isn't strictly necessary, and one can perform all of the same operations (gate application, computing expectation values, etc.) and get the same results without transforming to the Vidal gauge by just using message tensors found from performing belief propagation on the TNS in an arbitrary gauge [14].This is why we opt to name our tensor network state as a 'BP-approximated TNS' and will discuss further in the sections below.
The simple update procedure To apply two-site gates to the tensor network state we adopt the simple update procedure [25].The procedure is depicted diagramatically in Fig. 6.For better efficiency, we use the 'reduced tensor' variant of simple update [35,36] (not shown in Fig. 6).The simple update procedure can be performed on a TNS in an arbitrary gauge by working with the fixed point message tensors found from belief propagation on the TNS and updating the message tensor on the edge where the gate is applied with the bond matrix returned from the SVD procedure [14].This is similar to previous work in Ref. [12] where message tensors were used for energy optimization.

Measuring expectation values within the Vidal gauge
In order to measure a single-site observable ⟨O v ⟩ on site v of our TNS we absorb the neighboring bond tensors onto the on-site tensor Γ v and contract the result with its conjugate, inserting the single-site observable O v along the physical index which is being contracted over.For the example of a site with 3 neighbors this can be visualized as where we have defined the square of the bond tensors Equation ( 7) approximates ⟨O v ⟩ by treating the environment as a tensor product of environments coming from each of the neighbors of v.If the tensors of the network obey the Vidal gauge conditions (Eqs. 5 and 6), this approximation is equivalent to computing the expectation value in an arbitrary gauge by using the fixed point message tensors of belief propagation as approximations of the environments [14,37].For networks that are locally tree-like, this can provide very good approximations for local observables.We find evidence that this approximation holds very well for the model and lattice studied in this work.In the section below we detail how we use belief propagation to maintain the Vidal gauge during our simulations.This helps maintain accuracy in the simple update procedure and when taking expectation values.

Belief propagation on a tensor network state
At the heart of our method for maintaining the gauge properties of our tensor network state is belief propagation (BP).BP is a well-established technique for approximating the marginals of the probability distributions of graphical models [38] and has recently gained interest in the context of contracting tensor networks [11][12][13]30].
To perform belief propagation on a TNS |ψ⟩ in the Vidal gauge we first absorb the square roots of the bond tensors onto the Γ v tensors.For instance, for the example of a tensor with 3 neighbors we define .(9) We then form the closed network which respresents the square norm ⟨ψ|ψ⟩ of |ψ⟩.This 'norm' network consists of the on-site tensors , where the summation is over the external indices s v of the tensors T v and (T v ) * (though in practice we keep the tensors separate for efficiency).
We next define a series of 'message tensors' over the the norm network, with M v,vi denoting the message tensor directed along the edge from T v to its neighbor T vi .The indices of M v,vi match the indices shared by tensors T v and T vi .Here, the direction of the edge is important and, generally, M v,vi ̸ = M vi,v .A set of self-consistent equations is defined for the message tensors: where the product runs over all d neighbors of v excluding v i and multiplication of two tensors implies a contraction over any common indices they share.This equation can be expressed diagrammatically as Initializing the message tensors, one can iterate these equations in an attempt to converge them.The converged messages then form a rank-one approximation of the exact environment for a given tensor T v , and can be used to approximate expectation values of the TNS.The procedure for taking an expectation value in Eq. ( 7) using a state in the Vidal gauge is equivalent to taking the same state in an arbitrary gauge and approximating the environment surrounding the site with the BP message tensors [14,37].The more tree-like the network, the better the approximation.Importantly, the BP message tensors can also be used to directly define a gauge transformation which, when applied to a TNS, brings it into the Vidal gauge and guarantees the satisfaction of Eqs. ( 5) and (6) along all edges of the network.A detailed description of this method, including extensive benchmarking and discussion of its relation to other gauging methods, is given in Ref. [14].
In our simulations here we perform belief propagation gauging after every Trotter step.This is to maintain accuracy in these procedures.We should note that for shorter depth circuits, such as those simulated in Fig. 3), gauging after every Trotter step does not have a significant affect on accuracy.The circuit here is not deep enough to significantly alter the isometric condition in Eq. ( 6) and affect the simple update procedure.For longer depth circuits, such as the one run in Fig. 4b), gauging every Trotter step is important and not performing gauging during the simulation leads to a significant loss in accuracy.
We also always gauge the state before taking expectation values (which, as we have discussed, is equivalent to computing the expectation values with fixed point BP message tensors).Even for the shorter circuits in Fig. 3, not performing gauging prior to taking an expectation value can noticeably affect the accuracy of the result -unless the bond dimension used is large enough for the simulation to be exact and so the simple update procedure preserves the gauge properties of the TNS.

Estimating the error of BP for general tensor networks
In order to quantify the error that belief propagation makes when approximating contractions of the network we compute the χ 2 × χ 2 'edge environment' associated with cutting the norm network along a given edge.The edge environment has previously been used in other contexts, such as improving the performance of periodic MPS methods [39] and fixing the gauge and performing truncations of general tensor networks [40].Here we use it to define a measure of the error of belief propagation.Specifically, we choose an edge of the norm network of our tensor network state and split the corresponding index (which is formed from the product of the bra and ket indices of that edge) -contracting the network down to a single matrix.This process is pictured in Fig. 7 for an example network.
Belief propagation essentially corresponds to assuming this edge environment is rank-1 along every edge of the network, i.e. that the spectrum of the positive diagonal matrix of singular values Λ is {λ 1 , λ 2 , ..., λ χ 2 } = {1, 0, ..., 0}.We can thus estimate the error from belief propagation along an edge e as where the quantity inside the square root is known as the 'index of separability' and is an established measure of the separability of a matrix based on a singular value decomposition [41].Although the singular values are dependent on the gauge of the tensor network, we fix to the symmetric gauge [14] when calculating them for consistency.This is the quantity we compute in Fig. 2 for the pictured networks and given edges.We also include examples of spectra which correspond to a variety of separabilities in the insets.The BP error measure defined in Eq. ( 12) is in the same spirit as the 'cycle entropy' defined in Ref. [40].In general for this model we find the choice of edge is unimportant and does not qualitatively change the results we observe.
It is worth making a further comment on the BP error in Eq. (12).In general, computing it is very costly due to the need to contract the full norm network, even if it is contracted approximately.One can observe, however, that generically the error will decrease exponentially in the smallest loop size of the lattice.This is most straightforwardly seen by considering a norm network which is a single translationally invariant periodic ring (loop) of size l.Assuming the spectrum of the on-site transfer matrix T is gapped then the BP error of the ring will scale as O (exp(−cl)) for some constant c related to the correlation length.This scaling should be generic for a lattice with smallest loop size l and gapped loop correlations.

Computing resources and software packages
The code used to produce the numerical results in this paper was written using the ITensorNetworks.jlpackage [34] -a general purpose and publicly available Julia [42] package for manipulating (gauging, contracting, partitioning, evolving, etc.) tensor network states of arbitrary geometry.It is built on top of ITensors.jl[43], which provides the basic tensor operations.Code is available in the current version of ITensorNetworks.jl for performing belief propagation, gauging, and the simple update procedure on arbitrary tensor network states.An example script is also included for specifically simulating the model in this paper with our BPapproximated TNS approach.All data was produced using a single node of Flatiron Institute's Forming the edge environment from the norm network of a tensor network state.One of the edges e is split and all other indices are network are contracted over, reducing the cut network to a single matrix where a singular value decomposition can be performed.
Rusty computing cluster.Timings and memory usage quoted for Fig. 3, however, are based on the same code being run on a Macbook M1 pro.The tensor network diagrams in this paper were produced using the publicly available package Graph-Tikz.jl[44], a general-purpose Julia package for visualizing graphs, including tensor networks.Data for the 127-qubit simulations is currently available at: https://github.com/JoeyT1994/BP-TNS-Data.

I. APPENDIX A: COMPARISON OF METHODS FOR CALCULATING Z62
Here we compare data for the expectation value of ⟨Z 62 ⟩ after 20 Trotter steps from a range of classical methods and the quantum processor.The results are shown in Fig. 8.We include results at χ = 200 and χ = 500 from our BP-evolved TNS approach, as well as results from an extrapolation in 1/χ to χ → ∞.The shaded region shows the difference between our finite χ = 500 bond dimension data and the data extrapolated to infinite bond dimension, where we believe the true answer lies.For smaller bond dimensions the BP-evolved TNS method overshoots the true results in the region π/4 ≲ θ h ≲ 7π/16, an arti-fact which is resolved by extrapolating in the bond dimension.Our extrapolated BP-evolved TNS results are in close agreement with results obtained in [20], which uses a combination of a BP-evolved TNS and a BP-evolved tensor network operator (TNO) to perform a mixed Schrödinger and Heisenberg evolution.Moreover, a method based on truncated Pauli strings [19,20], which approximates the Heisenberg evolution of the system and relies on a different type of approximation than tensor network methods and BP, also shows close agreement with our BP-evolved TNS approach and the aformentioned BP-evolved mixed TNS-TNO approach, although it slightly overshoots the value of ⟨Z 62 ⟩ for 5π 32 ≲ θ h ≲ 7π 32 -which can be observed in Fig. 4d) where we have numerically accurate results from MPS calculations.
Results in Fig. 8 are also shown for a matrix product operator (MPO)-based method of simulating the Heisenberg evolution of the system [18].This approach, however, undershoots the value in the region 4π  16 ≲ θ h ≲ 5π 16 -mapping the system to a 1D ansatz suffers from the drawback of requiring much larger bond dimensions than the more general tensor network-based approaches.The limitations of methods based on 1D ansatzes for this problem, like MPS and MPO, is especially clear when considering the MPS data originally reported in [10], along with the MPS results reported in this work in Fig. 4. We also include data from calculations based on a tensor network operator (TNO) approach [17], extrapolated to χ → ∞, which is evolved using simple update (SU) -which is closely related to BP [14,37] -and then is contracted exactly to compute expectation values 2 .This method overshoots the true value of ⟨Z 62 ⟩ in the region π 8 ≲ θ h ≲ 2π 8 -which can be directly observed FIG. 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) [20] approach using a combination of Schrödinger and Heisenberg evolution and evolved with BP, Sparse Pauli Dynamics [19,20], a MPO [18] approach using Heisenberg evolution, a 31-qubit full state simulation (31 Qubit Sim.[21]), and a TNO (SU) [17] approach using Heisenberg evolution and evolved with simple update, which is closely related to the BP approximation [14,37], and contracted exactly to compute expectation values.
We also would like to point out work in Ref. [45], which used a dissipative mean-field theory to simulate the 127-qubit system out to short circuit depths (∼ 5 Trotter steps) and obtained relatively accurate results -although with some noticeable deviations from the true result, which can be computed to an accuracy ∼ 10 −14 with our BP-evolved TNS method.

APPENDIX B: MPS CALCULATIONS
We performed MPS calculations for benchmarking in the unverifiable regime in Fig. 4d-f).We implemented an ordering of the sites of the heavy-hex lattice which is distinct from that in Ref. [10].Specifically, Fig. 9 shows the qubits of the Eagle processor numbered as 1 through 127 which then map directly to the sites 1 through 127 of the MPS we use.Two numberings are shown: our own and the original ordering from Ref. [10].We find that our own ordering results in slightly lower truncation errors when simulating the 127-qubit kicked Ising model at fixed MPS bond dimension χ.We break the two-site term in the propagator U (θ h ) down into a product of, at most, 5 commuting matrix product operators (MPOs) of bond dimension 2. The MPS is then evolved by successively applying these MPOs and after each appli-cation truncation is performed of the state down to a maximum bond dimension χ.The single-site gates are applied at the start of each Trotter step exactly.If implementing all two-site gates our decomposition of U (θ h ) involves 5 MPOs.In practice, however, we employ light cone depth reduction (LCDR), which means that if there are n ′ Trotter steps left until we take a measurement on a specific site, we only include two-site gates in the MPOs which are within the remaining light cone of the simulation.This means that for a given θ h in Fig. 4 we perform 20 separate simulations to get the most accurate MPS results for the value of ⟨Z 62 ⟩ after each number of steps n.
In order to approximate the error of this simulation we calculate the sum of the singular values discarded during the application of a single MPO, which we call ϵ i , where i refers to the MPO being applied.The total approximate error during the n trotter step simulation (and plotted in the insets of Fig. 4d-f) is then [5] where the sum runs over all N MPO-MPS applications (which may be less than 5n due to LCDR) during the simulation up to n Trotter steps.The error per gate that we calculate is approximately the same as applying the gate exactly and taking the overlap with the truncated state.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 Fig. 4d-f), 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.

APPENDIX C: BOUNDARY MPS FOR THE HEAVY-HEX LATTICE
In order to approximate the error of BP using the separability of the edge environment (see the BP error estimate section) for lattices that are too large to contract exactly, for example the larger lattices in Fig. 2d), we employ a method similar to boundary MPS [31,46,47] but generalized to arbitrary network structures [48] to approximately contract the norm network down to a chosen edge to compute the edge environment.Specifically, we approximate the contraction of two rows of the heavy hex as an MPS of fixed bond dimension D -with the external indices of the given region of the TNS being mapped to the dangling indices of the MPS.This allows us to approximate the successive contraction of rows (from top to bottom and bottom to top) of the heavy-hex TNS as a sequence of MPS-MPO contractions.This process is depicted diagramatically in Fig. 10.
We also computed the difference between calculating the magnetization of the BP-evolved TNS using the belief propagation method versus using boundary MPS for the parameters considered in Fig. 2d).We present the results here in Fig. 11 and observe how close the values obtained are -with the largest lattices showing differences between the two methods which are on the order of 10 −4 .The difference between the two methods generally decreases as the system size is increased.FIG. 10.Approximating the edge environment of an edge e of a TNS on the heavy-hex lattice of bond dimension χ using a boundary MPS-style scheme [48].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.FIG.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.

FIG. 1 .
FIG. 1. Left: Structure of the Eagle quantum processor which consists of a 6 × 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.Onsite tensors Γv are coloured in blue and possess physical, uncontracted, indices of dimension 2 (represented by their dangling legs) and virtual indices of dimension χ (represented by the edges of the network) which are shared with neighboring tensors.Positive, diagonal bond tensors Λe live on the edges e between the site tensors and are coloured in grey.

Fig. 3
Fig.3we overlay the quantum simulation results with that of our BP-approximated TNS dynamics shown as cross symbols.Here expectation values are measured after 5 Trotter steps and exact results, based on brute force light cone simulation techniques, are available to allow us to directly assess the errors.The tensor network state (TNS) and gate evolution methods we use simulate the full 127-qubit system and result in highly accurate expectation values.In Fig.3(a) we compute the expected value of the average single-site magnetization and show that we can obtain an accuracy 1 of ∼ 10 −14 with a simulation that runs in less than 10 seconds on a laptop computer.Importantly, even for Figs.3(b) and Figs.3(c), where we consider higher weight observables which could be more strongly affected by loop correlations, we are still able to calculate these observables to a remarkable accuracy using our BP-approximated TNS.Specifically, we obtain values of these higher weight observables to orders of magnitude better accuracy than the quantum processor with a simulation that takes less than 4 minutes (for these observables) to run on a laptop and a state that takes up, at most,

FIG. 4 .
FIG. 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 6 steps of evolution.Here exact data is now available[17] and our BP data for χ = 500 is within 10 −4 of the exact result for all θ h .b) Weight-1 observable after 20 steps.The shaded region shows the difference between our finite χ = 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 θ h = 0.7 and θ 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/χ taken on the data, and we take A to be the predicted value of the observable in the limit χ → ∞.Even in the limit χ → ∞ 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 the Methods section).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 ⟨Z62⟩ using the BP-approximated TNS approach versus a MPS approach (with bond dimension 2500) with light cone depth reduction (orange) for θ 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[19,20]), purple diamonds (TNO[17]) and orange pentagons (MPO[18]).Inset shows average gate error from the MPS approach (pink circles) and absolute difference between the BP-approximated TNS and the MPS result (solid grey line).
always a single Pauli string and its expectation with respect to |ψ(θ h , n)⟩ can be obtained by further evolving |ψ(θ h , n)⟩ with U ( π 2 ) † and then measuring Z i .For instance, X 13,29,31 Y 9,30 Z 8,12,17,28,32 = U 5 π 2 Z 13 U 5 π 2 † FIG.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 Γv) and 3 (site tensor Γw) respectively.i) -ii) The bond tensors, site tensors, and gate are combined into a single composite tensor Θv,w.ii) -iii) A SVD is performed on Θv,w and singular values of the resulting bond tensor Λ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 Γv and Γw.