Can a quantum homotopy algorithm solve nonlinear flow PDEs efficiently?

Quantum operations are linear, flow equations are not

Quantum gates are linear and unitary, while the Navier-Stokes equations governing turbulence and other flow problems are nonlinear. That mismatch is the central obstacle the paper attacks: a solver built on quantum linear-systems algorithms needs the nonlinear PDE rewritten as a linear problem before any quantum speedup can be claimed.

The paper's route is homotopy analysis. Instead of discretizing the nonlinear equation directly, it builds a sequence of deformation equations that start from a solvable linear guess, such as the linear diffusion solution for Burgers, and progressively deform it toward the nonlinear solution. Each deformation term is computed from the previous ones, so the nonlinearity is handled in the construction of the embedding rather than inside the quantum circuit.

New variables close the system without truncation error

The deformation equations are recursive: the M-th term depends on products of earlier terms, which would normally force an extra approximation. The authors remove that recursion by dimension dilation, defining auxiliary variables for the nonlinear products and their derivatives. For the Burgers example at order M = 4, the lifted state contains v1 through v4 plus v5 = ū1², v6 = ū1ū2 and v7 = (∂xū1)², alongside the known base profile ū0 and its derivatives.

Because the base profile is retained in the lifted state, terms that would otherwise act as forcing become linear couplings from that base. The result is a closed, homogeneous linear system dv/dt = Av with no additional truncation error, which the authors contrast with Carleman and Koopman embeddings where such truncation is intrinsic. After discretization on Ng grid points, the M = 4 Burgers operator has dimension 8Ng and O(Ng) nonzero entries.

Finite differences feed a compact time-marching circuit

Spatial derivatives are replaced by central finite differences, so the first-derivative blocks scale as O(Δx⁻¹) and the second-derivative blocks as O(Δx⁻²). Time stepping uses either explicit Euler, v_{j+1} = (I + ΔtA)v_j, or implicit Euler, v_{j+1} = (I − ΔtA)⁻¹v_j. The implicit scheme is unconditionally stable for the dissipative linear part, while the explicit scheme requires a timestep restriction.

The quantum routine advances this linear system by decomposing the step operator into a linear combination of unitaries and either expanding the sum to the power τ or iterating the circuit τ times with a countdown register of log(τ) qubits. Matrix inversion in the implicit and one-shot variants is approximated by a truncated Neumann series, whose truncation error is bounded by O(κ − 1)^{P_min} and requires P_min = O(log(1/ε_N)/log(1/(κ − 1))) terms. The paper reports that its complexity estimates improve on existing approaches in time-marching system size, simulation time, accuracy parameters and post-selection parameters.

A Reynolds-like number caps the usable time window

The paper's most decision-relevant quantity is Re_H, a homotopy nonlinearity parameter built from the base-profile amplitude, the domain length and the ratio of derivative-operator norms, evaluated at the grid spacing needed to resolve the smallest physical scale. For one-dimensional Burgers flow, where that scale is the viscous shock width, Re_H scales as the ordinary Reynolds number Re.

Crucially, convergence does not require Re_H < 1. What must stay below one is the accumulated finite-time deformation, giving the condition Re_H < 2/(|ĥ|t*) and a nonlinear time window t_NS = (L/Ū)·2/(|ĥ|Re_H). In dimensional terms the window shrinks roughly as Re⁻¹. So a higher Reynolds number does not make the method fail immediately; it shortens the interval over which the homotopy truncation is guaranteed to converge geometrically, and the error falls as ρ̄_ĥ^{M+1}/(1 − ρ̄_ĥ) while that condition holds.

Benchmarks span cubic reactions and fourth-order stiffness

The embedding is illustrated numerically on Burgers, Fisher–KPP, damped Kuramoto–Sivashinsky and real Ginzburg–Landau/Allen–Cahn equations. Together these cover quadratic transport, cubic nonlinearity and saturation, and fourth-order dissipative stiffness, so the demonstration is not confined to a single nonlinear form.

For the real Ginzburg–Landau/Allen–Cahn case, the reaction time t_H = 1/α and interfacial length ℓ_H = √(ν/α) define Re_H^{RGL} = βU_H²/α, which is O(1) at the normalized saturated amplitude. With α = 2 and ĥ = 0.25 the estimated nonlinear time is t_NS ≈ 4, well beyond the final simulation time T = 1.2, so this benchmark sits inside the admissible boundary. The final-time mean-squared error decays approximately exponentially with homotopy order M, and at M = 50 the homotopy solution tracks the direct numerical simulation over the full window.

Accuracy, order and window trade off against each other

The practical reading is a three-way trade: for fixed ĥ and non-dimensional time t*, raising the homotopy order M reduces truncation error geometrically, but only while Re_H stays under 2/(|ĥ|t*). Push the integration window longer or the nonlinearity higher and the guaranteed convergent regime ends, even though Re_H itself may exceed one.

The paper frames this as potential rather than demonstrated advantage on hardware: it presents bounds on stability, accuracy, gate counts and query complexity, and positions the approach as a hybrid candidate for near-term and fault-tolerant devices. The numerical evidence is classical simulation of the embedding, not a quantum run, so the complexity claims rest on the derived bounds rather than measured quantum performance.

The priority the paper establishes is therefore not raw asymptotic speedup but the shape of the admissible region: a finite-time, accuracy-dependent window governed by Re_H, with the homotopy order as the knob that buys accuracy inside it. Readers comparing this with Carleman-based quantum solvers should compare that window and its scaling with Re, not just the truncation order.

Sources used in this answer