Dynamics of disordered quantum systems with two- and three-dimensional tensor networks

TL;DR

This work shows how two- and three-dimensional tensor networks can accurately and efficiently simulate the quantum annealing dynamics of Ising spin glasses on a range of lattices and demonstrates that tensor networks are a viable approach for simulating large scale quantum dynamics in two and three dimensions on classical computers.

Abstract

Quantum spin glasses form a good testbed for studying the performance of various quantum annealing and optimization algorithms. In this work we show how two- and three-dimensional tensor networks can accurately and efficiently simulate the quantum annealing dynamics of Ising spin glasses on a range of lattices. Such dynamics were recently simulated using D-Wave's Advantage$2$ system [A. D. King et al, Science, 10.1126/science.ado6285 (2025)] and, following extensive comparisons to existing numerical methods, claimed to be beyond the reach of classical computation. Here we show that by evolving lattice-specific tensor networks with simple belief propagation to keep up with the entanglement generated during the time evolution and then extracting expectation values with more sophisticated variants of belief propagation, state-of-the-art accuracies can be reached with modest computational resources. We exploit the scalability of our simulations and simulate a system of over $300$ qubits, allowing us to verify the universal physics present and extract a value for the associated Kibble-Zurek exponent which agrees with recent values obtained in literature. Our results demonstrate that tensor networks are a viable approach for simulating large scale quantum dynamics in two and three dimensions on classical computers, and algorithmic advancements are expected to expand their applicability going forward.

Figures (8)

• Figure 1: Fully classical approach to simulating the quantum dynamics of a locally interacting system. In this work, we apply this to simulate the correlated, glassy dynamics induced by the Hamiltonian in Eq. (\ref{['Eq:Hamiltonian']}) on cylindrical, dimerized cubic, and diamond cubic lattices. The wavefunction is encoded in a tensor network whose structure matches the system's underlying geometry. This state, with maximum bond dimension $\chi_{\rm BP}$, is time evolved via a belief propagation-based simple update scheme whose efficiency is fundamental for capturing the large entanglement growth and keeping the wavefunction on the correct trajectory given bounded computational resources. Measurements of local and non-local observables are taken intermittently at designated measurement times via corrections of belief propagation such as loop corrections or, if the underlying lattice is planar or nearly planar, MPS message passing. If necessary, a truncation is performed via belief propagation prior to this measurement to enable the efficient use of these more controlled contraction schemes.
• Figure 2: a) Error $\epsilon_{c}$ --- see Eq. \ref{['Eq:ErrorMetric']} --- from two-dimensional tensor network simulations of a glassy quantum annealing protocol on an $8 \times 8$ cylindrical lattice. The same $N = 20$ disorder realization are used as in Ref. king2024. Error bars correspond to double the standard error on the mean. To obtain our results (orchid data points) we run a BP-based evolution protocol of the Trotterized circuit with a maximum bond dimension $\chi_{\rm BP} = 32$ and truncate down, with BP, to a final state of bond dimension $\chi$ before using cylindrical MPS message passing with a MPS rank of $r = 2\chi$ to calculate $\langle \sigma^{z}_{i}\sigma^{z}_{j} \rangle \ \forall i, j$. The blue dotted line is the average error from the D-Wave Advantage2 quantum annealer, where the shaded region represents $\pm 2\sigma / \sqrt{N}$ and the white markers represent the average error from the 2D TNS simulation in Ref. king2024. The inset shows the value for the spin glass order parameter $\langle q^{2} \rangle$ versus annealing time for a single disorder instance computed with BP-TNS. b) Schematic visualization of the cylindrical MPS message passing method to efficiently compute the two point correlators $\langle \sigma^{z}_{i} \sigma^{z}_{j} \rangle$ associated with $\ket{\psi}$. The operator $\sigma^{z}_{i}$ is inserted into $\langle \psi \vert \psi \rangle$ and boundary MPS is run around the cylinder a finite number of times until convergence. All two point correlators $\langle \sigma^{z}_{i} \sigma^{z}_{j} \rangle$ with $i$ fixed can then be computed in $\mathcal{O}(L)$ time, where $L$ is the total number of spins.
• Figure 3: Error $\epsilon_{c}$ --- see Eq. \ref{['Eq:ErrorMetric']} --- from two- and three-dimensional tensor network simulations of a glassy quantum annealing protocol for annealing times $t_{a} = 7{\rm ns}$ and $t_{a} = 20{\rm ns}$. We use the same $N = 20$ disorder realizations as those used in Ref. king2024. Error bars correspond to the standard error on the mean. The tensor network is time evolved with a simple BP-based evolution protocol with a maximum bond dimension $\chi_{\rm BP}$. Expectation values are then obtained from the TNS with either cylindrical message passing or BP loop corrections with configurations involving a maximum number $l_{\rm max}$ of antiprojectors (see Eq. (\ref{['eq:LoopCorrectedSumApprox']})). a)$R \times R$ cylindrical lattice geometry with $\chi_{\rm BP} = 32$ and the final state truncated, under the BP approximation, to $\chi = 10$. Message passing is performed with MPS rank $r = 2\chi$. b-c) Three-dimensional diamond cubic and dimerized cubic lattice geometries with $\chi_{\rm BP} = 16$ and $\chi_{\rm BP} = 6$ respectively. The final states are not truncated (i.e. $\chi = \chi_{\rm BP}$). The circled data point is annotated with the average clock time for simulating the dynamics ($T_{\rm Evo}$) and the average clock time for measuring a single two-point $z-z$ correlator ($T_{\rm Measure}$) for the corresponding system on a single Intel Icelake CPU.
• Figure 4: Verifying our own ground truth. Two point correlator in a $50$ qubit diamond cubic lattice after a quench with $t_{a} = 7$ns. The state is evolved with simple belief propagation with a maximum bond dimension of $\chi_{\rm BP} = \chi$ before being contracted with loop corrections up to size $l_{\rm max}$ to obtain $\langle \sigma^{z}_{4} \sigma^{z}_{5} \rangle$. The loop corrections above $l_{\rm max} = 0$ collapse onto each other on the scale of the plot with increasing bond dimension. Matrix product state simulations using the time evolving block decimation (TEBD) method Vidal2004 with swap gates included to apply non-local gates are shown in dotted lines and convergence is not possible to obtain despite us using over three days of walltime for $\chi = 1024$.
• Figure 5: Scaling collapse of the correlation function following the dynamics induced by Eq. (\ref{['Eq:Hamiltonian']}) on a cylindrical lattice. A two-dimensional tensor network approach is used, with a belief propagation-based time evolution protocol of the Trotterized circuit implemented with a maximum bond dimension $\chi_{\rm BP} = 32$. The final state is truncated down, with BP, to bond dimension $\chi = 8$ and MPS message passing with a MPS rank of $r = 16$ is used to calculate $\langle \sigma^{z}_{i}\sigma^{z}_{j} \rangle \ \forall i, j$ for all pairs of spins that are aligned along the same column of the cylinder. The correlation function is then extracted via Eq. (\ref{['Eq:CorrelationFunction']}) and plotted in the top panels, along with fits to the compressed exponential $C(d) \sim a_{0}\exp(-a_{1} d^{\alpha})$. Bottom panels: We rescale the distances as $\tilde{d} = t_{a}^{\frac{1}{\mu}}d$ to obtain collapse. The annotated value of $\mu$ is the one which provides a best fit to the squeezed exponential.