Provably Efficient Quantum Algorithms for Solving Nonlinear Differential Equations Using Multiple Bosonic Modes Coupled with Qubits

TL;DR

The paper introduces a provably efficient continuous-variable quantum algorithm for nonlinear PDEs by lifting classical nonlinear dynamics to linear evolution in an enlarged KvN Hilbert space realized with multimode bosons and qubits. Dynamics are implemented as CPTP Kraus maps compiled into a binary-tree circuit, yielding logarithmic circuit depth per time step and a Kraus-rank–driven cost that scales as $\mathcal{O}\left(T(\log L + d r \log K)\right)$ for a $d$-dimensional grid, derivative order $K$, and polynomial nonlinearity of degree $r$. Demonstrations on the 1D Burgers’ equation and 2D Fisher–KPP equation validate mean-field updates and first-order fluctuations against analytic predictions, with extensive sampling ($N=10^4$ shots per grid point) supporting the claimed scaling and accuracy. The work also analyzes photon-loss resilience, introducing a counterterm to counteract dominant noise and outlining mitigation strategies, positioning KvN-based continuous-variable methods as a viable route to near-term quantum speedups in nonlinear dynamical simulations on analog hardware.

Abstract

Quantum computers have long been expected to efficiently solve complex classical differential equations. Most digital, fault-tolerant approaches use Carleman linearization to map nonlinear systems to linear ones and then apply quantum linear-system solvers. However, provable speedups typically require digital truncation and full fault tolerance, rendering such linearization approaches challenging to implement on current hardware. Here we present an analog, continuous-variable algorithm based on coupled bosonic modes with qubit-based adaptive measurements that avoids Hilbert-space digitization. This method encodes classical fields as coherent states and, via Kraus-channel composition derived from the Koopman-von Neumann (KvN) formalism, maps nonlinear evolution to linear dynamics. Unlike many analog schemes, the algorithm is provably efficient: advancing a first-order, $L$-grid point, $d$-dimensional, order-$K$ spatial-derivative, degree-$r$ polynomial-nonlinearity, strongly dissipative partial differential equations (PDEs) for $T$ time steps costs $\mathcal{O}\left(T(\log L + d r \log K)\right)$. The capability of the scheme is demonstrated by using it to simulate the one-dimensional Burgers' equation and two-dimensional Fisher-KPP equation. The resilience of the method to photon loss is shown under strong-dissipation conditions and an analytic counterterm is derived that systematically cancels the dominant, experimentally calibrated noise. This work establishes a continuous-variable framework for simulating nonlinear systems and identifies a viable pathway toward practical quantum speedup on near-term analog hardware.

Figures (6)

• Figure 1: Overview of the continuous-variable quantum algorithm. Left: A nonlinear field evolution $\dot z=V\left(z(t)\right)$ is lifted by the KvN map to a linear operation generated by $A$ (which encodes the differential equation information) on an enlarged space. Middle: Each short step is implemented as a local CPTP map $K_a=e^{-A\Delta t}$ acting on bosonic modes prepared in multimode coherent states that encode the classical field. Right:$K_a$ is compiled into a measurement-adaptive binary-tree circuit: a system–ancilla unitary $U$, measurement and reset at each layer, and post-selection on $\ket{0}$ yield logarithmic circuit depth in the Kraus rank.
• Figure 2: Circuit realization and unitary construction for a 1D Burgers' equation. Left: Measurement–adaptive binary-tree realization of an $N$-rank-Kraus-operator CPTP map using a single ancilla qubit Shen_2017. At depth $\ell$, the joint unitary $U_{\mathbf{b}^{(\ell-1)}}$ is applied, where $\mathbf{b}^{(\ell-1)}$ records the past measurement outcomes. After each interaction, the ancilla is measured and reset. Post-selection on $\ket{0}$ (empty-circle control) singles out the Kraus operator $K_a$. Right: Magnitude plots of all compiled block unitaries $U_{\mathbf{b}}$ for a four-mode bosonic system truncated to four levels per mode ($d_{\mathrm{sys}}=256$), arranged hierarchically from the root ($U_{\varnothing}$) to depth-4 leaves ($U_{0000}, \dots, U_{1111}$ with subscripts in the lexicographical order). For each $U_{\mathbf{b}}$, the left half encodes the analytically derived Kraus block $\langle b_{\ell+1} | U_{\mathbf{b}} | 0 \rangle$, while the right half shows the QR-completed orthonormal complement required to form a valid unitary.
• Figure 3: 1D Burgers' validation.Left: mean-field solution profiles at selected times. The results demonstrate the rightward advection with speed set by $u$, nonlinear steepening from the advective term $-u {\partial u}/{\partial x}$ together with diffusive broadening from ${1/R_e}\cdot{\partial^2 u}/{\partial x^2}$. Middle: spacetime map of the mean field $u(x,t)$ over the entire evolution window, showing advection-dominated rightward drift with weak viscous spreading. Right: Sample-mean bias using $N=10^4$ shots per grid point. The residuals are centered around zero and spatially unstructured, consistent with the prediction that $\mathrm{Var}_j(t)$ remains constant and time-independent to first order [Eq. \ref{['eq:var-burgers-main']}], so that sampling fluctuations are the only visible deviations.
• Figure 4: 2D Fisher-KPP validation. Row 1: analytic mean fields $u(x,y,t)$ at representative times $t=0.02,\,0.20,\,0.40,\,0.80$ s, showing the rotational advection imposed by the velocity field, the smoothing action of diffusion ($\mathrm{Pe}=200$), and the amplitude saturation from the Fisher reaction term ($\mathrm{Da}=1$). Row 2: Bias maps from $N=10^4$ trajectories. The residuals remain centered around zero with no coherent structure, confirming that the stochastic sampling is unbiased. The omitted panels for the predicted $1\sigma$ noise width and the empirical-to-analytic error-bar ratio behave as in Supplementary Materials \ref{['app:2D-fisher']}, i.e. Gaussian fluctuations with variance $V_{p,q}(t)$ that scale as $\sigma/\sqrt{N}$. Together, the results validate the analytic update formulas \ref{['eq:fisher2d-mean']}--\ref{['eq:fisher2d-var']} across the evolving vortex and front structures.
• Figure 5: The lid-driven cavity simulation with a $128\times128$ grid and $\mathrm{Re} = 1000$ in top row and a $256\times256$ grid and $\mathrm{Re} = 10000$ in bottom row. All panels show the flow at the steady state. Left: Contour plots of the stream function $\psi$. Middle: The velocity field $\mathbf{u} = (u,v)$. Right: The vorticity $\omega$. Simulation is done using bosonic simulator. The evolution follows the Euler update rule, and therefore, matches with DNS calculation.

Theorems & Definitions (2)

• Theorem 4.1: Kraus Rank and Circuit Depth for Linear Generators
• Theorem 4.2: Kraus Rank and Circuit Depth for Polynomial Generators