We report experimental digital quantum simulation of the one-dimensional Fermi-Hubbard model on a superconducting quantum processor at a scale beyond the reach of exact statevector simulation and challenging for state-of-the-art tensor-network methods. We encode this problem using up to 120 qubits through an efficient mapping that reduces circuit complexity, and we improve accuracy through error suppression to simulate dynamical evolution using up to 90 Trotter steps. From a vacancy defect introduced in the middle of an
L=31-site (62-qubit) Néel initial state, we directly observe spin-charge separation to
t=9 in natural units and quantitatively extract velocities
vc,
vs, which match classical simulations across model parameters, and analytics in appropriate regimes. We then extend experiments to
L=60 fermions (120 qubits) and evolution times to
t=6 using 30 Trotter steps. Quantum-processor outputs agree quantitatively with approximate classical simulations performed using a time-dependent variational principle (TDVP) solver and we observe that increasing the TDVP bond dimension through
χ=4096 expands the range of evolution times within which agreement has RMSE
∼1% before the approaches diverge. Owing to the large scale of the simulation and the use of overhead-free error-suppression techniques, for simulated evolution times at the limit of quantum/classical agreement (
t≳5 in natural hopping units), the wall-clock runtime of the quantum processor is nearly three orders of magnitude faster than TDVP simulations with
χ=4096 and implemented using the ITensor and TeNPy packages. These results establish contemporary digital quantum processors as a versatile, quantitatively accurate, and competitive platform for the study of fermionic many-body dynamics in regimes where leading classical methods can become prohibitively expensive.
Establishing quantum advantage requires comparison against the best achievable classical simulation. The Q-CTRL team recently simulated quench dynamics of the one-dimensional Fermi-Hubbard model on an IBM processor, completing a
L=60 evolution to time
t=6 in under three minutes and claiming a
3000× speedup over classical Time-Dependent Variational Principle (TDVP) simulation at bond dimension
χ=4096. Their classical benchmark required over 160 hours on a CPU cluster, failed to converge in the high-entanglement regime
t∈[5.2,6], and left the most challenging window of the experiment unverified. Here, we push the boundaries of classical simulation by exploiting the full
U(1)×SU(2) symmetry of the Fermi-Hubbard Hamiltonian combined with GPU-accelerated tensor contractions. Reaching bond dimensions up to
χ≈62,000 on four NVIDIA H200 GPUs -- among the largest ever achieved in TDVP simulations and fifteen times larger than Q-CTRL's classical baseline -- we achieve fully converged results across the entire simulation window, including rigorous certification of the previously unresolved high-entanglement regime
t∈[5.2,6]. We further advance the classical frontier to
t=7, which lies beyond the quantum hardware experiment and any previously verified classical evolution of the full wavefunction. At the bond dimension comparable to Q-CTRL's best classical run, our GPU implementation completes in
∼100 minutes, directly reducing the claimed
3000× quantum advantage to
∼36×. These results substantially narrow the quantum-classical performance gap and establish a new standard for tensor-network benchmarking of large-scale quantum simulations.
Classical simulations of quantum circuits play a vital role in the development of quantum computers and in taking the temperature of the field. Here, we classically simulate various physically-motivated circuits using 2D tensor network ansätze for the many-body wavefunction, which match the geometry of the underlying quantum processor. We employ a generalized boundary Matrix Product State contraction algorithm to controllably generate samples from the resultant tensor network states. Our approach systematically converges both the quality of the final state and the samples drawn from it to the true distribution defined by the circuit, with GPU hardware providing significant speedups over CPU hardware. Using these methods, we simulate the largest local unitary Jastrow ansatz circuit from recent IBM experiments to numerical precision. We also study a domain-wall quench in a 2D discrete-time Heisenberg model on large heavy-hex and rotated square lattices, reflecting IBM’s and Google’s latest processors, respectively. We observe a rapid buildup of complex loop correlations on the Google Willow geometry, which significantly impacts local properties of the system. Meanwhile, loop correlations build up extremely slowly on heavy-hex processors, with almost negligible impact even at large circuit depths. Our results underscore the role processor geometry plays in classical simulability.
Simulating the time-dynamics of quantum many-body systems was the original use of quantum computers proposed by Feynman, motivated by the critical role of quantum interactions between electrons in the properties of materials and molecules. Accurately simulating such systems remains one of the most promising applications of general-purpose digital quantum computers, in which all the parameters of the model can be programmed and any desired physical quantity output. However, performing such simulations on today's quantum computers at a scale beyond the reach of classical methods requires advances in the efficiency of simulation algorithms and error mitigation techniques. Here we demonstrate programmable digital quantum simulation of the dynamics of the 2D Fermi-Hubbard model -- one of the best-known simplified models of electrons in crystalline solids -- at a scale beyond exact classical state-vector simulation. We implement simulations of this model on lattice sizes up to
6×6 using 72 qubits on Google's Willow quantum processor, across a range of physical parameters, including different on-site electron-electron interaction strengths and magnetic flux values, and study phenomena including formation of magnetic polarons (charge carriers surrounded by local magnetic polarisation), dynamical symmetry-breaking in stripe-ordered states, attraction of charge carriers on an entangled background state known as a valence bond solid, and the approach to equilibrium through thermalisation. We validate our results against exact calculations in parameter regimes where these are feasible, and compare them to approximate classical simulations performed using tensor network and operator propagation methods. Our results demonstrate that meaningful programmable digital quantum simulation of many-body interacting electron models is now feasible on state-of-the-art quantum hardware.
A recent quantum simulation of observables of the kicked Ising model on 127 qubits implemented circuits that exceed the capabilities of exact classical simulation. We show that several approximate classical methods, based on sparse Pauli dynamics and tensor network algorithms, can simulate these observables orders of magnitude faster than the quantum experiment and can also be systematically converged beyond the experimental accuracy. Our most accurate technique combines a mixed Schrödinger and Heisenberg tensor network representation with the Bethe free entropy relation of belief propagation to compute expectation values with an effective wave function-operator sandwich bond dimension >16,000,000, achieving an absolute accuracy, without extrapolation, in the observables of <0.01, which is converged for many practical purposes. We thereby identify inaccuracies in the experimental extrapolations and suggest how future experiments can be implemented to increase the classical hardness.
The famous, yet unsolved, Fermi-Hubbard model for strongly-correlated electronic systems is a prominent target for quantum computers. However, accurately representing the Fermi-Hubbard ground state for large instances may be beyond the reach of near-term quantum hardware. Here we show experimentally that an efficient, low-depth variational quantum algorithm with few parameters can reproduce important qualitative features of medium-size instances of the Fermi-Hubbard model. We address 1 × 8 and 2 × 4 instances on 16 qubits on a superconducting quantum processor, substantially larger than previous work based on less scalable compression techniques, and going beyond the family of 1D Fermi-Hubbard instances, which are solvable classically. Consistent with predictions for the ground state, we observe the onset of the metal-insulator transition and Friedel oscillations in 1D, and antiferromagnetic order in both 1D and 2D. We use a variety of error-mitigation techniques, including symmetries of the Fermi-Hubbard model and a recently developed technique tailored to simulating fermionic systems. We also introduce a new variational optimisation algorithm based on iterative Bayesian updates of a local surrogate model.