A Quantum-Inspired Algorithm for Wave Simulation Using Tensor Networks

TL;DR

This work reformulates the isotropic acoustic wave equation as a first-order Schrödinger-like problem using a massless Dirac analogy, then develops a quantum-inspired classical algorithm based on Tensor Networks to simulate wave propagation. By diagonalizing spatial derivatives with the Quantum Fourier Transform and encoding states and operators as MPS/MPO, the method enables exact in Fourier space propagation and scalable approximate propagation via Trotterization, achieving efficiency gains on large grids. The authors demonstrate 2D results with a Ricker wavelet on grids up to about $10^{12}$ points and show promising 3D scaling, highlighting potential speedups over FFT-based approaches and outlining extensions to Maxwell’s equations. Overall, the paper advances quantum-inspired classical PDE simulation by combining Dirac-based reformulations with MPO/QFT techniques, offering a path toward high-resolution, tensor-network-based wave simulations and potential near-term computational advantages.

Abstract

We present an efficient classical algorithm based on the construction of a unitary quantum circuit for simulating the Isotropic Wave Equation (IWE) in one, two, or three dimensions. Using an analogy with the massless Dirac equation, second order time and space derivatives in the IWE are reduced to first order, resulting in a Schrödinger equation of motion. Exact diagonalization of the unitary circuit in combination with Tensor Networks allows simulation of the wave equation with a resolution of $10^{13}$ grid points on a laptop. A method for encoding arbitrary analytical functions into diagonal Matrix Product Operators is employed to prepare and evolve a Matrix Product State (MPS) encoding the solution. Since the method relies on the Quantum Fourier Transform, which has been shown to generate small entanglement when applied to arbitrary MPSs, simulating the evolution of initial conditions with sufficiently low bond dimensions to high accuracy becomes highly efficient, up to the cost of Trotterized propagation and sampling of the wavefunction. We conclude by discussing possible extensions of the approach for carrying out Tensor Network simulations of other partial differential equations such as Maxwell's equations.

Figures (6)

• Figure 1: Diagram representing either a quantum circuit or an MPO for Eq. \ref{['eq: Trotter MPS Evol']} in three spatial directions with $n_x,n_y,$ and $n_z$ qubits each. Here $\mathcal{D}_l = e^{-ic_l\Delta t\hat{\sigma}^l\hat{D}_l}$ and $Q_{l}$ is the Quantum Fourier Transform in the $l^{\text{th}}$ spatial direction. In two dimensions, one can simply drop the $Z$ direction and merge the two $\mathcal{D}_y$ operators. For $X$ and $Y$$c_l = 1/2$ and for $Z$, $c_l = 1$. In one dimension for the small angle approximation this reduces to the same approach by Wright et. al Wright2024, see the text for details.
• Figure 2: The two-dimensional Ricker Wavelett centered at $(0.5,0.5)$ with $\sigma=0.1$ represented across ${50}\times {50}$ qubit sites, or equivalently about $10^{30}$ grid points. This histogram is constructed from $10^7$ samples of the MPS state, grouped into $512\times 512$ bins. Preparing this state using ITensorMPS.jl "out of the box" with no attempt at optimization takes about 43 seconds on a laptop.
• Figure 3: The two dimensional Ricker wavelett for $\mu=0.5$ and $\sigma=0.1$, evolved exactly to $t_f=0.3$ using Eq. \ref{['eq: exact evolution']} on a $2^{12}\times 2^{12}$ grid. These results agree precisely with the direct RK4 integration of the real space equations of motion in Eq. \ref{['eq: direct EOM']}.
• Figure 5: Scaling behavior of the algorithm in two dimensions for a Ricker wavelett initial condition. The left y-axis shows the time associated with each operation named in the first four entries of the legend, where 'Prop.' refers to the time for propagating the Exact or MPS solution, equations \ref{['eq: exact evolution']} and \ref{['eq: Trotter MPS Evol']} respectively, and FFT/QFT refer to the Fast Fourier Transform and Quantum Fourier Transform algorithms, respectively. The right y-axis shows the absolute error between the MPS and exact solution: $|\psi_0^{\text{MPS}}-\psi_0^{\text{exact}}|$ with a maximum value of $0.86\%$. The bottom x-axis shows the number of qubits for a given spatial direction. The top x-axis shows the log base 10 of corresponding total number of grid points, i.e. $\text{log}_{10}(2^{2n})$. Past $n=15$ on a side we no longer include the exact calculation. The MPS algorithm begins to be faster than direct representation on the grid for grids of size $2^{14}$ on a side due to the FFT versus QFT scaling, although alternative MPS propagation strategies could potentially improve this. See the text for details.
• Figure 6: The tapered three dimensional Gaussian function for $\mu=0.5$ and $\sigma=0.1$, evolved exactly to $t_f=0.3$ using Eq. \ref{['eq: exact evolution']} on a $2^{8}\times 2^{8}\times2^{8}$ grid, sliced through the $z=0$, $xy$ plane. At the $z=0$ plane we expect the errors due to the approximate initialization of $\psi_1$ to be minimal as $\partial_z\psi_0=0$ at $t=0$ on this plane. These results agree precisely with the direct RK4 integration of the real space equations of motion in Eq. \ref{['eq: direct EOM']}.