The quantum-simulation neighborhood of the neutrino Vlasov algorithm
A survey of quantum simulation — using a quantum machine to compute the behaviour of a physical system. We map the field on two maps, both organised the same way: rows = the physical phenomenon, columns = the simulation technique, and a bold stamp on every paper saying whether it is only proposed or has actually been run. Every card links to that paper's full slide.
Map 1 · near neighbourhood — using a quantum computer to solve our kind of problem: the classical kinetic / Vlasov PDE behind the neutrino seed (plasma, dark matter, fluids). Almost all still proposed. Map 2 · wider landscape — quantum machines simulating quantum systems and analogs (gauge theory, chemistry, materials, black-hole physics); much of it has already run on real hardware.
Put Maps 1 and 2 side by side and the whole story appears: Map 1's hardware column is nearly empty (aspirational), while Map 2's is full (real today). The seed we build from is Miyamoto+24 — the next slide.
The seed: Miyamoto, Yamazaki, Mori & Yoshida (2024)
The anchor. A fault-tolerant quantum algorithm that evolves the linearized massive-neutrino Vlasov equation to predict the neutrino contribution to the matter power spectrum \(P_\nu(k)\) — with an exponential memory saving and a near-linear-in-resolution query cost.
Massive neutrinos free-stream out of gravitational wells and suppress small-scale structure by \(\Delta P/P \approx -8\,\Omega_\nu/\Omega_m\). Capturing this requires the full 6-D phase-space distribution \(f(\mathbf x,\mathbf v,t)\), whose grid cost scales as \(O(n_{\rm gr}^6)\) — the curse of dimensionality that makes a direct classical Vlasov solve hopeless at cosmological resolution.
(external \(\mathbf F_{\rm CDM}\))→ discretize
\(\dot{\mathbf f}=A(t)\mathbf f\)→ \(H=iA\)
Hermitian→ Schrödinger
+ amplitude encoding→ Ham-sim
+ QAE read-out
Neglecting neutrino self-gravity (\(\Omega_\nu/\Omega_m\lt 1\%\)) makes the equation linear; a linear, norm-preserving evolution is exactly a Schrödinger equation, which a quantum computer runs natively. The 2\(^n\)-dimensional state vector holds the whole phase space in \(n=O(\log)\) qubits — the exponential memory win.
Complexity: classical \(O(n_{\rm gr}^6 n_t)\) → quantum \(\widetilde O(n_{\rm gr}+n_t)\) queries with \(O(\log^{5/2}(n_{\rm gr}/\epsilon))\) qubits. It owes the survey's three taxes — and because it linearizes, only two bite: QRAM (Q) to load the CDM force, and read-out (R) capped at \(P_\nu(k)\) at \(1/\epsilon\); the nonlinearity tax (N) is zero. What it pays instead is maturity — deep fault-tolerant circuits. Those three taxes, Q · R · N (QRAM · read-out · nonlinearity), recur across the neighborhood as the spine of this survey. arXiv:2310.01832 · the deep 48-slide walkthrough is the Vlasov Tutorial deck.
Basic definitions — how to read every map
Three things fix where a paper sits and how real it is. Learn these and you can read all three maps.
Simply the system whose behaviour you want. Map 1: cosmology & dark matter · plasma kinetics · fluids & classical fields. Map 2: gauge theories & HEP · condensed matter & materials · chemistry & nuclear · gravity, black-hole & vacuum analogs.
Hamiltonian simulation — rewrite the problem as a Schrödinger equation and run \(e^{-iHt}\). Carleman / linear-systems — turn a nonlinear problem into linear algebra, then solve it. Variational / analog / other — shallow NISQ circuits, or a native/analog quantum system (cold atoms, spins, lattices). On Map 2, digital gate-based — Trotterised circuits on a gate quantum computer.
proposed needs a future fault-tolerant machine (theory + resource estimate) · NISQ / sim ran on a small device or a classical emulator · ✓ ran on hardware a real quantum-hardware result. On Map 1 only, three extra badges Q·R·N flag the costs a PDE algorithm pays: QRAM (loading data), Read-out (getting the answer), Nonlinearity.
Map 1 · the near neighbourhood — phenomenon × technique
Map 1 · the near neighbourhood. Using a quantum computer to solve our kind of problem — the classical kinetic / Vlasov PDE behind the neutrino seed. Rows = phenomenon, columns = technique, card colour = how real it got. Click any card for its full slide. Almost the whole map is proposed — algorithms aimed at a future quantum computer (green outline), classically validated at most; only QVLASOV and Sundar have actually run on one.
neutrinos · structure · fuzzy DM
Vlasov · Landau damping · instabilities
Navier–Stokes · waves · EM
Map 2 · the wider landscape — phenomenon × technique
Map 2 · the wider landscape. A little farther out: quantum machines simulating quantum systems and analogs — gauge theory, matter, chemistry, black-hole physics. Same layout (phenomenon × technique); card colour is how real, platform tag names the device. The headline: this map is almost all analog quantum and quantum computer — where quantum simulation is already real, the mirror image of Map 1.
Schwinger model · lattice gauge · scattering
Hubbard · magnets · time crystals
molecules · light nuclei (VQE)
Hawking · false vacuum · scrambling
Engel–Smith–Parker 2019 — plasma Landau damping
This is the paper that started the field: the first demonstration that a kinetic plasma equation can be recast as a Hamiltonian-simulation problem suitable for a future fault-tolerant quantum computer. Engel, Smith & Parker, “Quantum algorithm for the Vlasov equation,” Phys. Rev. A 100, 062315 (2019) [arXiv:1907.09418]. Every later slide in this family — and the Miyamoto cosmology seed — is a descendant of the reduction introduced here.
The target is a collisionless electron plasma described by the Vlasov–Poisson system in one spatial and one velocity dimension. The authors linearize about a Maxwellian background \(f_0(v)\propto e^{-v^2/2v_{\rm th}^2}\), writing \(f=f_0+\delta f\) with \(|\delta f|\ll f_0\). The physics they aim to reproduce is Landau damping: the collisionless decay of a small electrostatic wave as it phase-mixes among particles of different velocities. This is the canonical test of any kinetic solver — a purely kinetic effect with no fluid analogue — so choosing it as the benchmark is a deliberate statement that the quantum method captures genuine phase-space physics, not just advection.
Linearizing kills the quadratic \(E\,\partial_v f\) self-interaction, leaving a linear PDE for \(\delta f\) with the self-consistent field determined by Poisson’s equation \(\partial_x E = -\!\int\delta f\,dv\). They expand in Fourier modes in \(x\) and discretize velocity, so the continuous PDE becomes a finite linear ODE system
where \(A\) encodes streaming \(v\,\partial_x\), the field coupling, and the resonant velocity structure. The essential trick — inherited by the whole family — is to arrange the discretization so that \(A\) is (anti-)Hermitian: then \(H=iA\) is a legitimate Hamiltonian and the classical ODE becomes a Schrödinger equation \(i\,\partial_t|\delta f\rangle = H|\delta f\rangle\).
The \(N_v\)-component state vector is stored in the amplitudes of \(\lceil\log_2 N_v\rceil\) qubits — the exponential-memory compression that motivates the whole program. Because \(H=iA\) is sparse (streaming and field couplings connect only neighbouring grid/Fourier indices), one applies an optimal sparse-Hamiltonian-simulation routine to realize the propagator \(e^{-iHt}\) as a circuit. A physical observable — the field energy, or a low-order velocity moment of \(\delta f\) — is then read out by amplitude estimation. The unitary evolution conserves the state norm, which is exactly consistent with the linearized Vlasov dynamics being norm-preserving (phase mixing moves energy between resolved modes rather than destroying it).
They designed the full algorithm and verified it numerically on Landau damping, recovering the correct damping rate. The headline complexity contrasts classical cost \(O(N_v\,t)\) with the quantum query cost
where \(\delta\) is the target accuracy on the extracted observable. The \(\mathrm{polylog}(N_v)\) is an exponential improvement in the velocity-resolution dependence — the payoff of amplitude encoding. The authors are explicit that the \(1/\delta\) factor, coming from the amplitude-estimation readout, is the practical bottleneck: you pay it for every digit of accuracy on the answer, and it does not shrink with better grids.
This is the template the Miyamoto cosmology paper specialises: replace the plasma’s self-consistent electric field with an external CDM gravitational force, and the same “antisymmetric \(A\) ⇒ \(H=iA\) ⇒ Hamiltonian simulation ⇒ amplitude-estimation readout” chain carries over to neutrino structure formation. What is solid here: the reduction is exact for the linearized system, and the exponential grid compression is real. What is assumed/limited: the plasma is linearized (the genuinely hard nonlinear \(E\,\partial_v f\) term is dropped); Landau damping is captured only up to the recurrence time set by finite \(N_v\); readout costs \(1/\delta\); and the circuits presume a future fault-tolerant machine, not NISQ hardware. It is a proof of principle, honestly labelled as such. arXiv:1907.09418
Yamazaki et al. 2023 — self-gravitating collapse
Yamazaki, Uchida, Fujisawa, Miyamoto & Yoshida, “Quantum algorithm for collisionless Boltzmann simulation of self-gravitating systems,” Comput. Fluids 288 (2025) [arXiv:2303.16490]. From the same group as the Miyamoto seed (Koichi Miyamoto is a co-author), this paper takes the crucial step the seed deliberately avoids: it makes gravity self-consistent, extracting the potential from the quantum state itself at every timestep and feeding it back — closing the Poisson loop inside the quantum simulation.
The target is the collisionless Boltzmann equation (the Vlasov equation) for a system bound by its own gravity — star clusters, galaxies, and the cosmological growth of structure. Unlike the neutrino seed, where CDM supplies an external force, here the matter that streams is the same matter that sources the potential. The showcase instability is Jeans collapse: a density perturbation whose self-gravity overcomes pressure support and grows, the fundamental mechanism of gravitational structure formation. Reproducing its growth rate is the physics validation.
They solve the collisionless Boltzmann equation coupled to Poisson,
on a phase-space grid, building on Todorova & Steijl’s (2020) quantum free-streaming scheme. The free-streaming and force operators are discretized so the advection part is realized as unitary Hamiltonian evolution — the same antisymmetric-generator structure as the rest of this family. The new ingredient is that \(\Phi\) is no longer prescribed: it must be recomputed from \(f\) each step.
The key idea is to extract the potential without collapsing the whole \(2^n\)-amplitude state. Because a gravitational potential is dominated by its long-wavelength (major Fourier) modes, they estimate only those few dominant modes of the density \(\rho(\mathbf x)=\int f\,d^3v\) from the solution-encoding state — via quantum Fourier / amplitude estimation on the position register — reconstruct \(\Phi\) classically from those modes through \(\Phi_{\mathbf k}=-4\pi G\,\rho_{\mathbf k}/k^2\), and then apply the corresponding force operator for the next timestep. This is a hybrid inner loop: unitary streaming inside each step, a small classical Poisson solve between steps driven by a handful of measured modes. Critically, it improves the dependence on \(N_v\) (the velocity grid size) relative to classical solvers — the velocity axis is where the classical curse of dimensionality bites hardest, so compressing it is the main win.
Two concrete numerical tests, run on a state-vector emulator: (i) a 1+1-dimensional free-streaming calculation on a \(64\times64\) phase-space grid using 13 simulated qubits, verifying the advection scheme; and (ii) a Jeans-collapse run whose growth was compared against analytic linear-perturbation theory and found to agree. Together these show the closed gravity loop actually works: the self-consistent field, reconstructed from a few measured Fourier modes each step, reproduces the textbook instability.
This is the seed’s self-consistent sibling. The neutrino paper keeps the force external and therefore stays fully linear and fully coherent — no measurement until the very end — which is what buys its clean exponential speedup. Yamazaki et al. instead reinstate the nonlinear gravity coupling, and the price is honest and explicit: you must measure the state every timestep to get the density modes, and you must truncate to the “major” Fourier modes, an approximation that discards small-scale potential structure. Per-step measurement and re-preparation erode the coherent advantage and reintroduce a shot-noise-like cost. What is solid: the loop closes, the \(N_v\)-scaling improves, and Jeans collapse is reproduced. What is assumed: only a few modes matter, the emulation stands in for real hardware, and — as throughout the family — an eventual fault-tolerant machine is required. arXiv:2303.16490
Ameri et al. 2023 — collisional plasma kinetics
Ameri, Cappellaro, Krovi, Loureiro & Ye, “A quantum algorithm for the linear Vlasov equation with collisions” [arXiv:2303.03450, Phys. Rev. A 107, 062412 (2023)]. This paper adds the ingredient every real plasma has — collisions — which turn the evolution non-unitary, and it is unusually candid about the size of the quantum advantage: a careful quadratic speedup, not an exponential one, once you compare against an efficient classical spectral method.
The target is the 1D electrostatic linearized Vlasov equation, studied both collisionless and collisional. Collisions are modelled by a Fokker–Planck–type (Lenard–Bernstein / Dougherty) operator that relaxes the distribution toward a Maxwellian and, crucially, provides real dissipation in velocity space. This matters because collisionless Landau damping is only reversible phase mixing (it recurs), whereas physical collisions smooth fine velocity structure irreversibly — the regime relevant to plasma oscillations, transport, and stability.
Instead of a phase-space grid, they use a Fourier–Hermite expansion — Fourier modes in \(x\), Hermite functions \(\psi_m(v)\) in velocity:
This basis is exponentially more compact for smooth distributions than a naive velocity grid, and it diagonalises the collision operator: the Hermite index \(m\) is the “number of velocity wiggles,” and the Lenard–Bernstein operator damps mode \(m\) at rate \(\propto -\nu\,m\). The advection term becomes a tridiagonal ladder coupling \(m\leftrightarrow m\pm1\). The result is a linear ODE \(\dot{\mathbf C}=(H_{\rm adv}+C_{\rm coll})\mathbf C\) where the collisional part is anti-Hermitian-plus-negative-definite — i.e. genuinely non-unitary.
The two regimes need two different primitives. Collisionless: the generator is anti-Hermitian, so the distribution is amplitude-encoded and evolved by standard Hamiltonian simulation of \(e^{-iHt}\), exactly as in Engel–Smith–Parker. Collisional: the dissipative generator is not Hermitian, so \(e^{At}\) is contractive, not unitary; they invoke a quantum linear-ODE solver that block-encodes the non-unitary propagator (a linear-combination-of-unitaries / block-encoding construction) and recovers the correctly normalised output with amplitude amplification. The dissipation actually helps the encoding: a contractive map keeps the state inside the encodable subspace and suppresses the high-\(m\) recurrence that plagues the collisionless case.
They establish that a quadratic speedup in system size (the number of Fourier–Hermite modes) is attainable for both cases. The sharpest and most valuable result is a piece of intellectual honesty: they show that a classical algorithm working in the same Fourier–Hermite representation can already match the performance of a previously proposed quantum algorithm for this problem. In other words, once the classical method is allowed the same compact basis, the earlier “exponential” framing evaporates and the defensible quantum advantage is quadratic — the generic Grover-type square-root, not a Hilbert-space-dimension exponential.
Relative to the seed, this paper is the reality check of the family. It shares the linear-Vlasov, amplitude-encoding backbone, but (a) adds collisions, moving beyond the purely unitary picture the seed relies on, and (b) insists on benchmarking against a strong classical spectral baseline rather than a naive grid solver. The lesson for reading the seed’s exponential claim: the exponential is against a brute-force \(O(n_{\rm gr}^6)\) grid; against a smart classical basis the honest margin shrinks. What is solid: the non-unitary ODE handling and the quadratic bound. What is limited: strictly 1D and linear, no self-consistent nonlinearity, and — as always — a fault-tolerant device is assumed. arXiv:2303.03450
Toyoizumi et al. 2024 — plasma Landau damping (QSVT)
Toyoizumi, Yamamoto & Hoshino, “Hamiltonian simulation using quantum singular value transformation: complexity analysis and application to the linearized Vlasov–Poisson equation” [arXiv:2304.08937]. This paper swaps the simulation engine: instead of Trotterization or generic sparse-Hamiltonian methods, it drives the evolution with the quantum singular value transformation (QSVT) — the modern, query-optimal primitive — and delivers a rigorous end-to-end query count for simulating Landau damping.
The target is the 1D linearized Vlasov–Poisson system, and the physics probed is once more linear Landau damping — the collisionless decay of an electrostatic wave by phase mixing among resonant particles. Using the same well-understood benchmark as Engel–Smith–Parker is deliberate: it isolates the contribution of this paper, which is not new physics but a better and provably-costed quantum method for the same reduction.
The linearized Vlasov–Poisson operator is discretized (Fourier in space, a velocity discretization) into the familiar antisymmetric \(A\), giving \(H=iA\) and the Schrödinger form \(i\,\partial_t|f\rangle=H|f\rangle\). QSVT requires the operator to be presented as a block-encoding: a larger unitary \(U_H\) whose top-left block is \(H/\alpha\) for some normalisation \(\alpha\ge\|H\|\). Sparsity of \(A\) makes this block-encoding cheap to build. The whole simulation is then reduced to applying a polynomial of the block-encoded \(H\).
QSVT realizes the propagator by approximating the complex exponential with a polynomial — the Jacobi–Anger expansion
and implementing that polynomial of \(H\) directly on the singular values of the block-encoding via a tuned sequence of phase rotations. Because QSVT’s output is a sub-normalised block, the near-unitary evolution is restored with amplitude amplification. The authors carry out a detailed error-and-query analysis comparing two variants and conclude that oblivious amplitude amplification outperforms fixed-point amplification “in the sense of simulation time,” because the oblivious version boosts success without needing to reflect about the (unknown) target state. QSVT is known to achieve the optimal query complexity for Hamiltonian simulation — linear in \(t\) and logarithmic in \(1/\varepsilon\), saturating the no-fast-forwarding lower bound.
They successfully simulated linear Landau damping within the QSVT framework, backing it with an extensive numerical study (sixteen figures) of the error behavior and the query counts under the two amplification strategies. The deliverable is thus twofold: a working QSVT Vlasov–Poisson simulation, and a transferable, rigorously-costed recipe for how many oracle queries the evolution actually takes as a function of \(t\), \(\varepsilon\), and the block-encoding normalisation \(\alpha\).
This is the seed’s reduction executed with a sharper tool. Miyamoto’s neutrino algorithm leans on block-encoding plus (implicitly) Trotter-style evolution and amplitude estimation; Toyoizumi et al. show that QSVT gives the query-optimal way to realize the very same \(e^{-iHt}\), with the amplitude-amplification bookkeeping made explicit. What is solid: the complexity analysis is careful and the method is optimal for the linear problem. What is limited: strictly 1D and linearized (the Vlasov–Poisson coupling is kept but linearized, so no genuine nonlinearity), Landau damping is the only demonstration, and the results are emulated under a fault-tolerant assumption. It advances the engineering of the family, not its physical reach. arXiv:2304.08937
Higuchi et al. 2025 — two-stream plasma instability
Higuchi, Pedersen, Toyoizumi, Yoshikawa, Kiumi & Yoshikawa, “Quantum calculation for two-stream instability and advection test of Vlasov–Maxwell equations,” J. Plasma Phys. 91, E116 (2025) [arXiv:2408.11550]. This is the first paper in the family to tackle a genuine kinetic instability with a self-consistent field, and it does so pragmatically: evolve the Vlasov part on a quantum circuit, solve the fields classically, and hand data back and forth each step.
The target is the electrostatic Vlasov–Maxwell system, tested on 1D advection and on the 1D1V two-stream instability — two counter-streaming particle populations whose relative drift feeds a growing electrostatic wave, the archetypal kinetic instability and a standard PIC/Vlasov benchmark. Unlike Landau damping (which decays), this mode grows, so it stresses a solver’s ability to track exponentially amplifying, self-consistent field–particle feedback. Classically the obstacle is the full 6D phase space spanning wide spatial and temporal scales.
They partition the coupled system by physics: the Vlasov advection of \(f\) in phase space is the quantum workload, while Maxwell/Poisson field updates from the charge and current densities are the classical workload. On an \(N\)-point spatial grid over time \(T\), the classical direct cost is quoted as
the reduction coming from amplitude-encoding the phase-space advection so that the \(6\)-D grid enters only logarithmically.
Each timestep runs a three-part loop: (1) a QSVT-based Hamiltonian simulation advances the amplitude-encoded distribution \(f\) under the streaming/force operator; (2) the state is measured to reconstruct charge and current densities; (3) a classical Maxwell solver updates the fields, which are re-encoded to define the next step’s Hamiltonian. A notable advantage is the absence of a CFL condition on the quantum advection: classical explicit schemes must obey \(\Delta t \lesssim \Delta x/v_{\max}\), but the unitary Hamiltonian evolution is stable under much larger timesteps, so the quantum substep can leap over the Courant limit that throttles classical solvers. The circuits were run on a Qiskit-Aer-GPU state-vector emulator on an NVIDIA A100.
They validated the 1D advection test for transport accuracy, and then reproduced the growth of the 1D1V two-stream instability on the emulator — the field feedback correctly drives the exponential amplification of the electrostatic perturbation. This is the family’s first end-to-end demonstration of a self-consistent, unstable kinetic mode via the hybrid route, and it confirms that the no-CFL property holds in practice for the advection substep.
Where the seed keeps the force external to preserve full coherence and a clean speedup, Higuchi et al. accept a hybrid loop to gain a self-consistent field and an instability the seed cannot represent. The trade-off is the central honest caveat: the loop requires measuring the quantum state and re-preparing it every timestep, and each measurement demands many shots to beat statistical noise — precisely the per-step readout cost the seed was engineered to avoid. So the polynomial advantage is asymptotic, with an unclear practical crossover, and the frequent quantum↔classical communication can dominate real runtime. What is solid: the no-CFL advection and the reproduced instability. What is assumed: GPU emulation stands in for hardware, fault-tolerant gates underlie the complexity bounds, and readout overhead is not yet costed at scale. arXiv:2408.11550
Berntson et al. 2026 — weakly-nonlinear plasma
Berntson, Jennings, Lostaglio & Parker (PsiQuantum), “An end-to-end quantum algorithm for weakly nonlinear plasma physics with superquadratic speedup” [arXiv:2607.14308]. This is the most ambitious descendant in the family: it keeps the nonlinearity (via Carleman linearization), gives a full end-to-end logical resource analysis (logical gate and qubit counts with rigorous convergence proofs — not yet fault-tolerant T-counts) from problem statement to observable, and — closing a seven-year arc — carries Scott Parker’s name from the 2019 origin paper.
The target is a 3D electron–ion plasma with adiabatic electrons, kinetic ions, Debye screening, and a Krook collision (relaxation) model — a realiztic kinetic system sitting between full Vlasov–Poisson and reduced models. The dynamics are weakly nonlinear: the authors use the plasma free energy to build a Lyapunov transform under which the nonlinear terms stay bounded, defining a certified regime in which the whole analysis is rigorous. This is the family’s first serious move past strict linearization into controlled nonlinear plasma dynamics.
A Fourier–Hermite truncation (cutoffs \(N_F\) in the Fourier modes, \(N_H\) in the Hermite index) reduces the kinetic dynamics to a high-dimensional quadratic ODE,
where the quadratic term \(B\) is the plasma nonlinearity. The Lyapunov/free-energy structure guarantees this quadratic system stays weakly nonlinear over the simulated horizon — the condition that will make the linearization below converge.
The quadratic ODE is embedded into an (infinite, then truncated) linear system by Carleman linearization — stacking tensor powers \(\mathbf u,\ \mathbf u^{\otimes2},\ \mathbf u^{\otimes3},\dots\) so the quadratic map becomes linear on the enlarged space. They prove this Carleman truncation converges exponentially in the truncation order within the certified weakly-nonlinear regime — the load-bearing theorem of the paper. The resulting linear generator is dense, which would normally kill the query cost; they overcome this with a hierarchical block-encoding for dense matrices that exploits the spatial decay of the field, avoiding the generic dense-access penalty. Evolution then proceeds by linear-ODE / Hamiltonian-simulation on the Carleman system, and a bespoke extraction subroutine reads observables out of the Carleman history state, using the encoded nonlinear components to sharpen estimates of linear functionals of \(f\).
The end-to-end gate and qubit costs are
Two advantages follow. The logarithmic qubit count is an exponential memory saving over classical Fourier–Hermite solvers, which must store all \(\sim N_F N_H\) coefficients. And the gate count scales as \(\sqrt{N_H}\) rather than linearly in the Hermite dimension, which combined with the log factors yields a superquadratic time speedup — strictly better than the generic quadratic (Grover-type) improvement that limits Ameri et al. Everything comes with rigorous convergence guarantees, from the Carleman embedding through to the observable.
This paper is the family’s maturity point. Against the seed it advances on every front the seed left open: it keeps genuine nonlinearity (Carleman rather than external-force linearization), it delivers a full end-to-end logical complexity accounting instead of an asymptotic query bound, and its superquadratic speedup beats the honest quadratic ceiling the collision paper exposed. What is solid: the exponential-memory advantage, the certified-regime convergence proofs, and the dense-matrix block-encoding trick. What is honestly bounded: the certification holds only in the weakly nonlinear regime — turbulence, shocks, and strong nonlinearity are explicitly out of scope; the \(1/\varepsilon\) readout factor persists; and while “end-to-end,” the abstract stops short of surface-code T-counts and logical-qubit overheads, and there is no hardware demonstration — it establishes a controlled nonlinear benchmark, a theoretical milestone rather than a demonstrated quantum advantage. arXiv:2607.14308
Mocz & Szász 2021 — fuzzy dark matter
The seed paper (Miyamoto et al.) forces a classical kinetic PDE into Schrödinger form as a mathematical trick. Mocz & Szász start from the opposite end: for fuzzy / wave dark matter the Schrödinger–Poisson system is the actual physics, so a quantum computer is not a re-encoding of a classical solver — it is the natural substrate. Their contribution is a hybrid variational scheme that confronts the one thing quantum circuits hate: the self-gravity nonlinearity. [arXiv:2101.05821, ApJ 910, 29 (2021)]
Fuzzy dark matter (FDM) posits an ultralight boson, \(m\sim 10^{-22}\ \text{eV}\), whose de Broglie wavelength is kiloparsec-scale — so the entire galactic halo is a single coherent wavefunction. On scales above that wavelength it must reproduce cold-dark-matter clustering; below it, wave interference builds a solitonic core and suppresses small-scale structure, a candidate cure for CDM's "small-scale crises" (cusp–core, missing satellites). Simulating it means evolving a complex field \(\psi(\mathbf x,t)\) whose \(|\psi|^2\) is the mass density — genuinely a quantum-like object, unlike CDM particles.
FDM obeys the Schrödinger–Poisson (Gross–Pitaevskii-with-gravity) system,
The kinetic term is the ordinary Schrödinger Laplacian; the potential \(\Phi\) is sourced by the field's own density \(|\psi|^2\) through Poisson. That self-coupling is what makes the problem nonlinear — \(\Phi[\psi]\) depends on \(\psi\), so the Hamiltonian is state-dependent. Note the sharp contrast with the seed paper, where neutrino self-gravity was dropped (external CDM force) precisely to keep the equation linear; here the nonlinearity is the whole point and cannot be discarded.
Classically the SP system is solved by a split-step Fourier "kick–drift–kick": a half potential kick \(e^{-i m\Phi\,\Delta t/2\hbar}\) (diagonal in position), a full kinetic drift \(e^{i\hbar\nabla^2\Delta t/2m}\) (diagonal in Fourier space via an FFT), then another half kick. Each factor is unitary — ideal for a circuit — but the kick needs \(\Phi[\psi]\), i.e. reading \(|\psi|^2\) back out, which a quantum state forbids cheaply. Mocz & Szász adopt the Lubasch et al. (2020) variational quantum nonlinear processing framework: represent \(\psi\) by a shallow parameterized circuit \(|\psi(\boldsymbol\theta)\rangle\), and advance one timestep by minimizing a cost function
where the nonlinear operator \(\hat U_{\rm step}\) — including the density-dependent kick — is evaluated through overlap (Hadamard-test) circuits that compute the required expectation values without ever exposing all \(2^n\) amplitudes. The FFT becomes the quantum Fourier transform; nonlinearity is handled by re-using the same register in multiple copies inside the cost estimator rather than by measuring the density directly.
This is an early, exploratory paper: the authors present a proof-of-concept mock quantum simulation of the SP dynamics — a small-scale, classically-emulated run of the variational algorithm demonstrating that the interference and gravitational-condensation behavior of FDM can be captured by the variational circuit, not a fault-tolerant hardware run. The memory argument is the headline: \(n\) qubits encode \(2^n\) grid amplitudes, so a 3-D FDM box that costs \(N=n_{\rm gr}^3\) classical numbers needs only \(\log_2 N\) qubits — the same exponential-compression lever the seed paper rides, but here for a field that is physically a wavefunction.
Same target, opposite philosophy. Miyamoto keeps the force external to stay linear and reaches for fault-tolerant Hamiltonian simulation (block-encoding, QAE); Mocz & Szász keep the self-gravity nonlinear and reach for a NISQ-friendly variational loop. The honest costs: (i) the variational method trades exponential circuit depth for a classical outer optimization that can stall in barren plateaus as qubit count grows; (ii) every nonlinear cost evaluation needs many circuit repetitions and copies of the state, inflating shot counts; (iii) as always the output is a summary (density profile, soliton mass), not the full \(\psi\) — reading \(2^n\) amplitudes is impossible. It is a direction-setting paper: it shows FDM is the most "quantum-native" dark-matter target, while being candid that a decisive advantage awaits both better hardware and a solution to the nonlinearity/readout tax.
Cappelli et al. 2023 — dark-matter structure formation
Cappelli et al. sit exactly between the seed paper and Mocz: they take the cosmological Vlasov–Poisson problem, collapse its 6-D phase space to a 3-D Schrödinger–Poisson proxy via the Madelung correspondence, and then — unlike Miyamoto's fault-tolerant machinery — solve it with a variational real-time quantum algorithm aimed squarely at near-term devices. [arXiv:2307.06032 (preprint)]
Standard cold-dark-matter structure formation is a Vlasov–Poisson problem in \(6\text{D}+1\): the distribution \(f(\mathbf x,\mathbf v,t)\) gridded on \(n_{\rm gr}^6\) cells (the same curse the seed paper fights). The observable is again the density field and its power spectrum \(P(k)\). Cappelli et al.'s aim is a route that keeps the essential nonlinear gravitational clustering but pays only for 3 spatial dimensions, then hands that reduced problem to a quantum computer whose qubit count grows only logarithmically with resolution.
The Madelung transform writes a wavefunction \(\psi=\sqrt{\rho}\,e^{i\theta}\) whose modulus gives the density and whose phase gradient gives the fluid velocity, \(\mathbf u = (\hbar/m)\nabla\theta\). Substituting into Schrödinger–Poisson reproduces the cold-collapse (dust) dynamics up to a "quantum pressure" term \(\propto \hbar^2\nabla^2\sqrt\rho/\sqrt\rho\) that vanishes as \(\hbar/m\to 0\):
So a single 3-D complex field stands in for the full 6-D phase-space distribution: multi-streaming is encoded in the wavefunction's phase interference rather than in an explicit velocity axis. The price is a controlled proxy — faithful in the classical (small \(\hbar/m\)) limit, with the authors explicitly analyzing how the nonlinearity variance biases observables away from that limit.
Rather than compile \(e^{-iHt}\) exactly, they evolve a parameterized state \(|\psi(\boldsymbol\theta(t))\rangle\) by moving its parameters along the McLachlan/variational-principle flow that best matches true Schrödinger evolution over each step. The genuinely novel ingredient is a set of quantum circuits that link the solution of the Poisson equation to the Schrödinger step: at each timestep the density-sourced potential \(\Phi[\psi]\) is recomputed and fed back, closing the nonlinear loop on-device. Because the ansatz is shallow and the Poisson solve is expressed as circuits rather than a QRAM oracle, the scheme is NISQ-oriented — no fault-tolerant block-encoding, no amplitude estimation.
The concrete finding is an empirical logarithmic relationship between qubit count and SP-system scale: pushing spatial resolution up costs only \(\sim\log_2 n_{\rm gr}\) extra qubits, the exponential-compression win, now demonstrated for the variational proxy. They run proof-of-concept simulations (small grids, emulated) tracking the collapse toward the classical limit and quantify how the nonlinearity's variance perturbs measured observables. It is not a cosmological-scale run and does not claim a wall-clock speedup — it is a feasibility-and-scaling study.
Same target (DM structure / P(k)), three ways it differs from Miyamoto: (i) dimensional reduction 6-D→3-D via Madelung instead of a full 6-D Vlasov grid; (ii) a NISQ variational engine instead of fault-tolerant Hamiltonian simulation; (iii) the self-gravity kept nonlinear and on-device instead of linearized to an external force. Honest caveats: the Madelung proxy only equals Vlasov in the classical limit and can misrepresent genuine multi-stream crossings; VarQTE inherits barren-plateau and measurement-overhead problems that worsen with qubit count; and, as everywhere on this family, the accessible output is a low-dimensional summary. Its value to the survey is as the bridge case — it shows the Schrödinger-proxy idea (Mocz) and the cosmological Vlasov target (Miyamoto) can be married on near-term hardware, at the cost of an approximation the seed paper never had to make.
Engel–Smith–Parker 2021 — nonlinear Vlasov plasma
The seed paper avoided nonlinearity by fiat (external CDM force). Engel, Smith & Parker confront it head on: they ask how a genuinely nonlinear Vlasov–Poisson plasma could ever be put on a quantum computer, whose gates are strictly linear. Their answer is Carleman linearization — trade nonlinearity for infinite dimension, then truncate. [arXiv:2012.06681, Phys. Plasmas 28, 062305 (2021)]
Electrostatic Vlasov–Poisson for a plasma is the self-consistent version of the seed problem: the force is sourced by the plasma's own charge density, \(\mathbf E = -\nabla\Phi\), \(\nabla^2\Phi=-\rho[f]/\epsilon_0\), so the Vlasov term \(\mathbf E\!\cdot\!\partial_{\mathbf v} f\) is quadratic in \(f\). This nonlinearity drives the physics that matters — Landau damping's nonlinear saturation, two-stream instability, phase-space trapping — and it is exactly what a quantum circuit cannot represent directly. The motivation is identical to the seed's: kinetic plasma/fluid PDEs need enormous grids, so a future error-corrected machine might beat classical supercomputers if the nonlinearity can be tamed.
Write the discretized dynamics as a quadratic ODE \(\dot{\mathbf x} = F_1\mathbf x + F_2\,(\mathbf x\otimes\mathbf x)\). Carleman's idea: adjoin the higher monomials as new variables. Track not just \(\mathbf x\) but \(\mathbf x^{\otimes 2},\mathbf x^{\otimes 3},\dots\); each obeys a linear equation coupling only to its neighbors in degree. The toy case \(\dot x = x^2\) with \(y_r=x^r\) illustrates it exactly:
The nonlinear \(F_2\) term becomes an off-diagonal block coupling degree \(r\) to \(r{+}1\); the system is now exactly linear but infinite-dimensional. The paper generalizes this into a common framework and details three concrete embeddings: the standard Carleman (occupation-number/monomial basis), a coherent-state embedding using bosonic ladder operators \(a^\dagger,a\), and a position-space embedding with canonical \(\hat x,\hat p\) whose restriction gives Hermite expansions.
Truncate the tower at degree \(m = b + s(g-1)\) (\(b\)=degree of the output quantity, \(g\)=polynomial order, \(s\)=a convergence index). The truncated linear ODE is then solved by a standard quantum linear-system / linear-ODE routine. The qubit cost is the log of the truncated dimension \(\binom{N+m}{m}\):
logarithmic in the number of grid variables \(N\) — the exponential-compression promise, provided \(m\) stays modest.
Writing the system as \(\dot{\mathbf z}=A\mathbf z + \eta\,G[\mathbf z]\) with \(\eta\) the nonlinearity strength, they show the truncation error of an output quantity obeys
So the embedding converges — and the quantum speedup survives — only if the nonlinearity is sufficiently weak relative to the linear (dissipative/streaming) part. For Vlasov–Poisson they note a useful structural fact: under centered differences with periodic boundaries the discretized force satisfies \(\mathbf x\!\cdot\!F(\mathbf x)=0\), which keeps the position-space embedding's evolution norm-preserving. The paper is careful that this is a prospects analysis: it maps the machinery and states the condition, rather than claiming an achieved advantage.
This is the "what if we hadn't linearized?" counterpart to the seed. Miyamoto's escape hatch — replace self-gravity with an external field — made the problem exactly linear and Carleman unnecessary. Engel et al. keep the true nonlinearity and pay for it with an infinite tower and a hard convergence condition. The honest limit, flagged here and confirmed brutally by the next paper, is that the weak-nonlinearity requirement is exactly what interesting plasma physics violates: strong turbulence, large-amplitude instabilities, and phase-space structure all live at \(\eta\sim O(1)\). It is a foundational, fault-tolerant-oriented paper that defines the technique the survey needs — and sets up the sobering assessment that follows.
Vaszary et al. 2024 — nonlinear Vlasov–Krook plasma
Every survey needs one paper that says "we tried, and it doesn't work yet." Vaszary et al. take the Carleman route of the previous slide, apply it carefully to a real nonlinear Vlasov problem, and reach a negative result: no clear quantum advantage, and the convergence that Carleman needs demands unphysical dissipation. This is not a failure to include — it is the survey's honesty anchor. [arXiv:2411.19310]
They study the electrostatic Vlasov equation in one space + one velocity dimension, with a Krook collision operator \(-\nu(f-f_0)\) supplying physical relaxation — a standard, well-posed plasma test problem. The nonlinearity is the same self-consistent \(E[f]\!\cdot\!\partial_v f\) quadratic term. Crucially this is a concrete, discretized system, not a toy: the goal is to measure the actual query and gate complexity a quantum Carleman solver would incur and compare it honestly to the classical cost.
Applying the Carleman embedding produces the same block-bidiagonal linear tower as Engel et al. But when instantiated for this system the evolution matrix is large in norm, poorly sparse, and high-dimensional — the three properties that quantum linear-ODE algorithms are most sensitive to (their cost scales with sparsity, matrix norm, and condition number). The abstract virtue of "logarithmic qubits" is real, but the gate and query complexity — which depend on those matrix properties, not just the dimension — balloon.
Carleman convergence (the \(\eta^{s+1}\) bound of the previous slide) requires the linear/dissipative part to dominate the nonlinear part — formally, the negative real part of the linear operator's spectrum must exceed a threshold set by the nonlinearity. Vaszary et al. show that for typical plasma parameters the physical damping available — the Krook collision rate \(\nu\) at realistic values — is far too small to satisfy that condition:
To make the Carleman series converge one would have to add dissipation "far exceeding" what the physics provides — i.e. simulate a different, over-damped plasma. The weak-nonlinearity requirement that Engel et al. flagged is not a mild technicality here; it is violated by the very regime plasma physicists care about.
Carrying the analysis through, they find the quantum algorithm's query and gate complexities scale polynomially larger than the time complexity of the corresponding classical algorithm for this problem — the opposite of a speedup. Three compounding causes: (i) the Carleman matrix's poor sparsity and large norm inflate the Hamiltonian-/ODE- simulation cost; (ii) the truncation order needed for accuracy grows because the nonlinearity is not weak; (iii) the unphysical-dissipation fix, if adopted, changes the physics rather than the cost. This is a ~40-page, carefully worked negative result, not a hand-wave.
Read fairly, this paper sharpens rather than refutes the field. (1) It is problem-specific: it targets nonlinear electrostatic Vlasov–Krook via Carleman, and finds that particular pairing wanting — it does not prove all quantum plasma simulation is hopeless, and it says nothing against the seed paper's linear algorithm, whose external-force trick sidesteps exactly the nonlinearity that sinks Carleman here. (2) It is constructive in the negative: by pinning the failure to the dissipation-vs-nonlinearity mismatch, it tells the community precisely which knob (better-conditioned embeddings, different linearizations, or genuinely dissipative target problems) must change. (3) For the survey it is essential ballast — it prevents the "map any PDE to Carleman and claim exponential speedup" overclaim, and it motivates the alternative reformulations on this family (second quantization, spin models) that try to dodge the Carleman matrix pathology altogether. A negative result of this rigor is a contribution, not a gap.
May & Qin 2025 — nonlinear Vlasov–Poisson
If Carleman's way of linearizing nonlinear Vlasov leads to a hostile matrix (previous slide), is there a different exact linearization? May & Qin borrow the trick physics has used for a century — second quantization — to render the nonlinear Vlasov–Poisson system "linear, finite-dimensional, and discrete," at the cost of moving into a Fock space. [arXiv:2506.01895]
Same self-consistent kinetic system as the last two slides — the electrostatic Vlasov equation with its quadratic \(E[f]\!\cdot\!\partial_v f\) nonlinearity and Poisson coupling — benchmarked here on nonlinear plasma oscillations (a Bertrand-type initial condition) driven into the nonlinear regime. The ambition is to keep the nonlinear dynamics exactly, not in a weak-nonlinearity approximation, while still handing the circuit a strictly linear object.
Discretize phase space into cells \(k\) and promote the classical distribution \(f\) to occupation-number operators. Introduce ladder operators \(\hat a_k,\hat a_k^\dagger\) with \([\hat a_k,\hat a_j^\dagger]=\delta_{kj}\) and number operators \(\hat n_k=\hat a_k^\dagger\hat a_k\); the classical field is recovered as the expectation \(f_k=\langle\hat n_k\rangle\). The magic is that the quadratic nonlinearity in \(f\) becomes a bilinear operator in \(\hat a,\hat a^\dagger\) acting on Fock states — and operator evolution on Fock space is strictly linear (Heisenberg/Schrödinger dynamics of a Hamiltonian):
Nonlinearity in the physical variable \(f\) is traded for a linear equation in a higher-dimensional Fock space — the same bargain as Carleman (nonlinearity ↔ dimension), but through a many-body operator algebra rather than a monomial tower. They derive explicit conditions for correspondence between the original ("pre-quantized") and the second-quantized equations — essentially the discretization fineness and occupancy regime under which \(\langle\hat n_k\rangle\) reproduces classical \(f\).
Because the result is a strictly linear Hamiltonian system, it is amenable to quantum Hamiltonian simulation — no variational outer loop, no Carleman truncation-convergence condition — with qubits encoding the occupation numbers, \(\sim\log_2\) of the Fock dimension for the register. The paper does not compile a specific circuit; instead it establishes explicit correspondence conditions between the pre-quantized (classical VP) and second-quantized equations and supplies asymptotic scaling relations — notably a simulation cost growing as \(\sim(\dim)^4\) in the Fock dimension — arguing the resource growth is favorable relative to classical 6-D grids. The contribution is the mapping and its scaling, not a gate-level implementation.
They present numerical simulations demonstrating the quantized linear system reproduces genuinely nonlinear behavior — following nonlinear plasma oscillations from a Bertrand-type initial condition as higher-amplitude oscillations are driven into the nonlinear regime and couple to additional modes, with agreement to the classical Vlasov–Poisson evolution over the nonlinear phase (and expected divergence once fine phase-space structure outruns a coarse discretization). These are classical emulations / analyses of the mapping, not hardware runs, establishing correctness of the second-quantized correspondence and the scaling case rather than an experimental demonstration.
A third distinct answer to the nonlinearity the seed paper sidestepped. Where Miyamoto linearizes by physics (external force), Engel/Vaszary linearize by Carleman monomials, and Mocz/Cappelli keep nonlinearity in a variational loop, May & Qin linearize by second quantization — arguably the most "exact" of the reformulations and one that plugs into clean Hamiltonian simulation instead of a fragile truncation. Honest costs: (i) the Fock space is large — occupation cutoffs and the correspondence conditions bound how faithfully \(\langle\hat n_k\rangle\) tracks classical \(f\), and fidelity degrades once driven deep into the nonlinear regime; (ii) the paper stops at correspondence conditions and asymptotic scaling — no concrete circuit or resource count, so the gate-level cost of the Poisson coupling is not yet established; and (iii) readout remains a summary of \(\langle\hat n_k\rangle\), many shots per observable. It is a promising, still-theoretical direction that answers Vaszary's Carleman pathology with a different linearization — but, like everything on this family, its decisive test awaits fault-tolerant hardware.
Sundar et al. 2025 — plasma-wave propagation
Every other paper on this family is theory or emulation. Sundar et al. is the survey's one genuine hardware demonstration: plasma-wave propagation actually executed on a superconducting quantum processor. The move that makes it fit on today's noisy hardware is to give up on the full nonlinear kinetic problem and simulate the linear wave sector through a local spin model. [arXiv:2507.09479, Phys. Rev. Applied]
They simulate linear electromagnetic-wave propagation in a plasma with a spatially varying density — the model problem for a laser pulse scattering off an inhomogeneous plasma (relevant to inertial-fusion and laser–plasma experiments). The governing wave equation is
with a position-dependent plasma frequency \(\omega_p(x)\) that plays the role of the inhomogeneous medium: waves refract, reflect at sharp density jumps, and scatter off Gaussian density bumps. This is the linear face of plasma physics — deliberately not the nonlinear Vlasov of the previous three slides — which is what makes a real hardware run feasible now.
They map the wave equation onto a nearest-neighbor spin-1/2 chain whose single-excitation sector reproduces the plasma dispersion:
A single spin-wave amplitude \(b_j\) then obeys \(i\hbar\,\partial_t b_j = \tfrac{iJ}{2}(b_{j-1}-b_{j+1}) + \Delta_j(-1)^j b_j\), whose dispersion \(\hbar\omega_k=\pm\sqrt{\Delta^2+J^2\sin^2(ka)}\) matches the plasma \(\omega_k=\pm\sqrt{\omega_p^2+c^2k^2}\) at long wavelength — with the staggered field \(\Delta_j\leftrightarrow \omega_p(x)\) encoding the density inhomogeneity and the hopping \(J\leftrightarrow c\) the wave speed. Each qubit is one spatial grid point, and the mapping needs a lower gate count than alternative encodings because the coupling is strictly local (nearest-neighbor).
Hardware: a 9-qubit sublattice of Rigetti's Ankaa-3 superconducting chip (transmons, tunable-coupler \(\sqrt{i\text{SWAP}}\)-type FSIM native gates; measured error-per-layered-gate \(\approx1.5\%\)). They implement first-order Trotterized real-time evolution of \(\hat H\) with timestep \(\delta\!\approx\!0.8\), then apply an error-mitigation stack — pseudo-twirling (turning coherent FSIM errors into stochastic Pauli noise), Clifford data regression, and observable rescaling that enforces the conserved total \(\sigma^z\) — with classical-shadow readout (\(\sim\!10^6\) shots acquired in seconds).
Two experiments on the 9 qubits. (1) Dispersion: initialize one qubit in \((|0\rangle+|1\rangle)/\sqrt2\), evolve, read \(\langle\sigma^x_j\rangle+i\langle\sigma^y_j\rangle\), extract the mode frequencies — the measured spectrum agrees excellently with the exact eigenfrequencies and shows the predicted gap opening when \(\Delta\!\gt\!0\) (massive-like plasma mode) versus a gapless mode at \(\Delta\!=\!0\). (2) Wave-packet scattering: launch a two-site packet and watch it (a) propagate ballistically in "vacuum" (\(\Delta\!=\!0\)), (b) largely reflect off a sharp density jump, and (c) scatter at a Gaussian density profile — the mitigated data track the noiseless simulation qualitatively across roughly half the lattice before decoherence dominates. It is a real, if modest, demonstration of the target physics on hardware.
The reality check for the whole family. Same broad arena (kinetic/plasma waves) as the seed, but the honest scale gap is stark: 9 qubits = 9 grid points, against \(10^6\)–\(10^8\) points in production classical full-wave codes — this is a physics demonstration, explicitly not a speedup, and the authors say so. Three fair caveats: (i) it simulates the linear wave sector, deliberately avoiding the nonlinear Vlasov dynamics that defeat the Carleman papers — the ZZ terms needed for nonlinear/many-body plasma effects are flagged as future work; (ii) results survive only \(\sim\)half the lattice / limited depth before noise wins, and lean hard on error mitigation; (iii) the per-qubit-per-grid-point encoding is linear, not exponential, in spatial compression. Its value to the survey is precisely its honesty: it shows what today's superconducting hardware can genuinely do (clean linear dispersion and scattering on \(\sim\!10\) sites), and by contrast how far the fault-tolerant, exponentially-compressed algorithms of Miyamoto, Engel and May remain from execution.
Costa–Jordan–Ostrander 2019 — classical wave & EM fields
The cleanest proof-of-concept that the Vlasov trick is a general trick. Costa, Jordan & Ostrander take the classical wave equation — nothing quantum about it — and show it becomes, exactly, a Schrödinger equation that a quantum computer evolves natively. [arXiv:1711.05394, Phys. Rev. A 99, 012323 (2019)]
The scalar wave equation \(\partial_t^2\phi=c^2\nabla^2\phi\) (with Dirichlet or Neumann boundaries) governs acoustics, seismology, and electromagnetism. Classically one grids space into \(N\) points and time-steps an FDTD scheme at cost \(O(N)\) memory and \(O(N\,n_t)\) work. The paper asks: can a quantum computer evolve the field in \(\operatorname{polylog}(N)\) space and output a functional of it faster?
The wave equation is second order in time, so it is not yet Schrödinger-shaped. The move is to write the discrete Laplacian as a square of a first-order operator. On a grid graph let \(D\) be the edge–node incidence (discrete gradient) matrix; then
a Laplacian that is automatically negative semidefinite. Now define a two-component state stacking the field's time-derivative and its gradient, \(\psi=\big(\dot\phi,\ c\,D\phi\big)^{\mathsf T}\). The single second-order PDE becomes a first-order system whose generator is Hermitian and off-diagonal:
This is the same structural pivot as the Vlasov paper: rearrange a classical PDE until its generator is Hermitian, then the solution is \(\psi(t)=e^{-iHt}\psi(0)\). The authors extend it verbatim to the Klein–Gordon equation (add a mass term) and to Maxwell's equations, which in vacuum are already first order — the Riemann–Silberstein vector \(\mathbf F=\mathbf E+i\mathbf B\) obeys \(i\partial_t\mathbf F=c\,\nabla\times\mathbf F\), a Schrödinger equation for a "photon wavefunction."
Because \(H\) is \(s\)-sparse (a stencil couples \(O(1)\) neighbours) it is Hamiltonian-simulable. Using the Laplacian factorization to sharpen truncation error and improve state preparation, the gate cost is
i.e. near-linear in the evolution time \(t\), only \(\log\) in the grid size \(N\), and \(\log(1/\epsilon)\) in accuracy — an exponential compression of the \(N\)-dependence versus classical FDTD. The improvement over a generic linear-ODE solver comes precisely from exploiting \(\nabla^2=-D^{\mathsf T}D\).
- Read-out. You do not get the \(2^n\) field values; you get an expectation \(\langle\psi|M|\psi\rangle\) (an energy in a region, a detector overlap). The full waveform is inaccessible.
- State preparation. The initial field/source must be loadable in \(\operatorname{polylog}\) depth — else loading alone costs \(O(N)\) and erases the win.
- Fault tolerance. The block-encoding + Hamiltonian-simulation circuits assume error correction, not NISQ.
This is the same engine, different physics. Miyamoto's Vlasov reduction (real antisymmetric \(A\Rightarrow H=iA\) Hermitian \(\Rightarrow e^{-iHt}\)) is one instance of the paradigm Costa et al. crystallize: a linear PDE is secretly a Schrödinger equation. If you understand why the wave equation maps to \(e^{-iHt}\), you understand why the collisionless Vlasov equation does too — and that the entire class (wave, Klein–Gordon, Maxwell, Vlasov) shares one quantum backbone: sparse block-encoding + optimal Hamiltonian simulation, with read-out and state-prep as the shared soft spots.
Gaitan 2020 — compressible Navier–Stokes flow
A different route to a fluid: don't turn the PDE into a Schrödinger equation — turn it into a big system of ODEs and solve the ODEs with a quantum integrator whose engine is amplitude estimation. Gaitan's is the first end-to-end quantum CFD demonstration. [npj Quantum Information 6, 61 (2020)]
Solve the compressible Navier–Stokes / Euler equations for a real engineering flow: steady, inviscid, compressible air through a convergent–divergent (de Laval) nozzle — the geometry inside rocket and supersonic jet engines — and recover the full profiles of Mach number, pressure, density and temperature, including shock capture when the back-pressure forces a standing shock in the diverging section.
Discretise the spatial derivatives (method of lines) so the PDE becomes an initial-value problem for a state vector \(\mathbf y\),
integrated forward to steady state. The IVP is solved with Kacewicz's quantum ODE algorithm, an adaptive Taylor/Runge–Kutta scheme in which the one expensive ingredient — the mean of the vector field over each subinterval, \(\bar{\mathbf f}=\tfrac1h\int_t^{t+h}\mathbf f\,dt'\) — is evaluated by quantum amplitude estimation rather than classical Monte Carlo. Amplitude estimation returns a mean to additive error \(\epsilon\) with
a quadratic speedup in the integration precision that propagates to the whole ODE solve. The nonlinearity of Navier–Stokes is handled by the classical outer Taylor/RK structure — the quantum part only estimates integrals.
For a vector field with \(r\) continuous derivatives, Kacewicz's algorithm attains an \(\epsilon\)-accurate solution with quantum query cost scaling like \((1/\epsilon)^{1/(q+1)}\) — where \(q=r+\rho\) is the Hölder smoothness of the vector field — versus the best classical randomized bound \((1/\epsilon)^{1/(q+1/2)}\). The gain is polynomial in the precision \(\epsilon\); crucially, only in the rough-solution limit (\(q,\gamma\ll1\)) does it reach the headline \(O(1/\epsilon)\)-vs-\(O(1/\epsilon^2)\) advantage of amplitude estimation — for smooth fields (\(q\gg1\)) the quantum and classical exponents coincide and there is no quantum advantage. Gaitan's simulation matched the exact nozzle solution, shock included.
- Not an exponential speedup. The quantum win is only in the integration subroutine's \(\epsilon\) scaling; the state is read out and re-loaded each step, so there is no exponential compression of the grid as in the Schrödinger route.
- Low-dimensional. The demonstrated system is a small ODE (quasi-1D). The method does not by itself defeat the curse of dimensionality that motivates the Vlasov work.
- Oracle + fault tolerance. Requires a coherent oracle for \(\mathbf f\) and deep phase-estimation circuits — fault-tolerant, not NISQ.
Two shared tools. First, amplitude estimation is the exact primitive Miyamoto uses to read the neutrino power spectrum \(P_\nu(k)\) to accuracy \(\epsilon\) in \(O(1/\epsilon)\) — Gaitan puts the same primitive at the heart of the solver instead of the read-out. Second, it illustrates the rival reduction: "PDE → ODE → quantum ODE solver," a genuinely nonlinear-capable path that trades the exponential grid win for the ability to keep the nonlinearity. Reading Gaitan next to the seed makes explicit what Miyamoto buys by first linearizing to a Schrödinger form — and what it gives up (self-consistent nonlinearity).
Budinski 2021 — incompressible Navier–Stokes flow
The kinetic cousin of the Vlasov route. Lattice Boltzmann already describes a fluid by a distribution function over discrete velocities; Budinski maps its stream–collide cycle onto quantum operations, giving the first quantum LBM solver for Navier–Stokes. [arXiv:2103.03804, Int. J. Quantum Information 19, 2150039 (2021)]
Solve incompressible Navier–Stokes in the stream-function / vorticity formulation (which recasts the dynamics as coupled advection–diffusion equations), demonstrated on 2-D lid-driven cavity and Taylor–Green flows on IBM's Qiskit simulator. LBM is attractive because its update is local and explicit — a good fit for a circuit.
Lattice Boltzmann evolves populations \(f_i(\mathbf x,t)\) along a small set of discrete lattice velocities \(\mathbf c_i\) via the BGK stream–collide rule
Budinski amplitude-encodes the populations into a quantum state and implements the two half-steps as quantum primitives:
- Collision (the relaxation toward \(f^{\rm eq}\), local at each site) via a standard-form / linear combination of unitaries (LCU) block-encoding of the collision operator.
- Streaming (the shift \(\mathbf x\to\mathbf x+\mathbf c_i\Delta t\), conditioned on the velocity register) via a quantum-walk shift operator — a controlled permutation of grid sites.
The advection–diffusion (stream-function/vorticity) form keeps the mapped operator effectively linear, so the quadratic moment-coupling of full LBM is sidestepped in this demonstration.
This is a proof of concept, not a proven speedup. The critical obstruction is that collision is non-unitary (relaxation dissipates the \(L^2\) norm), so it must be realized by LCU with post-selection, whose success probability decays and which forces measurement and re-preparation of the state each time step. That per-step read-out/re-load is the dominant cost and no end-to-end exponential advantage is claimed.
- Non-unitarity of collision ⇒ post-selection, shrinking success amplitude, repeated re-initialisation.
- Nonlinearity not fully captured — only the linearized advection–diffusion form is mapped cleanly.
- Macroscopic read-out (velocity, vorticity fields = moments of \(f_i\)) is itself expensive.
LBM and Vlasov are both discrete-velocity kinetic schemes: they evolve a distribution over \((\mathbf x,\mathbf v)\), and their streaming operator is the very same "shift along velocity" structure as Miyamoto's block \(A_x=-D_{\rm per}\otimes\dots\otimes E_u\). The instructive contrast is the collision step: Budinski's dissipative, non-unitary collision breaks the norm-preserving structure and forces per-step measurement — exactly the property Miyamoto engineers away by choosing central differences so that \(A\) is antisymmetric and \(H=iA\) generates a unitary, measurement-free evolution. Seeing LBM's struggle is the sharpest way to appreciate why the Vlasov paper fights so hard for antisymmetry.
Liu–Childs et al. 2021 — nonlinear dissipative ODEs
The landmark theory paper for the nonlinear branch: the first quantum algorithm with a provable efficiency guarantee for nonlinear differential equations — and a sharp statement of exactly when it works. This is the theoretical backbone every nonlinear quantum-PDE claim is measured against. [arXiv:2011.03185, PNAS 118, e2026805118 (2021)]
Solve a quadratically-nonlinear ODE system (the generic form after discretising a nonlinear PDE — Burgers, reaction–diffusion, Navier–Stokes),
where \(F_1\) carries linear dissipation (eigenvalues with negative real part) and \(F_2\) the nonlinearity. Quantum mechanics is linear, so nonlinearity is the fundamental obstruction — this paper is how you get around it with a guarantee.
Carleman linearization embeds the nonlinear system into an infinite-dimensional linear one by tracking the tower of tensor powers \(\mathbf u,\ \mathbf u^{\otimes2},\ \mathbf u^{\otimes3},\dots\). Because \(\tfrac{d}{dt}\mathbf u^{\otimes k}\) couples only to levels \(k\!-\!1,\,k,\,k\!+\!1\), the lifted system is an infinite block-tridiagonal linear ODE:
Truncate at level \(N\) and solve the resulting finite linear ODE with the forward-Euler + quantum linear-system approach. Convergence is governed by one dimensionless number — the ratio of nonlinearity (and forcing) to dissipation,
Their new convergence theorem shows the truncation error decays like \(R^{N}\), so for \(R\lt1\) a truncation \(N\sim\log(1/\epsilon)\) already suffices.
For \(R\lt1\), the algorithm produces a state \(\propto\mathbf u(T)\) with gate/query complexity
an exponential improvement over prior approaches whose cost grew exponentially in the evolution time \(T\).
- \(R\lt1\) is restrictive: dissipation must dominate nonlinearity. Turbulence (high Reynolds number) sits at \(R\gt1\) and is not covered.
- Proven hardness: for \(R\ge\sqrt2\) they show the problem is intractable — strongly nonlinear dynamics would enable state discrimination/cloning, contradicting quantum mechanics. This is a hard wall, not a gap in the proof.
- Inherits QLSA caveats: condition number, state preparation, and the read-out-a-summary-only limit all carry over.
Miyamoto sidesteps nonlinearity physically: because neutrinos are \(\lt1\%\) of the matter, the CDM force is external and the Vlasov equation is exactly linear — so no Carleman, no \(R\), no truncation error. Liu et al. is the map of the territory the seed carefully avoids: the moment a Vlasov problem becomes self-consistent (plasma Vlasov–Poisson, or neutrinos sourcing their own gravity), the nonlinearity returns and Carleman with the \(R\lt1\) condition becomes the only route to a guarantee — one that provably fails in the turbulent/strongly-clustered regime. It defines both the reach and the ceiling of the whole nonlinear quantum-PDE program.
Martinez et al. 2016 — Schwinger-model pair creation
The first real-time quantum simulation of a lattice gauge theory on a digital quantum computer. Martinez, Muschik, Schindler, Nigg, Erhard, Heyl, Hauke, Dalmonte, Monz, Zoller & Blatt watched electron–positron pairs boil out of the vacuum of 1+1-dimensional quantum electrodynamics — the Schwinger model — on a four-ion trap, turning an abstract QFT process into a measured, coherent circuit. [arXiv:1605.04570 · Nature 534, 516 (2016)]
The Schwinger model is QED in one space and one time dimension: charged fermions coupled to a U(1) gauge field. Its defining feature is that the vacuum is unstable — an electric field can spontaneously create particle–antiparticle pairs, the field-theory analogue of a spark jumping a gap. The dynamics are governed by the Hamiltonian density combining a fermion mass term, a gauge-field energy \(\propto E^2\), and a minimal coupling. Because the coupling is strong and the process is genuinely non-perturbative, it is exactly the kind of real-time evolution that classical lattice methods, built for imaginary time, struggle to reach.
The continuous gauge field is eliminated: in 1+1D Gauss's law lets one solve for the electric field in terms of the fermion charges, so the theory reduces to interacting fermions alone. A Jordan–Wigner transformation then turns those lattice fermions into spins, and the eliminated gauge field reappears as long-range spin–spin interactions — the Coulomb interaction of 1D. The result is a spin Hamiltonian \(H=\sum_{i}\big(\text{mass}\big)\sigma^z_i+\text{hopping}\big(\sigma^+_i\sigma^-_{i+1}+\text{h.c.}\big)+\text{long-range }\sigma^z\sigma^z\) that a trapped-ion machine implements natively, since ions couple through phonon-mediated long-range forces.
A trapped-ion processor at Innsbruck: four \(^{40}\)Ca\(^+\) ions encoding four qubits, representing two matter sites and two antimatter sites of the lattice. Single- and two-qubit entangling gates were driven by laser pulses using the Mølmer–Sørensen interaction, and the circuit ran the digital (Trotterized) time evolution of the spin Hamiltonian step by step.
Starting from the bare vacuum (all matter sites empty), they Trotter-evolved and tracked the particle-number density as a function of time, watching it rise as pairs were produced and then oscillate — the tell-tale signature of coherent pair creation and recombination. They mapped the production as the fermion mass was tuned, and separately followed the growth of entanglement between the produced particles, demonstrating that the created pairs are quantum-correlated, not classical noise. The runs used a handful of Trotter steps before gate errors accumulated.
This was the proof of concept that a full gauge theory — not just a spin toy — could be encoded and time-evolved on a programmable digital quantum device, with gauge invariance built in by construction through the Gauss-law elimination. Honest limits: only four qubits (two physical sites), a strictly 1+1D Abelian theory, and a shallow circuit whose depth was capped by decoherence, so no scaling advantage was claimed. Its importance is foundational — it opened the entire experimental program of digital LGT simulation that Klco, Nguyen and others extended.
Klco et al. 2018 — Schwinger dynamics on cloud qubits
The first Schwinger-model dynamics run on a cloud-accessible superconducting quantum computer, showing that publicly available hardware — not just bespoke lab machines — could reach into quantum field theory. Klco, Dumitrescu, McCaskey, Morris, Pooser, Sanz, Solano, Lougovski & Savage combined classical symmetry analysis with short quantum circuits to simulate a two-site Schwinger model on IBM's processors. [arXiv:1803.03326 · Phys. Rev. A 98, 032331 (2018)]
The target is again the Schwinger model — 1+1D quantum electrodynamics — and specifically its real-time dynamics: how an initial state evolves under the interplay of the fermion mass, the electric-field energy, and the hopping that creates and annihilates charges. The physics of interest is the coherent oscillation of the field between vacuum-like and particle-populated configurations, a clean test bed because the two-site model is small enough to know the exact answer yet rich enough to exhibit genuine gauge dynamics.
As in the trapped-ion case, Gauss's law removes the gauge field and a Jordan–Wigner map turns the staggered lattice fermions into qubits. The authors' distinctive move is to first exploit the model's symmetries — total charge, parity, and rotational structure — classically, to identify the small physical subspace the dynamics actually explore. This symmetry reduction shrinks the qubit requirement by roughly a factor of five, so that only the quantum fluctuations within a fixed symmetry sector need to be placed on hardware while the sector labels are fixed by hand.
Superconducting transmon qubits on IBM's cloud machines (the publicly accessible QX processors), using only a small number of qubits after the symmetry reduction — a two-spatial-site Schwinger model. Circuits were compiled to the device's native gate set and submitted remotely, making this an early demonstration of QFT simulation as a cloud workload.
They Trotterized the time evolution within the chosen symmetry sector and measured the resulting real-time dynamics of the two-site system, extracting observables such as the persistence of the initial state and its oscillation in time. Because the reduced problem was tiny, the measured trajectories could be compared directly against exact classical results, and the paper is candid about the role of error mitigation: raw device output was noisy, and post-processing was needed to recover physically meaningful curves.
The contribution is methodological: symmetry-informed qubit reduction plus error mitigation lets even noisy, remotely accessed superconducting hardware reproduce gauge-theory dynamics. Honest limits: a two-site model, heavy classical pre-processing that does much of the work, and results dominated by mitigation rather than raw fidelity — no advantage over classical computation. Its place in the landscape is as the bridge from Martinez's dedicated ion trap to the era of democratized, cloud-based lattice-gauge experiments.
Yang et al. 2020 — gauge invariance in a 71-site simulator
A large-scale analog cold-atom realization of a U(1) lattice gauge theory in which Gauss's law itself — the local gauge constraint that defines a gauge theory — was experimentally quantified. Yang, Sun, Ott, Wang, Zache, Halimeh, Yuan, Hauke & Pan engineered a 71-site Bose–Hubbard superlattice so that its natural dynamics enforce local gauge symmetry between matter and gauge fields. [arXiv:2003.08945 · Nature 587, 392 (2020)]
A gauge theory is not defined by its Hamiltonian alone but by a local constraint at every site: Gauss's law, \(\nabla\!\cdot\!\mathbf E=\rho\), tying the electric field's divergence to the local charge. In 1+1D lattice QED this means each link's gauge field is locked to the matter occupation on its neighbouring sites. The deep experimental challenge is that nothing forces a generic quantum simulator to respect this constraint — gauge invariance must be built into the microscopic dynamics, and then demonstrated rather than assumed.
The trick is an analog encoding: bosonic atoms hopping in an optical superlattice are arranged so that matter sites and gauge-field links alternate. By tuning the lattice so that the atomic interactions and tilts reproduce the required correlated hopping, a matter atom can only move if the adjacent gauge link changes accordingly — the Bose–Hubbard dynamics are mapped onto the U(1) gauge-theory dynamics, with Gauss's law emerging from energetic constraints rather than being imposed by hand. There is no gate circuit here; the Hamiltonian is engineered directly.
An ultracold \(^{87}\)Rb Bose–Einstein condensate loaded into a tilted optical superlattice forming a 71-site chain (thousands of atoms across the array), at USTC. This is an analog quantum simulator: the physics lives in the continuously evolving many-body atomic state, read out by site-resolved occupation measurements, and its scale dwarfs the handful of qubits available to digital devices of the era.
They prepared the system and let it evolve, then extracted the probabilities of gauge-invariant atomic configurations from correlated occupation measurements, thereby directly quantifying how much Gauss's law was violated — and showing the violation was small, i.e. the dynamics stayed within the physical gauge sector. Tuning the coupling, they also observed a quantum phase transition in the gauge model, confirming that the emulator reproduced non-trivial gauge dynamics across its full parameter range.
This established that gauge invariance can be experimentally certified at a scale (71 sites) far beyond digital circuits, proving analog platforms viable for large lattice gauge theories. Honest limits: it is analog, so gauge symmetry is approximate (enforced energetically, with residual violation) rather than exact; the theory is 1+1D and Abelian U(1); and the device is special-purpose, not programmable. It anchors the analog branch of the field — the large-scale counterpoint to the small-but-digital line of Martinez and Nguyen.
Kokail et al. 2019 — self-verifying variational Schwinger model
A variational quantum simulation of the lattice Schwinger model that computes its own error bars — the first "self-verifying" gauge-theory simulation. Kokail, Maier, van Bijnen, Brydges, Joshi, Jurcevic, Muschik, Silvi, Blatt, Roos & Zoller used up to twenty trapped ions to variationally find the ground state of 1+1D QED and, crucially, to certify how close they got. [arXiv:1810.03421 · Nature 569, 355 (2019)]
Rather than real-time dynamics, the target here is the ground state of the lattice Schwinger Hamiltonian — the interacting vacuum of 1+1D QED, dressed with virtual pairs. Knowing this state gives access to static properties like the chiral condensate and the mass gap, quantities that in real QCD require enormous classical Monte Carlo effort. The physics challenge is that the true vacuum is a highly entangled superposition, and finding it means minimizing energy over an exponentially large space.
Gauss's law again eliminates the gauge field and Jordan–Wigner maps the fermions to spins. The state is not built by a fixed Trotter circuit but by a parametrized variational ansatz: a sequence of tunable entangling operations whose angles \(\boldsymbol\theta\) define a trial state \(|\psi(\boldsymbol\theta)\rangle\). A classical optimizer adjusts \(\boldsymbol\theta\) to minimize the measured energy \(\langle\psi|H|\psi\rangle\) — the variational quantum eigensolver strategy — with the ansatz designed to respect the model's symmetries so the search stays in the physical sector.
A programmable trapped-ion analog quantum simulator at Innsbruck using up to 20 \(^{40}\)Ca\(^+\) ions (20 qubits) — a large register for 2019. Global laser-driven entangling operations generated the variational ansatz, and individual-ion readout supplied the energy and its variance to the classical optimization loop.
They ran the closed variational loop to converge on the Schwinger-model ground state, then did the paper's signature step: they measured the variance of the Hamiltonian, \(\langle H^2\rangle-\langle H\rangle^2\), on the prepared state. That variance bounds how far the trial state is from a true eigenstate, yielding a rigorous, device-independent error bar on the energy — the simulation verifies itself without needing a classical benchmark of the answer. This was demonstrated on systems up to the full 20-qubit size.
The breakthrough is verification: it answers the deepest worry about quantum simulation — how do you trust an answer no classical computer can check — by having the device certify its own accuracy from measured variances. Honest limits: variational optimization can stall in local minima and the number of measurements needed grows quickly; it is still 1+1D Abelian QED at modest size; and the reachable precision is set by gate noise. It remains a landmark for the self-certification paradigm now sought across the whole NISQ program.
Jordan-Lee-Preskill 2012 — scattering in scalar QFT
The theoretical foundation of the whole field: the first quantum algorithm proven to compute relativistic scattering amplitudes in an interacting quantum field theory with a provable speedup. Jordan, Lee & Preskill gave an explicit, efficient procedure for simulating \(\phi^4\) scalar field theory end to end — prepare, scatter, measure — establishing that quantum computers can attack QFT in a regime where classical methods provably fail. This is a proposal: an algorithm with complexity proofs, not a hardware experiment. [arXiv:1111.3633 · Science 336, 1130 (2012)]
The target is a scattering experiment inside a scalar field theory with a quartic self-interaction, \(\mathcal L=\tfrac12(\partial\phi)^2-\tfrac12 m^2\phi^2-\tfrac{\lambda}{4!}\phi^4\), in up to four spacetime dimensions. The physics to be computed is the S-matrix: send in wavepackets, let them interact, and read off the probability of the outgoing states. At strong coupling this is genuinely hard — perturbation theory diverges and the fastest known classical algorithms scale exponentially — so it is a natural proving ground for quantum advantage in fundamental physics.
The continuous field is discretized on a spatial lattice, and the field value at each site — a continuous variable — is represented in a truncated basis on a register of qubits. The algorithm then has three stages: (i) prepare the interacting vacuum and inbound wavepackets by starting from the free theory and adiabatically turning on the coupling \(\lambda\); (ii) simulate the time evolution under the full Hamiltonian by Trotterization; and (iii) measure the outgoing particle content. Careful treatment of the adiabatic state preparation is what makes the whole chain efficient.
There is no device — this is an algorithm and a complexity theorem. The key result is that the number of qubits and gates needed scales polynomially in the number of scattered particles, their energy, and the target precision, at both weak and strong coupling. Against the best known classical algorithms, this represents an exponential speedup in the strong-coupling and high-precision regime. The paper also carefully analyzes the error sources — lattice spacing, field truncation, adiabatic ramp speed, Trotter step — and bounds each.
They delivered a complete, rigorously costed recipe for a canonical QFT computation, proving that simulating quantum field theory is in the class of problems a quantum computer solves efficiently (BQP) even where classical field theory is intractable. The demonstration is entirely mathematical: convergence and resource proofs rather than measured data, with the \(\phi^4\) model chosen as the simplest interacting relativistic theory that already exhibits the essential difficulty.
This is the intellectual origin point that legitimized quantum simulation of high-energy physics: it turned "could a quantum computer do particle physics?" into a proven yes with an explicit algorithm. Honest limits: it targets a scalar theory (no gauge fields, no fermions in the original construction), assumes a fault-tolerant machine far beyond any existing hardware, and the polynomial has large constants — it is an existence proof of advantage, not a near-term protocol. Every later experimental paper in this family is, in spirit, chasing the vision this proposal set out.
Nguyen et al. 2022 — symmetry-protected Schwinger dynamics
A digital trapped-ion simulation of the Schwinger model that used the theory's own gauge symmetry as an error-mitigation tool, reaching real-time dynamics far longer than previously possible. Nguyen, Tran, Zhu, Green, Huerta Alderete, Davoudi & Linke turned Gauss's law from a constraint to be respected into an active shield against hardware noise. This is the Abelian U(1) Schwinger model. [arXiv:2112.14262 · PRX Quantum 3, 020324 (2022)]
The target is real-time dynamics of the 1+1D lattice Schwinger model, and specifically non-perturbative pair production — the same vacuum instability Martinez first saw, but followed to much later times. Long-time evolution is where the interesting non-equilibrium physics lives (thermalization, string breaking, repeated pair creation/annihilation), and it is exactly where accumulated gate errors normally destroy the signal, so reaching it is both a physics and a hardware challenge.
Gauss's law eliminates the gauge field, Jordan–Wigner maps the fermions to qubits, and the dynamics are Trotterized into a gate circuit. The novelty is in handling errors. The physical Hilbert space is only the gauge-invariant sector — states satisfying Gauss's law — and coherent hardware errors tend to push the state out of that sector. Two ideas exploit this: a symmetry-protection protocol that adds control terms energetically penalizing gauge-violating leakage, suppressing coherent errors during evolution; and a symmetry-based post-selection that discards any measured shot violating the conserved symmetries, cleaning the final data.
A trapped-ion digital quantum computer (the University of Maryland \(^{171}\)Yb\(^+\) system) using up to six qubits, with all-to-all connectivity via Mølmer–Sørensen gates. Full programmability let the authors insert the symmetry-protection control terms directly into the compiled circuit.
They ran Trotterized evolution and measured the particle-number density (pair-creation signal) as a function of time, comparing runs with and without symmetry protection. With the protection and post-selection applied, the simulation faithfully tracked non-perturbative pair creation for substantially longer evolution times / more Trotter steps than unprotected runs, where noise otherwise washed the oscillations out — a concrete, measured demonstration that gauge symmetry can extend usable circuit depth.
It demonstrated a broadly reusable principle: the conservation laws of the target theory are a free, physics-based resource for error mitigation on NISQ hardware. Honest limits: still six qubits and 1+1D Abelian U(1); post-selection throws away shots and so costs statistics; and symmetry protection suppresses only the coherent, symmetry-violating part of the error budget, not incoherent decoherence. It marks the maturing of digital LGT simulation from "can we run it at all" (Martinez) toward "can we run it long enough to matter."
Banuls et al. 2020 — review of lattice gauge simulation
The definitive community review of simulating lattice gauge theories with quantum technologies — the map of the whole territory these experimental papers occupy. Bañuls, Blatt, Catani, Celi, Cirac, Dalmonte, Fallani, Jansen, Lewenstein, Montangero, Muschik, Reznik, Rico, Tagliacozzo, Van Acoleyen, Verstraete, Wiese, Wingate, Zakrzewski & Zoller synthesized the theoretical proposals, tensor-network results, and hardware demonstrations into one coherent framework. This is a review, not a single experiment. [arXiv:1911.00003 · Eur. Phys. J. D 74, 165 (2020)]
The review's subject is the simulation of gauge theories themselves — from the Abelian Schwinger model up toward non-Abelian theories relevant to Quantum Chromodynamics. It frames the central difficulty: real-time dynamics, finite fermion density, and topological terms all induce a sign problem that cripples classical Monte Carlo, motivating quantum and tensor-network approaches. It also stresses the organizing principle every method must honour — local gauge invariance / Gauss's law — and the choices for representing continuous gauge fields on finite quantum registers.
The review catalogues the mappings the field uses: eliminating gauge fields via Gauss's law, Jordan–Wigner and Kogut–Susskind formulations, quantum-link / gauge-magnet models that give gauge fields a finite-dimensional Hilbert space suited to atoms and spins, and truncation schemes for the gauge degrees of freedom. It contrasts how each maps onto qubits and gates (digital) versus engineered Hamiltonians (analog), and how gauge invariance is either imposed exactly or protected energetically.
Being a review, it has no device of its own; instead it surveys three complementary pillars: (i) tensor-network methods — a classical, sign-problem-free approach for Abelian and non-Abelian theories; (ii) analog quantum simulators in cold atoms, trapped ions, Rydberg arrays and superconducting circuits; and (iii) digital quantum simulators with platform-specific gate proposals. It spans roughly 45 pages and places the era's landmark experiments — Martinez's ion trap, cold-atom gauge realizations — within this taxonomy.
The review's contribution is synthesis rather than data: it collects the theoretical toolkit, the state of the art in each platform circa 2019–2020, the resource requirements and error challenges, and — most usefully — the open problems: reaching higher dimensions, implementing genuine non-Abelian gauge fields, controlling gauge-symmetry breaking, and scaling beyond a few sites. It thereby defines the roadmap the subsequent hardware papers set out to follow.
It is the field's shared reference frame, written by the very groups building the experiments, and it clarified that no single platform wins — digital, analog and tensor-network methods are complementary tools for different regimes. Honest limits: as a 2020 snapshot it predates the largest recent demonstrations, and it charts a path rather than proving new physics. Its value is orientation: it is the document that lets any of the six experimental papers above be understood as one move in a common, well-mapped program toward simulating the Standard Model.
Bernien et al. 2017 — Rydberg quantum magnetism
The paper that opened programmable neutral-atom quantum simulation at scale: a chain of 51 individually trapped atoms tuned into a strongly interacting Ising magnet, used to watch symmetry-breaking order form and to uncover anomalously long-lived coherent revivals now understood as quantum many-body scars. Bernien et al., “Probing many-body dynamics on a 51-atom quantum simulator,” Nature 551, 579 (2017) [arXiv:1707.04344 · Nature 551, 579]. It set the template every later Rydberg-array result on this map builds on.
Each atom is a two-level pseudospin — electronic ground state \(|g\rangle\) versus a highly excited Rydberg state \(|r\rangle\) — so the array realises a transverse-field Ising-type spin model,
where \(n_i=|r\rangle\langle r|_i\), \(\Omega\) is the laser Rabi drive (the transverse field), and \(\Delta\) the detuning (a longitudinal field). The van-der-Waals tail \(C_6/r^6\) is the interaction. Sweeping \(\Delta\) from negative to positive drives the system from a trivial paramagnet into spatially ordered antiferromagnetic phases that break a discrete translational symmetry — genuine quantum many-body physics, not single-particle dynamics.
The mechanism that makes this Ising physics is the Rydberg blockade: two atoms closer than the blockade radius \(R_b\), defined by \(C_6/R_b^{6}=\hbar\Omega\), cannot both be excited because the interaction shift pushes the doubly-excited state out of resonance. Within \(R_b\) the constraint \(n_i n_{i+1}\!\to\!0\) is effectively hard — neighbouring “spins” cannot both point up — which is exactly the strong antiferromagnetic coupling of an Ising chain. By choosing the lattice spacing relative to \(R_b\), the blockade can be made to span one, two, or three neighbours, programming \(\mathbb{Z}_2\), \(\mathbb{Z}_3\), \(\mathbb{Z}_4\) ordered ground states respectively.
The device is an analog simulator: up to 51 individual \(^{87}\)Rb atoms held in a one-dimensional array of optical tweezers, rearranged atom-by-atom into a defect-free chain, then coherently driven to the \(70S\) Rydberg level. Interactions and fields are set by laser parameters and geometry rather than by gates — the whole Hamiltonian above is engineered directly and evolved continuously in time, with single-site readout by fluorescence imaging. Fifty-one coherent, strongly interacting spins was, at the time, a record-scale controllable quantum system.
Two headline experiments. First, adiabatic sweeps of \(\Delta\) prepared and verified the \(\mathbb{Z}_2\) antiferromagnet with high per-site fidelity, mapping the quantum phase transition into the ordered phase. Second, a sudden quench from the \(\mathbb{Z}_2\)-ordered product state \(|{\cdots}grgrgr{\cdots}\rangle\) into the resonant point \(\Delta\!=\!0\) produced persistent, slowly-decaying oscillations of the local Rydberg density — coherent revivals that survived far longer than the thermalization time naive expectations allow. Generic many-body systems should scramble such a state promptly; this one refused to.
This established that tweezer arrays are a bona fide programmable quantum simulator: system size, geometry, and interaction range are all tunable, and states are prepared and read out site-resolved. The unexplained revivals seeded an entire subfield — they were soon identified as quantum many-body scars, a weak, non-thermal subspace embedded in an otherwise chaotic spectrum, a new form of ergodicity breaking distinct from integrability or localization. Honest limits: it is a 1D, analog, van-der-Waals Ising model, and “scars” is the interpretation attached afterward, not proven in the paper itself. But as the origin point for Rydberg quantum magnetism its influence is hard to overstate. arXiv:1707.04344
Arute et al. 2020 — Fermi-Hubbard charge/spin dynamics
A digital, gate-based quantum computer used to simulate the real-time dynamics of the one-dimensional Fermi-Hubbard model and to resolve one of its signature strong-correlation effects: charge–spin separation, where a local disturbance splits into charge and spin excitations that travel at different speeds. Arute et al. (Google Quantum AI), “Observation of separated dynamics of charge and spin in the Fermi-Hubbard model” [arXiv:2010.07965 · Google Quantum AI preprint]. It is a benchmark for how far digital fermionic simulation had come on NISQ hardware.
The Fermi-Hubbard model is the minimal description of interacting electrons on a lattice,
with hopping \(J\), on-site repulsion \(U\), and \(\sigma=\uparrow,\downarrow\). In one dimension a fermion is not a stable quasiparticle: an added electron fractionalises, and its charge and spin degrees of freedom propagate as independent collective excitations with distinct velocities — the hallmark of Tomonaga–Luttinger physics, well beyond the single-quasiparticle picture. Capturing this requires resolving genuinely correlated, time-dependent many-fermion dynamics.
This is a digital simulation. The fermionic operators are mapped to qubit Pauli operators by a Jordan–Wigner (fermionic-swap) encoding, so each lattice site×spin becomes a qubit. The evolution \(e^{-iHt}\) is Trotterised into a repeated circuit of two-qubit gates implementing the hopping and interaction terms. The decisive engineering contribution was a fast, accurate gate-calibration procedure — fast enough to track drift of gate parameters between runs — combined with layered error-mitigation (post-selection on particle number, echo/rescaling) to claw back coherence lost to gate noise.
The processor is Google’s Sycamore superconducting-transmon chip, of which 16 qubits were used — an 8-site Hubbard chain with two spin species. Compared with the analog Rydberg arrays elsewhere on this map, the qubit count is small, but each qubit is fully gate-programmable, so the same hardware can in principle simulate any Hamiltonian expressible in these gates, not just one engineered interaction.
Starting from an initial state carrying localized charge and spin perturbations, they ran the Trotterised Hubbard evolution and tracked the spreading of the charge density and the spin density separately as functions of site and time. In the strongly excited, highly correlated regime they measured the two profiles spreading at different velocities — a direct real-time observation of charge–spin separation — and demonstrated that the calibration-plus-mitigation stack was essential to keep the signal above the noise floor across the required circuit depth.
It showed a digital quantum computer resolving a genuinely non-trivial strong-correlation phenomenon rather than a toy instance, and it made a lasting methodological point: on NISQ hardware, careful calibration and error mitigation can be as decisive as qubit count. Honest limits: only 16 qubits and an 8-site chain, results lean heavily on mitigation, the circuits are shallow-depth Trotter approximations, and — flagged on this map — the work remained a preprint rather than the Science paper it is sometimes cited as. It sits as a proof of principle for digital fermionic dynamics, complementary to the analog spin simulators. arXiv:2010.07965
Mi et al. 2022 — discrete time-crystal order
A superconducting processor used to realise and certify a genuine discrete time crystal — a phase of matter whose defining property is that it breaks the time-translation symmetry of a periodic drive, ticking at twice the drive period, robustly and indefinitely. Crucially the authors demonstrated eigenstate order, establishing the effect as a true phase of the many-body spectrum rather than a fragile transient. Mi et al. (Google Quantum AI), “Time-crystalline eigenstate order on a quantum processor,” Nature 601, 531 (2022) [arXiv:2107.13571 · Nature 601, 531].
A discrete time crystal (DTC) is a periodically driven (Floquet) many-body system whose observables respond at a subharmonic of the drive — here period-doubling — that is rigid against changes in drive parameters and initial states. For this to be a real phase and not a slow transient, the entire many-body spectrum must be many-body localized (MBL) so that the system never heats to a featureless infinite-temperature state. The test of “eigenstate order” is that the subharmonic persistence must hold from generic initial states across the spectrum, not just from one fine-tuned state.
This is a digital Floquet simulation. Each period applies a two-layer unitary,
built from calibrated single- and two-qubit gates. The transverse rotation angle \(g\) is deliberately set near but not exactly \(\pi\) (a spin flip), so that time-crystalline rigidity — not fine-tuning — must be what holds the period-doubling in place. Disorder in the couplings \(J_i\) and fields \(h_i\) supplies the MBL that prevents heating. Running the same circuit from many random initial bitstrings probes the whole spectrum.
The experiment used a 20-qubit chain on Google’s Sycamore superconducting processor. The gate-based platform is essential here: the DTC requires programmable disorder, tunable interaction strength, and an adjustable drive angle — parameters that are set gate-by-gate — plus the ability to average over hundreds of distinct random initial states, which a fixed analog Hamiltonian could not sweep as flexibly.
They evolved the 20-qubit chain for many Floquet periods and measured the site-resolved polarization \(\langle\sigma^z_i\rangle\), observing stable period-2 oscillations whose amplitude did not decay over the accessible drive cycles. To prove eigenstate order rather than a lucky initial state, they averaged over a large ensemble of random initial bitstrings and confirmed the subharmonic response was ubiquitous. They further separated the DTC phase from a thermal one by varying \(g\), and used a time-reversed echo protocol to distinguish intrinsic time-crystalline stability from mere decoherence.
This is among the cleanest demonstrations that a non-equilibrium phase of matter with no equilibrium analogue can be engineered and certified on a programmable quantum computer, with the spectrum-wide averaging that eigenstate order demands. Honest limits: it is a 1D chain of finite length observed over finite time, so “indefinite” stability is an extrapolation; the drive is idealized digital gates; and finite coherence still bounds the horizon. It stands as the benchmark for programmable Floquet-phase engineering and MBL certification on hardware. arXiv:2107.13571
Semeghini et al. 2021 — Rydberg topological spin liquid
A 219-atom Rydberg simulator used to prepare and probe a topological quantum spin liquid of the toric-code (\(\mathbb{Z}_2\)) type — a phase with no local order parameter, defined instead by long-range entanglement and topological string observables. Semeghini et al., “Probing topological spin liquids on a programmable quantum simulator,” Science 374, 1242 (2021) [arXiv:2104.04119 · Science 374, 1242]. It is the analog counterpart to the digital toric-code preparation elsewhere on this map.
A quantum spin liquid is a magnet that, even at zero temperature, never freezes into an ordered pattern: frustration keeps the spins fluctuating in a massively entangled superposition. The target here is the toric-code / \(\mathbb{Z}_2\) spin liquid — the same universality class as Kitaev’s exactly solvable model — whose order is topological, invisible to any local measurement and detectable only through non-local string operators and topological degeneracy. It hosts anyonic excitations and is a candidate substrate for fault-tolerant qubits, which is why realising it in a controlled system is a long-standing goal.
Atoms are placed not on lattice sites but on the links of a kagome lattice, and driven by the same Rydberg Hamiltonian \(H=\tfrac{\Omega}{2}\sum_i\sigma^x_i-\Delta\sum_i n_i+\sum_{i\lt j}V_{ij}n_i n_j\). Tuned into the blockade regime, an excited atom (a Rydberg “dimer”) forbids its neighbours from exciting, so the allowed configurations map onto the close-packed dimer coverings of the kagome links. The quantum superposition of these coverings, generated coherently by \(\Omega\), is precisely the \(\mathbb{Z}_2\) spin-liquid wavefunction — frustrated, locally disordered, topologically ordered.
An analog array of 219 \(^{87}\)Rb atoms in optical tweezers, arranged on the links of a kagome lattice and driven collectively into Rydberg states. The state is prepared by an adiabatic sweep of \(\Delta\) and \(\Omega\), then interrogated in place — no gates, the many-body ground state is engineered directly through the geometry and the drive.
They adiabatically prepared the candidate spin-liquid state and then measured topological string operators — closed and open Wilson-loop-like products of Rydberg occupations that diagnose \(\mathbb{Z}_2\) order. The closed-loop parity operators sat near their topologically ordered values while the corresponding open strings vanished, the signature that the correlations are non-local and topological rather than conventional. They further showed the observables were robust to deformations of the loop, consistent with genuine topological order rather than a symmetry-breaking artefact.
This is one of the first controlled experimental realisations of a toric-code-type topological spin liquid, prepared and diagnosed through the non-local observables that define it — a route toward encoding topologically protected quantum information in an atomic array. Honest limits: the state is prepared adiabatically over finite time with finite fidelity, so it is an approximate rather than exact spin liquid; the diagnostics are string expectation values rather than a full demonstration of anyon braiding or topological degeneracy; and it is a specific \(\mathbb{Z}_2\) instance. Still, it moved topological spin liquids from theory into programmable hardware. arXiv:2104.04119
Satzinger et al. 2021 — toric-code order on qubits
A superconducting processor used to prepare the toric-code ground state by an explicit quantum circuit and then verify its topological order two ways: measuring the topological entanglement entropy and performing anyon interferometry to read out braiding statistics. Satzinger et al. (Google Quantum AI), “Realizing topologically ordered states on a quantum processor,” Science 374, 1237 (2021) [arXiv:2104.01180 · Science 374, 1237]. The digital sibling of the Rydberg spin-liquid result — same target physics, opposite methodology.
Kitaev’s toric code is the paradigmatic exactly-solvable model of topological order,
with mutually commuting vertex (\(A_v\)) and plaquette (\(B_p\)) stabilizers. Its ground state has long-range entanglement and no local order; its excitations are anyons (\(e\) and \(m\) charges) whose mutual statistics are neither bosonic nor fermionic — braiding one around the other multiplies the wavefunction by \(-1\). These features are the foundation of topological quantum error correction.
This is a digital preparation. Rather than cooling into a ground state, the authors compile a deterministic quantum circuit that drives the qubits into the toric-code state — a shallow sequence of Hadamard and CNOT gates that sets every stabilizer \(A_v,B_p\) to \(+1\). Anyons are then created, moved, and braided by applying local Pauli operators along chosen paths, and the resulting phases are measured by interferometry. The full power of gate control is used: state preparation, excitation, transport, and readout are all explicit circuit operations.
The experiment used 31 superconducting qubits of Google’s Sycamore processor to host the toric-code lattice, plus ancillas for stabilizer readout and the interferometry protocol. Being fully gate-programmable is what lets them do things the analog route cannot: deterministically create individual anyons at chosen locations and drag them along arbitrary loops.
Three results. They prepared the toric-code ground state and measured a topological entanglement entropy near the ideal \(-\ln 2\), the entropic fingerprint of \(\mathbb{Z}_2\) topological order. They performed anyon interferometry, braiding \(e\) around \(m\) and detecting the expected mutual \(-1\) statistics. And they explored practical corollaries — logical-state injection into the code and the decay of a non-local order parameter — connecting the demonstration to quantum error correction.
By preparing topological order with a known circuit and then confirming it through both entanglement entropy and braiding statistics, this went a step beyond string-observable diagnostics: it manipulated the anyons directly. Together with the Rydberg spin-liquid paper it shows the same topological phase reached by two independent platforms. Honest limits: the state is circuit-prepared, so it is not a thermodynamic ground state and its stability is bounded by gate fidelity and circuit depth; 31 qubits is a small code; and the “braiding” is of Abelian anyons, not the non-Abelian kind needed for topological computation. It remains a landmark for topological order and error correction on hardware. arXiv:2104.01180
Ebadi et al. 2021 — quantum phases & Ising criticality
A 256-atom two-dimensional Rydberg simulator used to map out the quantum phase diagram of a programmable Ising magnet and to probe its quantum critical dynamics — including a direct test of Kibble–Zurek scaling across a (2+1)-dimensional Ising transition. Ebadi et al., “Quantum phases of matter on a 256-atom programmable quantum simulator,” Nature 595, 227 (2021) [arXiv:2012.12281 · Nature 595, 227]. The scale-up of Bernien’s 51-atom chain into two dimensions.
The target is quantum magnetism in 2D and the quantum phase transition separating a disordered paramagnet from antiferromagnetically ordered phases at zero temperature. Near the critical point the physics is governed by the (2+1)-dimensional Ising universality class: correlation length and relaxation time diverge, and static and dynamic critical exponents take universal values independent of microscopic detail. Driving the system through such a transition at finite speed is the setting for Kibble–Zurek physics, where the density of frozen-in defects scales as a power of the sweep rate.
The same Rydberg Hamiltonian \(H=\tfrac{\Omega}{2}\sum_i\sigma^x_i-\Delta\sum_i n_i+\sum_{i\lt j}V_{ij}n_i n_j\) is now realised on a two-dimensional square array. By setting the ratio of blockade radius \(R_b\) to lattice spacing, the blocked region spans different numbers of neighbours, programming distinct broken-symmetry ground states — checkerboard, striated, and star antiferromagnetic orders. Sweeping \(\Delta\) at controlled speed carries the system across the paramagnet-to-antiferromagnet transition, so the geometry selects the phase and the sweep probes the criticality.
An analog simulator of deterministically loaded \(^{87}\)Rb atoms in a 2D optical-tweezer array, scanned across system sizes from 64 up to 256 atoms. Larger sizes are what make sharp criticality and Kibble–Zurek scaling accessible: universal exponents only emerge as the system approaches the thermodynamic limit, so the atom count is itself a scientific resource here, not just a headline number.
They benchmarked the platform by preparing the checkerboard and other antiferromagnetic phases with high fidelity and mapping the 2D phase diagram. Then, sweeping the detuning at a range of speeds through the transition, they measured quantum critical dynamics consistent with the (2+1)D Ising universality class and extracted the Kibble–Zurek scaling of the correlation length against sweep rate. Using system sizes up to 256 atoms let them see the finite-size approach to the universal critical behaviour.
This demonstrated quantitative many-body physics — universal critical exponents, defect scaling — from a programmable analog quantum simulator at a scale (256 atoms) far beyond exact classical simulation, establishing Rydberg arrays as a tool for studying quantum criticality, not just preparing states. Honest limits: it is a specific van-der-Waals Ising model in 2D; adiabaticity and finite coherence bound how close to the true critical point one can sweep; and universal-exponent extraction from finite sizes carries fitting uncertainty. It is the reference point for programmable quantum simulation of phases and criticality at scale. arXiv:2012.12281
King et al. 2021 — annealer scaling advantage, frustrated magnets
A quantum annealer used to simulate the dynamics of geometrically frustrated quantum magnets, and — most notably — to demonstrate a scaling advantage: the annealer reached equilibrium faster than a leading classical method (path-integral Monte Carlo), with the margin growing with system size and inverse temperature. King et al. (D-Wave), “Scaling advantage over path-integral Monte Carlo in quantum simulation of geometrically frustrated magnets,” Nat. Commun. 12, 1113 (2021) [arXiv:1911.03446 · Nat. Commun. 12, 1113].
Geometric frustration occurs when a lattice’s bonds cannot all be satisfied at once — the classic case being antiferromagnetic couplings on triangular motifs, where no spin arrangement minimizes every bond. Adding quantum fluctuations via a transverse field gives a transverse-field Ising model on a frustrated lattice,
whose competition between frustration and quantum tunnelling produces rich ordering, order-by-disorder, and slow equilibration — a hard problem for classical samplers, and the setting where path-integral Monte Carlo (PIMC) is the standard classical workhorse.
A quantum annealer natively implements exactly this transverse-field Ising Hamiltonian: the couplings \(J_{ij}\) are programmed into the qubit interconnects and the transverse field \(\Gamma\) is the annealing control. The frustrated lattice (a distorted triangular / “Villain”-type geometry) is embedded into the device’s fixed qubit-connectivity graph, and the machine is used not to find a ground state but as an analog quantum simulator of dynamics — quenching \(\Gamma\) and watching how the frustrated magnet relaxes toward equilibrium in real time.
The device is the D-Wave 2000Q annealer built from roughly 2000 (~2048) superconducting flux qubits; the frustrated-magnet simulations used up to about 1440 qubits, with dynamics resolved at microsecond timescales. This is an analog platform with a fixed sparse connectivity graph — not a universal gate machine — so the class of Hamiltonians it can express is limited to what embeds into that graph.
They quenched the transverse field and measured how the frustrated magnet equilibrated as a function of system size and temperature, comparing the annealer’s dynamics against path-integral Monte Carlo run on classical CPUs. The annealer showed a dynamical scaling advantage: it approached equilibrium faster, and the advantage grew with both system size and inverse temperature, amounting to a reported speedup of order a million-fold over an efficient classical implementation for the largest, coldest instances.
This is one of the more credible claims of quantum-simulation advantage from an annealer, and it reframed the device as a robust statistical-physics simulator for frustrated quantum matter rather than merely an optimizer. The advantage is carefully scoped: it is a dynamical/scaling advantage against a specific classical algorithm (local-update PIMC), not a universal computational-complexity separation, and smarter classical update schemes (e.g. cluster methods) could narrow the gap. Honest limits: fixed sparse connectivity, embedding overhead, and a restricted Hamiltonian class. Within those bounds, it stands as a landmark analog simulation of frustrated magnetism at scale. arXiv:1911.03446
O'Malley et al. 2016 — molecular hydrogen, the first scalable VQE
This is the paper that put the variational quantum eigensolver (VQE) on real superconducting hardware: the full ground-state energy curve of the H\(_2\) molecule, computed on two qubits without the exponentially expensive pre-compilation that had crippled earlier demonstrations. O'Malley, Babbush, Kivlichan et al., "Scalable Quantum Simulation of Molecular Energies," Phys. Rev. X 6, 031007 (2016) [arXiv:1512.06860 · Phys. Rev. X 6, 031007]. It is the template every later chemistry-VQE slide in this family specialises.
The target is the Born–Oppenheimer electronic ground-state energy of H\(_2\) as a function of the inter-nuclear separation \(R\) — the potential-energy curve whose minimum is the bond length and whose asymptote is the dissociation limit. In a minimal (STO-3G) basis, H\(_2\) has two electrons in two spatial orbitals; the electronic Hamiltonian is the standard second-quantized form \(H=\sum_{pq}h_{pq}\,a_p^\dagger a_q+\tfrac12\sum_{pqrs} h_{pqrs}\,a_p^\dagger a_q^\dagger a_r a_s\). Getting the curve right to chemical accuracy (\(\approx1\) kcal/mol \(\approx1.6\) mHa) is the benchmark that separates a toy from a chemically meaningful calculation.
The fermionic Hamiltonian is reduced to two spin-orbitals and mapped to qubit operators via a Bravyi–Kitaev / parity-style encoding, collapsing to a compact two-qubit Hamiltonian that is a sum of tensored Pauli terms. VQE then prepares a parametrised trial state \(|\psi(\theta)\rangle=U(\theta)|\psi_0\rangle\) using a unitary coupled-cluster (UCC) single-excitation ansatz with a single variational angle \(\theta\), and minimises the energy expectation
with the quantum device measuring each Pauli term and a classical optimiser adjusting \(\theta\). The variational principle guarantees \(E(\theta)\ge E_0\), so the answer is a rigorous upper bound.
Two Xmon transmon superconducting qubits (the Martinis group, then at UC Santa Barbara / Google), driven by microwave single-qubit rotations and a dc-flux two-qubit entangling gate. For contrast the authors also ran the same problem with quantum phase estimation (PEA) on the same device, making this a head-to-head of the two paradigms on identical hardware.
They mapped the entire H\(_2\) dissociation curve point-by-point in \(R\). The key finding is a methodological one: VQE reached the correct dissociation energy to within chemical accuracy, whereas the phase-estimation run on the same noisy hardware did not — because VQE's variational feedback loop self-corrects against systematic control errors that PEA has no way to absorb. The circuits were short enough to run without the classical pre-diagonalisation that had made prior "quantum chemistry" demonstrations circular.
This established VQE as the practical route to chemistry on noisy, pre-fault-tolerant hardware, and its error-resilience argument became the field's founding rationale. Honest limits: a two-orbital, two-qubit problem is classically trivial — the achievement is the protocol, not the answer; the single-parameter ansatz does not generalise cheaply; and readout requires many repeated measurements per Pauli term. It is a proof of principle, correctly labelled as such, and the direct ancestor of Kandala 2017 and Arute 2020. arXiv:1512.06860
Kandala et al. 2017 — molecular ground-state energies up to BeH₂
Kandala, Mezzacapo, Temme et al., "Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets," Nature 549, 242–246 (2017) [arXiv:1704.05018 · Nature 549, 242]. This paper pushed VQE from a two-orbital toy up to BeH\(_2\) — the largest molecule simulated on a quantum computer at the time — by abandoning chemically-motivated circuits in favour of an ansatz built from the processor's native gates.
The targets are the electronic ground-state energies of H\(_2\), LiH and BeH\(_2\) along their bond-stretch coordinates, plus the ground state of a four-spin antiferromagnetic Heisenberg model. Each molecule is a larger electronic-structure problem than the last: more orbitals, more electrons, and a qubit Hamiltonian with over a hundred Pauli terms for BeH\(_2\). Recovering the correct binding curves is the physics goal.
The molecular Hamiltonians are reduced (frozen core, parity mapping with qubit tapering) to few-qubit qubit Hamiltonians. The decisive move is the hardware-efficient ansatz: rather than a UCC circuit whose depth explodes under the device's limited connectivity, the trial state is
alternating layers of single-qubit Euler rotations \(U_{\rm rot}\) with the processor's native entangling interaction \(U_{\rm ent}\). This keeps circuits shallow enough to survive decoherence, at the cost of no longer respecting particle-number or spin symmetries — a trade the paper makes explicitly.
Up to six qubits of a purpose-built seven-qubit IBM superconducting processor. H\(_2\) used 2 qubits, LiH and BeH\(_2\) up to 6; the Heisenberg-magnet demonstration used four. Verified: six qubits is the maximum in this work.
They variationally optimised each molecule's energy on hardware, tracing the H\(_2\), LiH and BeH\(_2\) dissociation curves. The energies followed the correct shape but sat above chemical accuracy: hardware noise, and the symmetry-breaking freedom of the hardware-efficient ansatz, left a residual error the raw runs could not close. The paper pairs the demonstration with an explicit account of these error sources rather than claiming chemical precision.
It set the device-native baseline for NISQ chemistry and became one of the most-cited VQE experiments. Honest limits, stated in the paper: the six-qubit results are not yet chemically accurate; the expressive hardware-efficient ansatz suffers barren-plateau trainability problems as it scales; and symmetry-breaking can drive the optimiser to unphysical states. It marks the moment VQE became a genuine hardware programme — and exposed the noise wall that Arute 2020's error mitigation would target next. arXiv:1704.05018
Arute et al. 2020 — Hartree–Fock of hydrogen chains on Sycamore
Google AI Quantum & collaborators (Arute, Arya, Babbush et al.), "Hartree–Fock on a superconducting qubit quantum computer," Science 369, 1084–1089 (2020) [arXiv:2004.04174 · Science 369, 1084]. This was the largest quantum-chemistry calculation yet run on a quantum processor — hydrogen chains up to H\(_{12}\) and a real reaction mechanism — using roughly twice the qubits and ten times the gates of any prior chemistry experiment.
The targets are the Hartree–Fock (mean-field) binding-energy curves of the linear hydrogen chains H\(_6\), H\(_8\), H\(_{10}\) and H\(_{12}\), and the isomerization pathway of diazene (N\(_2\)H\(_2\), the cis–trans rearrangement) — an actual chemical-reaction mechanism, not just a static energy. Reproducing the correct energy along each geometry is the physics validation.
Hartree–Fock is a non-interacting-fermion problem, so its exact ground state is a single Slater determinant reachable by a network of Givens rotations — two-qubit number-conserving gates that implement an arbitrary single-particle basis change:
Variationally optimising the \(\theta_{pq}\) rotates the occupied orbitals into the true Hartree–Fock solution. Because the state is classically tractable, the experiment is a rigorous, checkable benchmark of the hardware at unprecedented scale rather than a claim of quantum advantage.
The Google Sycamore superconducting processor, using up to 12 qubits (one per spin-orbital for H\(_{12}\)). Verified: 12 qubits is the maximum, roughly double the largest previous chemistry demonstration.
They optimised the Givens-rotation ansatz on hardware to recover the H\(_6\)–H\(_{12}\) binding curves and the diazene isomerization energetics. The headline is the error-mitigation result: post-processing that enforces fermionic \(N\)-representability — physical constraints on the measured one-particle density matrix (correct particle number, positivity, idempotency) — dramatically raised the effective fidelity, turning noisy raw data into chemically sensible curves. This validity-projection strategy is the paper's transferable contribution.
It demonstrated that careful error mitigation lets NISQ hardware scale chemistry to a dozen qubits and tackle a genuine reaction. The candid limit: Hartree–Fock is classically easy — the experiment deliberately simulates a problem it can check, so it proves hardware capability and mitigation, not a quantum-chemical speedup; capturing the electron correlation that makes chemistry hard remained out of reach. It set the NISQ-chemistry scale record and the error-mitigation playbook for what came after. arXiv:2004.04174
Dumitrescu et al. 2018 — deuteron binding energy, first nuclear VQE
Dumitrescu, McCaskey, Hagen et al., "Cloud Quantum Computing of an Atomic Nucleus," Phys. Rev. Lett. 120, 210501 (2018) [arXiv:1801.03897 · Phys. Rev. Lett. 120, 210501]. This is the first nuclear-physics calculation on a quantum computer: the binding energy of the deuteron — the proton–neutron bound state — obtained via VQE on cloud-accessed superconducting qubits.
The deuteron is the simplest nucleus, a single proton–neutron pair bound by the strong force, with an accepted binding energy of about \(-2.22\) MeV. Reproducing this number from a first-principles nuclear Hamiltonian is the physics target — a genuine strong-interaction bound state rather than an electronic-structure problem.
They build the nuclear Hamiltonian from leading-order pionless effective field theory in a discrete harmonic-oscillator basis of \(N\) states, then map it to qubits by Jordan–Wigner, giving an \(N\)-qubit Hamiltonian \(H_N\). The trial state is a low-depth unitary coupled-cluster ansatz with only one or two variational angles, and VQE minimises
with a classical extrapolation in the basis size \(N\to\infty\) to recover the continuum binding energy.
Cloud-accessed superconducting processors — IBM QX5 and Rigetti 19Q — using just 2 qubits for the \(N=2\) truncation and 3 qubits for \(N=3\). Verified: two-to-three qubits, run entirely over the cloud, no on-site hardware.
They computed \(E_N\) for the \(N=2\) and \(N=3\) truncations on hardware and extrapolated. With error mitigation the extrapolated binding energy landed within a few percent of the physical \(-2.22\) MeV. The demonstration's point is the full pipeline — nuclear EFT \(\rightarrow\) qubit Hamiltonian \(\rightarrow\) cloud VQE \(\rightarrow\) extrapolated observable — executed end-to-end on commercial quantum hardware.
It opened nuclear structure as a quantum-computing application and showed a complete VQE workflow could run on minimal, publicly available qubits. Honest limits: two/three basis states is a drastic truncation of a classically trivial system; the result leans on classical basis extrapolation and error mitigation; and the two-qubit deuteron is many orders of magnitude from a nucleus where quantum computing could beat classical ab initio methods. A milestone of scope — first nucleus — not of computational advantage. arXiv:1801.03897
Steinhauer 2016 — analog Hawking radiation from a sonic horizon
Jeff Steinhauer, "Observation of quantum Hawking radiation and its entanglement in an analogue black hole," Nature Physics 12, 959–965 (2016) [arXiv:1510.00621 · Nat. Phys. 12, 959]. Hawking's 1974 prediction that black holes radiate is essentially untestable at astrophysical scales; this experiment realises its acoustic analogue in a Bose–Einstein condensate and reports entangled, thermal Hawking pairs across a sonic horizon.
Hawking radiation arises when quantum vacuum fluctuations straddle an event horizon: one partner escapes as thermal radiation, the other is trapped, and the pair is entangled. Unruh's insight (1981) is that a fluid flowing faster than its own sound speed forms an acoustic ("sonic") horizon for phonons that is mathematically the same curved-spacetime problem — sound cannot escape upstream from the supersonic region, exactly as light cannot escape a black hole.
In a flowing condensate the phonon field obeys a massless-scalar equation in an effective acoustic metric set by the local flow velocity \(v\) and sound speed \(c\); the horizon sits where \(|v|=c\). The analog Hawking temperature is fixed by the analog surface gravity — the velocity gradient at the horizon,
so a steep, controlled flow step produces a measurable thermal phonon flux.
An elongated (quasi-1D) \(^{87}\)Rb Bose–Einstein condensate of a few thousand atoms, with a sweeping step-like potential that drives part of the cloud supersonic, creating the analog horizon. Verified species: rubidium-87.
From density–density correlations in absorption images, Steinhauer read out correlated phonon pairs emitted on opposite sides of the horizon — the Hawking partner and its infalling twin. Two quantitative results: the spectrum is approximately thermal (a Hawking temperature can be assigned), and for the high-energy pairs the correlations violate a Cauchy–Schwarz-type classical bound, certifying genuine entanglement; the low-energy pairs were not found to be entangled.
This is the first laboratory evidence that Hawking-type emission is real and quantum (entangled), not merely a classical scattering artefact — a striking cross-domain confirmation of Hawking's mechanism. Honest limits: it is an analogue, an acoustic metric in a cold gas, so it tests the kinematics of horizon physics, not real gravity; the signal is delicate and the thermality approximate; and the entanglement claim, while landmark, drew active scrutiny over statistics and interpretation. A defining result of the analog-gravity programme. arXiv:1510.00621
Zenesini et al. 2024 — false-vacuum decay by bubble nucleation
Zenesini, Berti, Cominotti et al., "False vacuum decay via bubble formation in ferromagnetic superfluids," Nature Physics 20, 558–563 (2024) [arXiv:2305.05225 · Nat. Phys. 20, 558]. This is the first laboratory observation of false-vacuum decay — the bubble-nucleation process behind cosmological phase transitions and the long-term fate of our own universe's vacuum — realised in a ferromagnetic atomic superfluid. It is the real experiment that our VACUA work reproduces numerically.
A quantum field trapped in a metastable ("false") minimum of its potential, separated by a barrier from a lower true minimum, decays by quantum/thermal nucleation of bubbles of true vacuum. Bubbles above a critical size lower the total energy by growing and coalescing, converting false vacuum to true — the Coleman–De Luccia mechanism underlying early-universe phase transitions. It had never been directly observed.
Following the Fialko/Billam cold-atom proposals, two coupled superfluid components form a ferromagnetic spin mixture; the relative phase \(\theta\) between them behaves like a relativistic scalar field whose effective potential — shaped by the coherent coupling and a detuning that tilts the wells — has a false minimum and a true minimum:
a tilted washboard. Preparing the mixture in the metastable well seeds a genuine false vacuum; a spin bubble is a localised patch that has rolled to the true minimum.
A \(^{23}\)Na (sodium) Bose–Einstein condensate — a two-state superfluid spin mixture — at the BEC group in Trento. Verified species: sodium-23.
They prepared the false vacuum and directly imaged spin bubbles nucleating and growing as the system decayed toward the true vacuum. By repeating the experiment they extracted the nucleation rate and its dependence on the coupling strength and on temperature, comparing against theoretical (instanton / thermal-activation) predictions and supporting numerical simulations. This turns an abstract field-theory process into measured, statistically characterised bubble events.
It is the first experimental realisation of false-vacuum decay, making a keystone of cosmology and quantum field theory accessible in a tabletop quantum gas. Honest limits: the mapping to a relativistic scalar field is an effective, time-averaged description of a non-relativistic superfluid; boundary effects and finite temperature complicate the clean instanton picture; and it is an analogue, so it probes the mechanism, not literal cosmological vacuum decay. Directly cross-referenced by our VACUA reproduction (slide 9). arXiv:2305.05225
Landsman et al. 2019 — verified quantum scrambling via teleportation
Landsman, Figgatt, Schuster et al., "Verified quantum information scrambling," Nature 567, 61–65 (2019) [arXiv:1806.02807 · Nature 567, 61]. Scrambling — the spreading of local information across a many-body system's entanglement — is the process behind black-hole thermalisation and the traversable-wormhole teleportation picture. This experiment gives the first unambiguous detection of scrambling, using teleportation as the probe.
Scrambling disperses information encoded in a few local degrees of freedom into non-local many-body correlations, so it can no longer be read out locally — the microscopic mechanism of thermalisation, and, in holographic duality, of how black holes process infalling information. Its standard diagnostic is the decay of the out-of-time-ordered correlator (OTOC), but OTOC decay alone is ambiguous: ordinary decoherence mimics it. Distinguishing true scrambling from noise is the challenge.
They implement the Yoshida–Kitaev protocol, in which successful teleportation of a qubit state through the circuit is a positive signature of scrambling that decoherence cannot fake — decoherence destroys the teleportation, scrambling enables it. A tunable unitary interpolates between scrambling and non-scrambling dynamics, and the teleportation fidelity reports which regime the circuit is in, cleanly separating scrambling from mere information loss.
A 7-qubit trapped-ion processor of \(^{171}\)Yb\(^+\) ions (Monroe group, University of Maryland), running a tunable 3-qubit scrambling unitary embedded in the teleportation circuit. Verified: seven qubits, ytterbium-171.
They ran the teleportation circuit for scrambling and non-scrambling settings and measured teleportation fidelities of about 80%, well above the classical threshold — a positive, decoherence-immune confirmation that the scrambling unitary genuinely disperses information. From the same data they bounded the OTOC, tying the teleportation signal quantitatively to the standard scrambling diagnostic.
This is the first experiment to detect scrambling unambiguously, resolving the OTOC-vs-decoherence confound, and it connects tabletop quantum computing to black-hole information and traversable-wormhole ideas. Honest limits: seven qubits with a three-qubit scrambler is a small system — it validates the protocol, not black-hole-scale dynamics; the "wormhole" language is an analogy via holographic duality, not gravity; and fidelities near 80% still carry appreciable hardware error. A conceptual landmark linking quantum information and gravity. arXiv:1806.02807
Viermann et al. 2022 — cosmological particle production in an expanding-universe analog
Viermann, Sparn, Liebster et al. (Oberthaler group), "Quantum field simulator for dynamics in curved spacetime," Nature 611, 260–264 (2022) [arXiv:2202.10399 · Nature 611, 260]. This experiment builds a controllable expanding-universe analogue in a cold gas and observes cosmological particle production — the pair creation of quanta driven by a time-dependent spacetime, the same mechanism that seeds primordial structure.
A quantum field in an expanding universe cannot keep its vacuum: the stretching metric mixes positive- and negative-frequency modes, so an initial vacuum evolves into a state populated with pairs of particles (Sakharov / Gibbons–Hawking particle production). In cosmology this converts vacuum fluctuations into the density perturbations that grow into galaxies. Reproducing this — including its characteristic Sakharov oscillations in the fluctuation spectrum — is the physics target.
In a 2D condensate the phonon field propagates in an effective Friedmann–Robertson–Walker-like metric whose scale factor is set by the sound speed \(c_s\propto\sqrt{g\,n}\). Because the interaction \(g\) is tunable by a Feshbach resonance, sweeping it in time makes \(c_s(t)\) — and hence the analog scale factor \(a(t)\) — follow a programmed expansion law. The condensate shape can also encode spatially curved (positive/negative) metrics, so both time-dependence and curvature are engineered.
A 2D "pancake" Bose–Einstein condensate of about 23,000 \(^{39}\)K (potassium) atoms, in Heidelberg, with the interaction magnetically tuned to drive the analog expansion. Verified: potassium, roughly 23,000 atoms in a quasi-2D geometry.
They imprinted phonon wavepackets and swept the sound speed to enact expanding (and curved) spacetimes, then measured the resulting phonon power spectrum. They observed cosmological particle production — the amplification of fluctuations expected from the time-dependent metric — including the predicted Sakharov oscillations, quantitatively matching quantum-field-theory-in-curved-spacetime calculations.
It is a programmable quantum-field simulator for curved spacetime, demonstrating cosmological particle production and Sakharov oscillations in the lab and adding a tunable expansion/curvature knob absent from earlier analog experiments. Honest limits: it is an analogue — a non-relativistic superfluid engineered to mimic an FRW metric, testing the field-theory kinematics rather than real cosmology; it is 2D; and the accessible expansion histories and mode range are bounded by atom number and interaction control. A leading result in analog cosmology. arXiv:2202.10399
VACUA (ours) — false-vacuum decay by a classical differentiable spectral solver
VACUA is our own numerical reproduction of cold-atom false-vacuum decay: a classical, differentiable spectral solver (implemented in JAX) that integrates the coupled-superfluid field equations and reproduces the bubble-nucleation phenomenology of the Fialko/Billam proposal and the Zenesini et al. 2024 experiment (slide 6). It is explicitly a classical simulation — status numerical, not quantum hardware — offered as a validation and design tool, not a quantum-advantage claim.
The target is the same as Zenesini's: a field held in a metastable false minimum that decays to the true vacuum through the nucleation and growth of critical bubbles. We aim to reproduce, from first-principles classical field dynamics, the qualitative bubble life-cycle and the scaling of the nucleation rate with the control parameters (coupling strength, detuning/tilt, effective temperature) that the experiment measures.
We evolve the two coupled superfluid components under the mean-field Gross–Pitaevskii equations with a coherent coupling and detuning; after time-averaging, the relative phase \(\theta\) obeys an effective relativistic-scalar (sine-Gordon-type) equation in the tilted double-well potential \(V(\theta)\sim-\lambda\cos\theta-\epsilon\theta\), exactly the metastable landscape of the experiment. Fluctuations that seed nucleation are injected semiclassically — a truncated-Wigner-style stochastic sampling of the initial field — standing in for the quantum/thermal vacuum noise.
No quantum device: a classical split-step Fourier (pseudo-spectral) integrator written in JAX, run on GPU. Being written in JAX, the whole solver is differentiable end-to-end, so gradients of observables with respect to the physical control parameters are available by automatic differentiation.
We integrated the coupled fields from a prepared false-vacuum state and observed spin/phase bubbles nucleating and growing, then accumulated statistics over many stochastic realisations to extract the nucleation rate and its dependence on coupling and effective temperature — reproducing the trends reported by Zenesini et al. The differentiability lets us run gradient-based parameter inference, fitting solver controls to a target nucleation behaviour rather than scanning by brute force.
Our contribution is a fast, differentiable, reproducible classical benchmark for the false-vacuum-decay platform — useful for pre-designing experiments, fitting parameters, and sanity-checking the analog mapping against the real \(^{23}\)Na experiment. The ceiling is explicit and honest: it is a classical mean-field plus semiclassical-noise model, so it does not access genuine quantum tunnelling or the true entanglement structure of vacuum decay — nucleation is seeded stochastically, not by a real quantum instanton — and it therefore offers no quantum advantage and does not supersede the Zenesini measurement. It is a numerical companion to that experiment, cross-referenced to it, and clearly labelled as our classical reproduction. [ours · classical JAX spectral solver · cf. arXiv:2305.05225]
HHL 2009 — quantum linear systems, the root of the whole genre
Every quantum PDE/ODE algorithm — the Vlasov seed included — is a descendant of one result: Harrow, Hassidim & Lloyd's algorithm for solving \(A\mathbf x=\mathbf b\) in time polylogarithmic in the dimension. It is also the source of every caveat the seed inherits. [arXiv:0811.3171, Phys. Rev. Lett. 103, 150502 (2009)]
Given an \(N\times N\) sparse matrix \(A\) and a vector \(\mathbf b\), produce the solution of \(A\mathbf x=\mathbf b\) — but as a quantum state \(|x\rangle\propto A^{-1}|b\rangle\), from which one reads an expectation \(\langle x|M|x\rangle\). Classically, conjugate gradient costs \(\sim O(N s\sqrt{\kappa})\); HHL targets \(\operatorname{polylog}(N)\).
Expand \(|b\rangle=\sum_j\beta_j|u_j\rangle\) in the eigenbasis of \(A\) (\(A|u_j\rangle=\lambda_j|u_j\rangle\)). Run phase estimation on \(e^{iA t}\) to write each \(\lambda_j\) into an ancilla register; apply a controlled rotation that multiplies each branch by \(1/\lambda_j\); uncompute the eigenvalue register and post-select the ancilla. What survives is
Phase estimation needs Hamiltonian simulation of \(A\) — which is why sparsity and \(e^{-iAt}\) sit at the core.
The runtime is
with \(s\) the sparsity, \(\kappa=\lambda_{\max}/\lambda_{\min}\) the condition number, \(\epsilon\) the accuracy — exponential in \(N\) versus the best classical solvers.
- State preparation. \(|b\rangle\) must be loadable efficiently — typically requires QRAM (slide 8).
- Sparsity / simulability. \(A\) must be sparse (or block-encodable) so \(e^{-iAt}\) is cheap.
- Condition number. Cost \(\propto\kappa\) (originally \(\kappa^2\)); ill-conditioned \(A\) destroys the speedup.
- Read-out. You get \(|x\rangle\) in amplitudes, not \(\mathbf x\). Extracting a component needs tomography \(\sim O(N)\); only a global \(\langle x|M|x\rangle\) is cheap. (And when \(A\) is low-rank with sampling access, Tang's dequantization matches HHL classically — so the exponential speedup is conditional.)
Time-evolving a discretised PDE is, at bottom, repeatedly solving linear systems — the linear-ODE quantum solvers the seed relies on are built directly on HHL and its successors. More importantly, the shape of the Vlasov caveats is the shape of HHL's fine print: state prep = the QRAM assumption, sparsity = the central-difference stencil, condition number = the stability of \(A(t)\), and "output only \(P_\nu(k)\), never the full \(f\)" = HHL's read-out limit. Learn HHL and you have learned, in miniature, exactly where every quantum-PDE speedup can leak.
Hamiltonian simulation, Lloyd → Trotter → qubitization — the engine of \(e^{-iHt}\)
Once a PDE is Schrödinger-shaped, the entire task is to build the circuit for \(e^{-iHt}\). Three decades of progress took this from crude product formulas to a provably optimal algorithm — and the seed's cost bound rides directly on the optimal result.
Given a description of a Hamiltonian \(H\) (a sum of terms, or a sparse-matrix oracle) and a time \(t\), output a circuit that applies \(e^{-iHt}\) to accuracy \(\epsilon\). This is the primitive that evolves the neutrino state in the Vlasov algorithm and the field in the wave/Maxwell algorithms.
Lloyd 1996 (Trotter). For \(H=\sum_j H_j\) with each \(H_j\) simulable, split the exponential:
Simple, but the \(1/\epsilon\) scaling is poor. Berry–Ahokas–Cleve–Sanders 2007 extended this to general sparse \(H\) via graph-coloring decompositions and high-order Suzuki formulas; LCU / truncated-Taylor methods (Berry et al. 2015) then achieved exponentially better accuracy, \(O(\log(1/\epsilon))\). Low & Chuang 2017 (QSP / qubitization). The optimal construction builds a block-encoding — a unitary \(U\) that hides \(H/\alpha\) in a corner,
whose associated walk operator has eigenphases \(\arccos(E_k/\alpha)\). Quantum signal processing then applies a polynomial of degree \(d=O(\alpha t+\log(1/\epsilon))\) to those eigenphases to synthesise \(e^{-iHt}\) exactly to order \(\epsilon\).
Qubitization achieves query cost
additive in \(t\) and \(\log(1/\epsilon)\), matching the no-fast-forwarding lower bound — you cannot do better. Here \(\alpha\) is the block-encoding normalization (\(\alpha\ge\|H\|\)).
- Trotter error accumulates and needs \(H\) as a sum of few simulable terms; the older methods scale badly in \(\epsilon\).
- The cost is \(\alpha t\), not \(\|H\|t\). A loose block-encoding (large \(\alpha\)) inflates everything — building a tight encoding is the real work.
- Ancillas + fault tolerance. QSP/qubitization need extra qubits and deep, error-corrected circuits.
This is Part 6 of the seed. Each factor \(e^{-i\Delta_t H_{i_t}}\) in the Vlasov time-evolution is implemented by block-encoding the sparse antisymmetric \(A(t)\) (its central-difference structure gives sparsity; a QRAM supplies the CDM force entries) and then applying qubitization. The optimal \(O(\alpha t+\log(1/\epsilon))\) cost is precisely what lets Miyamoto claim \(\widetilde{O}(n_{\rm gr}+n_t)\) queries rather than something polynomial in the grid. Every algorithm on this deck that says "then Hamiltonian-simulate" is cashing out this exact lineage.
Amplitude estimation — the quadratic read-out speedup, \(O(1/\epsilon)\) not \(O(1/\epsilon^2)\)
A quantum algorithm can prepare a state encoding the answer — but a naïve measurement wastes it. Amplitude estimation is the primitive that extracts a number from an amplitude-encoded state quadratically faster than sampling, and it is exactly how the seed reads out \(P_\nu(k)\). [Brassard, Høyer, Mosca & Tapp 2002; Montanaro, Proc. R. Soc. A 471, 20150301 (2015)]
Given a circuit \(\mathcal A\) that prepares \(|\psi\rangle=\mathcal A|0\rangle\) with a "good" component of amplitude \(\sqrt{a}\) (so a measurement yields the marked outcome with probability \(a\)), estimate \(a\) — or more generally a mean \(\mathbb E[g(X)]\) — to additive error \(\epsilon\). Classical Monte Carlo needs \(O(1/\epsilon^2)\) samples (Chernoff); can we do better?
Write \(a=\sin^2\theta\). The Grover/amplitude-amplification operator \(Q=-\mathcal A\,S_0\,\mathcal A^{\dagger}\,S_\chi\) rotates in the 2-D "good/bad" plane with eigenphases \(\pm2\theta\). Running phase estimation on \(Q\) extracts \(\theta\) — hence \(a\) — using \(O(1/\epsilon)\) applications of \(\mathcal A\):
Montanaro 2015 generalises this to a full quantum speedup of Monte Carlo: the mean of a bounded-variance random variable is estimated to error \(\epsilon\) in \(\widetilde O(1/\epsilon)\) queries — a near-quadratic improvement over the classical \(O(\sigma^2/\epsilon^2)\).
A quadratic speedup in the precision \(\epsilon\): \(O(1/\epsilon)\) coherent queries versus \(O(1/\epsilon^2)\) independent samples. This is polynomial (not exponential), but decisive whenever a target accuracy \(\epsilon\) is small — halving the error costs \(2\times\), not \(4\times\), the queries.
- Coherent access required. You need \(\mathcal A\) and \(\mathcal A^\dagger\) as circuits and must run them inside phase estimation — deep, fault-tolerant, not NISQ. (Iterative/maximum-likelihood AE variants cut ancillas but keep the \(1/\epsilon\) scaling.)
- Only quadratic. It never turns a hard problem easy — it sharpens read-out, nothing more.
- One number at a time. Estimating many quantities (a whole spectrum) multiplies the cost.
This is Part 7 of the seed made general. The neutrino power spectrum \(P_\nu(k)\) is a set of expectation values / overlaps of the evolved state \(|f(T)\rangle\); amplitude estimation reads each to accuracy \(\epsilon\) in \(O(1/\epsilon)\) queries — the \(O(1/\epsilon)\)-not-\(O(1/\epsilon^2)\) factor is what keeps the seed's total query count at \(\widetilde{O}(n_{\rm gr}+n_t)\). The very same primitive powers Gaitan's ODE integrator (slide 2). It is the universal answer to "how do you get a number out of a state whose \(2^n\) amplitudes you can never look at?"
Carleman linearization + QRAM honesty — the field's two soft undersides
Every quantum-PDE pipeline has two ends where the promised speedup is most fragile: getting the nonlinearity in (Carleman) and getting the data in (QRAM). This slide names both honestly — the second is where most "exponential speedups" quietly leak. [Carleman: Liu et al. 2020; QRAM: Jaques & Rattew, arXiv:2305.10310]
Because quantum evolution is linear, any nonlinear term must first be linearized. Carleman does this by lifting the state to its tensor-power tower \(\mathbf u,\mathbf u^{\otimes2},\dots\), giving an infinite block-tridiagonal linear ODE that is truncated at level \(N\) and solved with a quantum linear-ODE solver. Convergence hinges on the dissipation ratio
For \(R\ge\sqrt2\) it is provably intractable — so turbulent / strongly nonlinear regimes are out of reach (the seed avoids this entirely by physically linearizing to an external CDM force).
Quantum RAM answers a superposition of addresses: \(\sum_i\alpha_i|i\rangle|0\rangle\mapsto \sum_i\alpha_i|i\rangle|d_i\rangle\). Loading \(N\) classical numbers into amplitudes is the entry point for HHL, for state preparation, and for the seed's stored CDM force. The bucket-brigade architecture advertises query depth \(O(\log N)\) — the headline that makes "polylog" algorithms look free.
The \(O(\log N)\) query time hides \(O(N)\) query cost:
- Hardware. A QRAM over \(N\) cells needs \(\Theta(N)\) physical components, all coherent during the query.
- Error correction. An active (circuit) QRAM must error-correct all \(O(N)\) components every query — so the total work is \(O(N)\), and "most asymptotic quantum advantage disappears."
- Opportunity cost. The same \(O(N)\) controllers could just run a massively parallel classical algorithm and finish just as fast.
- Passive proposals (no per-query energy) rely on physically dubious assumptions; the survey finds "all non-circuit QRAM proposals fall short."
QRAM is not a solved primitive but an assumption. Any speedup that depends on loading \(O(N)\) arbitrary classical numbers into a state should be read as "fast given a device that must itself do \(O(N)\) coherent, error-corrected work" — which may erase the advantage it was invoked to provide. This is the field's soft underbelly.
The Vlasov seed leans on QRAM twice: it stores the externally-computed CDM force field (\(O(n_{\rm gr}^3)\) entries) in a QRAM queried at every time step, and it uses QRAM-backed state preparation for the Fermi–Dirac initial condition. This is the seed's single most load-bearing and least-physical assumption — exactly the target of the Jaques–Rattew critique. And the Carleman half maps the road not taken: had the neutrino force been self-consistent (nonlinear), \(R\lt1\) would have decided feasibility. Together these two primitives bound the whole enterprise — nonlinearity at one end, data-loading at the other — and honesty about both is what separates a real speedup claim from a hopeful one.
The thesis — one reduction, many dressings
Step back from seventeen papers and one pattern remains. Nearly every result in this neighborhood is the same move in different clothing — and seeing that move is the point of the whole survey.
(linearize / Carleman / discretize)→ Schrödinger / linear system→ Ham-sim or HHL→ amplitude estimation
read-out
Papers differ in which classical structure they exploit (linearity, sparsity, dissipation) and which quantum primitive they land on — but the skeleton is invariant. Learn it once and every paper becomes a variation you can place on the map.
The dressing determines the cost profile; the taxes determine whether any advantage survives. That is why the survey's spine is "one reduction, three taxes": the reduction tells you what a paper is, the taxes tell you whether it works.
Prediction from the map: the first genuinely useful quantum PDE result will be linear, QRAM-light, and read-out-light — a problem where the input is cheap to prepare and the wanted answer is a handful of numbers. The neutrino power-spectrum problem is a strong candidate if QRAM is solved. Strongly nonlinear / full-field problems (turbulence, full self-gravity) are the last, not the first, to fall.
- QRAM that doesn't erase the speedup — fault-tolerant data loading whose overhead doesn't cancel the advantage it enables (Jaques–Rattew).
- Nonlinearity beyond weak coupling — a method that survives \(R\ge\sqrt2\) (real turbulence, full self-gravity), not just Carleman's \(R\lt1\) regime.
- Read-out beyond low-dimensional observables — extracting more than a handful of numbers (e.g. the full field) without paying back the exponential.
- Fault-tolerant \(T\)-counts, not just logical resource estimates — the step Berntson's logical analysis still leaves open.
- Hardware past ~10 sites — carrying Sundar's demonstration (and our IonQ kernel) to a scale where the physics is non-trivial.
The comparison table — all 17, at a glance
The whole survey on one page. Facts (technique, maturity, speedup, result) are checked against each paper's source; the Q·R·N taxes are one consistent analytic overlay (as on the map). Click any paper to open its slide.
| Paper | Technique | Maturity | Q·R·N | Speedup claim | Key result |
|---|---|---|---|---|---|
| ★ the seed | |||||
| Miyamoto 2024 | Ham-sim + QAE | FTQC theory | QRN | exp. memory; Õ(n_gr+n_t) queries | neutrino P(k) from the linearized Vlasov equation |
| Family 1 · same process & technique (kinetic + Hamiltonian simulation) | |||||
| Engel–Smith–Parker 2019 | Ham-sim | FTQC theory | QRN | exp. in N_v, poly in t | first: linear Vlasov recast as Hamiltonian simulation |
| Ameri 2023 | Ham-sim (+ LCU) | FTQC theory | QRN | quadratic in system size | adds collisions; advantage is honestly only quadratic |
| Toyoizumi 2024 | QSVT | FTQC theory | QRN | query-optimal Ham-sim | QSVT linearized Vlasov–Poisson; Landau damping |
| Higuchi 2025 | QSVT + classical | FTQC + sim | QRN | O(N⁶T²) → poly, no CFL | hybrid QSVT+Maxwell; two-stream instability |
| Yamazaki 2023 | Ham-sim + Fourier | FTQC + sim | QRN | improved N_v scaling | closes self-gravity via measured density modes |
| Berntson 2026 | Carleman + block-enc. | FTQC (logical) | QRN | superquadratic in time | end-to-end logical estimate, weakly-nonlinear plasma |
| Family 2 · same target, different machinery | |||||
| Mocz & Szász 2021 | variational | NISQ (PoC) | QRN | exp. memory (log qubits) | proof-of-concept: fuzzy-DM Schrödinger–Poisson |
| Cappelli 2023 | variational (VarQTE) | NISQ (emul.) | QRN | log qubits vs resolution | emulated VarQTE on a Schrödinger–Poisson proxy |
| Engel 2021 | Carleman | FTQC theory | QRN | exp. only if nonlinearity weak | framework: Carleman embedding for nonlinear Vlasov |
| Vaszary 2024 ⚠ | Carleman | FTQC theory | QRN | none — polynomially worse | negative result: Carleman needs unphysical dissipation |
| May & Qin 2025 | second quantization | FTQC theory | QRN | encouraging scaling (∝ dim⁴) | second-quantization linearizes nonlinear Vlasov–Poisson |
| Sundar 2025 ⚛ | analog / spin model | ran on hardware | QRN | n/a — hardware demo | ran on a 9-qubit Rigetti chip: linear plasma waves |
| Family 3 · other PDEs, the same engine | |||||
| Costa–Jordan–Ostrander 2019 | Ham-sim | FTQC theory | QRN | exp. in grid size N | wave / Klein–Gordon / Maxwell become a Schrödinger equation |
| Gaitan 2020 | quantum ODE solver | FTQC theory | QRN | quadratic in ε (rough limit) | first quantum CFD: de Laval nozzle with a shock |
| Budinski 2021 | lattice-Boltzmann | PoC (simulator) | QRN | none proven (proof-of-concept) | quantum lattice-Boltzmann Navier–Stokes |
| Liu–Childs 2021 ★ | Carleman / QLSA | FTQC theory | QRN | exp. in t, only if R<1 | provable nonlinear-ODE solver; R≥√2 intractable |
One pattern to carry away: advantage is cleanest at the top (linear kinetics, nonlinearity tax = none) and most contested at the bottom (Carleman / nonlinear, where Vaszary finds no win) — and only two rows ever left a simulator.
Glossary — the four shared primitives
Every algorithm in the survey is built from these four building blocks. Learn their promises and their catches and you can read any paper on the map. Each links to its full slide.
| Primitive | What it provides | The catch |
|---|---|---|
| HHL 2009 | solves \(A\mathbf x=\mathbf b\) as a state \(|x\rangle\) in polylog(N) | read-out is expectation-only; cost grows with condition number \(\kappa\) and sparsity |
| Hamiltonian simulation | near-optimal \(e^{-iHt}\): \(O(\alpha t+\log 1/\epsilon)\) via qubitization | the cost is \(\alpha t\) (block-encoding norm), not \(\|H\|t\) — a looser bound |
| Amplitude estimation | a mean/amplitude to error \(\epsilon\) in \(O(1/\epsilon)\) | only quadratic over sampling; needs coherent \(A\) and \(A^\dagger\) |
| Carleman + QRAM | Carleman linearizes weak nonlinearity; QRAM loads classical data | \(R\ge\sqrt2\) provably intractable; QRAM's polylog query hides \(O(N)\) corrected cost |
The Miyamoto pipeline is exactly: Carleman/QRAM-free (it linearizes physically) → Hamiltonian simulation of \(e^{-iHt}\) → amplitude estimation to read \(P_\nu(k)\) — the HHL-lineage idea that a PDE is linear algebra, made unitary. The whole neighborhood is permutations of these four moves.