Source-linked AI summary
Strategies for solving the Fermi-Hubbard model on near-term quantum computers
Chris Cade, Lana Mineh, Ashley Montanaro, Stasja Stanisic
TL;DR
The paper studies how to make VQE practical for preparing Fermi-Hubbard ground states despite exponential classical-simulation costs and hardware-mapping overhead. It develops efficient HV and NP ansätze with swap-network implementations, then evaluates their depth and state-preparation performance on small systems. The results support high-fidelity ground-state preparation with relatively shallow circuits, while the conclusions remain based on limited lattice sizes and shallow numerical experiments.
Problem
The Hubbard model’s classical exact-solvers face exponentially growing subspaces, motivating variational quantum approaches to ground-state preparation.
Method
The paper develops efficient HV and NP ansätze using commuting Hamiltonian decompositions and fermionic swap networks for hardware-aware implementation.
Results
The study finds that its efficient variational circuits can prepare Hubbard-model ground states with high fidelity at relatively low depth, including hardware-mapped implementations.
Takeaways & Limitations
The results suggest that variational circuits scaling at roughly N layers may suffice for the studied Hubbard instances.
Takeaways & Limitations
The numerical evidence is limited to relatively shallow circuits and small lattice sizes, so extrapolation to larger hardware-accessible systems remains conditional.
Abstract
from arXiv · showhide
The Fermi-Hubbard model is of fundamental importance in condensed-matter physics, yet is extremely challenging to solve numerically. Finding the ground state of the Hubbard model using variational methods has been predicted to be one of the first applications of near-term quantum computers. Here we carry out a detailed analysis and optimisation of the complexity of variational quantum algorithms for finding the ground state of the Hubbard model, including costs associated with mapping to a real-world hardware platform. The depth complexities we find are substantially lower than previous work. We performed extensive numerical experiments for systems with up to 12 sites. The results suggest that the variational ansätze we used -- an efficient variant of the Hamiltonian Variational ansatz and a novel generalisation thereof -- will be able to find the ground state of the Hubbard model with high fidelity in relatively low quantum circuit depth. Our experiments include the effect of realistic measurements and depolarising noise. If our numerical results on small lattice sizes are representative of the somewhat larger lattices accessible to near-term quantum hardware, they suggest that optimising over quantum circuits with a gate depth less than a thousand could be sufficient to solve instances of the Hubbard model beyond the capacity of classical exact diagonalisation.
I. THE VARIATIONAL METHOD
The paper frames VQE as a hybrid approach for preparing Hubbard-model ground states and identifies encoding, ansatz, measurement, optimisation, and hardware implementation as key design choices.
- VQE uses a classical optimiser to adjust quantum circuits whose measured energy estimates approximate the Hamiltonian expectation value.
- The implementation must specify the variational ansatz, energy-measurement procedure, classical optimiser, and real-hardware circuit details.
- A. Fermionic encoding: Jordan-Wigner encoding maps each fermionic mode to one qubit without qubit-count overhead, but some interactions become long Pauli strings.
- A. Fermionic encoding: The snake-shaped Jordan-Wigner layout supports fermionic swap networks and low-overhead Hamiltonian measurements, despite differing from physical qubit connectivity.
- A. Fermionic encoding: Vertical hopping terms are hardest to implement because vertically adjacent physical sites may be separated in the Jordan-Wigner ordering.
B. Variational ans¨atze
The paper develops variational ansätze based on Hamiltonian decomposition and number-preserving gates, balancing efficient state preparation against optimisation flexibility.
- B. Variational ans¨atze: The HV ansatz starts from the non-interacting ground state and alternates evolutions associated with decomposed Hubbard-Hamiltonian components.
- B. Variational ans¨atze: The Hamiltonian is split into onsite, two vertical-hopping, and two horizontal-hopping components that contain commuting terms for efficient evolution.
- B. Variational ans¨atze: The HV ansatz uses at most 5 variational parameters per layer but requires nontrivial preparation of its non-interacting initial state.
- B. Variational ans¨atze: The NP ansatz permits simpler initial states and more freedom in gate choice, but its larger parameter count makes optimisation more challenging.
- B. Variational ans¨atze: The Number Preserving ansatz generalises HV by replacing hopping and onsite terms with more general number-preserving operators.
- B. Variational ans¨atze: The HV ansatz is the NP special case that preserves spin and constrains many parameters to be identical or zero.
C. Efficient implementation of HV and NP ans¨atze
The efficient EHV and NP implementations use fermionic swap networks to make otherwise nonlocal vertical interactions Jordan-Wigner-local, reducing ansatz circuit depth.
- C. Efficient implementation of HV and NP ans¨atze: Fermionic swap networks repeatedly permute columns so vertical hopping interactions become local in the Jordan-Wigner ordering while preserving fermionic parity.
- C. Efficient implementation of HV and NP ans¨atze: The EHV approach is designed to reduce the overhead caused by long Jordan-Wigner strings in vertical hopping terms.
- C. Efficient implementation of HV and NP ans¨atze: The circuit construction uses alternating column-swapping operators U_R and U_L, with vertical interactions implemented at newly Jordan-Wigner-adjacent positions.
- C. Efficient implementation of HV and NP ans¨atze: For unrestricted connectivity, all vertical hopping terms require depth 2n_x for even n_x and 2n_x+1 for odd n_x.
- C. Efficient implementation of HV and NP ans¨atze: Combining horizontal hopping operations with swap-network gates gives one ansatz-layer depth of 2n_x+1 for even n_x and 2n_x+2 for odd n_x.
- C. Efficient implementation of HV and NP ans¨atze: The efficient HV ansatz differs from standard HV by applying vertical hopping terms in a different order.
D. Measurement
The measurement scheme reduces Hubbard-energy estimation to at most five commuting measurement sets, while accounting for basis transformations, Jordan–Wigner strings, sampling cost, and optimizer requirements.
- Measurement grouping: The Hamiltonian terms can be grouped into at most five commuting sets: one onsite set and four hopping sets.The lattice structure permits parallel measurement of terms within each set.
- Onsite terms: Onsite energies are obtained from computational-basis measurements, where each site's contribution is the probability that both spin qubits are measured as 1.The mapped onsite operator is |11⟩⟨11| for the two qubits associated with a site.
- Horizontal hopping terms: A two-qubit unitary diagonalizes each horizontal hopping operator, whose expectation equals the probability of 01 minus the probability of 10.Adjacent overlapping hopping pairs cannot be transformed simultaneously, requiring two ansatz preparations when nx > 2.
- Vertical hopping terms: Vertical hopping measurements use the same basis transformation while separately accounting for Pauli-Z strings through parity measurements.Multiple vertical hopping pairs can be transformed simultaneously in selected groups.
- Sampling and optimization cost: An energy measurement requires five circuit evaluations, and an energy estimate uses m such measurements.The resource analysis assumes approximately 10^4 circuit evaluations per term for energy accuracy around 10^-2, motivating a budget near 10^5 function evaluations.
- Optimizer choice: With exact energy measurements, L-BFGS was usually effective, whereas realistic measurement noise required stochastic methods such as SPSA or coordinate descent.Deterministic NLopt optimizers often failed or found lower-quality minima under realistic measurements.
1. Simultaneous perturbation stochastic approximation
SPSA estimates gradients from random perturbations using two function evaluations, while the paper also generalizes exact coordinate-wise minimization to ansätze with integer-eigenvalue generators.
- SPSA: SPSA estimates a gradient along a random direction rather than computing the full gradient, using two function evaluations per gradient estimate.The perturbation vector has independent Bernoulli ±1 components.
- SPSA: The SPSA update uses ck = c/(k + 1)^γ and ak = a/(k+1+A)^α, with γ = 0.101, α = 0.602, A = 100, a = 0.15, and c = 0.2 in the experiments.The parameters were tuned by a joint parameter sweep; smaller values reduced convergence speed but improved accuracy.
- SPSA: A modified SPSA procedure uses two coarse runs with 10^2 and 10^3 measurements, followed by a finer run with 10^4 measurements and two averaged gradient evaluations.The three stages use a 10 : 3 : 1 ratio for their numbers of steps.
- Coordinate descent: Coordinate descent exactly minimizes the energy with respect to one parameter at a time, exploiting the trigonometric-polynomial dependence of energy on parametrized gates.For gates generated by integer-eigenvalue Hamiltonians, the paper gives a generalization applicable to the number-preserving and Hamiltonian Variational ansätze.
F. Handling noise
The VQE procedure faces both measurement-induced statistical noise and circuit errors, but the number-preserving ansätze support Hamming-weight-based error detection without additional qubits or execution time.
- Noise sources: VQE must contend with statistical measurement noise and errors in the quantum circuit.Statistical noise can be reduced by taking more measurements.
- Error detection: The Number Preserving and Hamiltonian Variational ansätze preserve fermionic occupation number, equivalent to Hamming weight after the Jordan–Wigner transform.A final state with support on different Hamming weights can therefore signal an error.
- Error detection: Hamming weight can be measured during the energy-measurement procedure without additional cost.Onsite and hopping measurements jointly provide the information needed to determine pair and total Hamming weights.
II. NUMERICAL VALIDATION
Numerical tests evaluated ansatz expressivity and optimisation behavior for Hubbard grids up to 12 sites. The EHV ansatz showed exponentially decreasing infidelity with depth, while NP generally achieved high fidelity at lower ansatz depth.
- Experimental setup: The experiments simulated VQE with HV, EHV, and NP ansätze, increasing circuit depths and grid sizes to test ground-state expressivity.The tests initially used exact energy measurements before considering realistic measurements and noise.
- Ground-state representation: At most 18 EHV layers achieved fidelity at least 0.99 for every tested grid with at most 12 sites.The corresponding tests used t = 1, U = 2 and L-BFGS optimisation.
- Ground-state representation: The EHV ansatz generally required fewer physical gates per layer than HV, offsetting HV’s lower layer count in most cases.For 1 × n grids, HV and EHV are equivalent.
- Ground-state representation: The NP ansatz generally required lower depth than HV and EHV, but its larger parameter count substantially increased runtime.For a 1 × 11 grid, NP required approximately 105 function evaluations and 16.5 hours, versus fewer than 9000 evaluations and 1.5 hours for EHV.
- Depth scaling: Infidelity decreased exponentially with EHV circuit depth for the largest tested 12-site grids.The 2 × 6 grid appeared more challenging than the 3 × 4 grid.
B. Optimisation with realistic measurements
The study tested CD and SPSA optimisation under realistic measurements and evaluated depolarising noise with error detection. Both optimisers reached relatively high fidelity, while error detection generally preserved or modestly improved results, though broader validation remains necessary.
- Optimisation with realistic measurements: CD and SPSA were compared on four representative grids under matched measurement budgets and circuit-evaluation limits.The limits were approximately 1.2 × 10^7 energy measurements and 6 × 10^7 circuit evaluations.
- Optimisation with realistic measurements: Both optimisers achieved relatively high fidelity despite realistic measurement noise, while SPSA performed noticeably better on 1 × 6 and 3 × 3 grids.Each energy measurement involved at most 10^4 circuit runs, corresponding to an estimated error of approximately 10^-2.
- Optimisation with realistic measurements: SPSA used substantially fewer energy measurements to reach high fidelity, partly because CD required 19 measurements per iteration on a 3 × 3 grid versus 2 for SPSA.The passage presents this as one possible reason for SPSA’s measurement efficiency.
- Optimisation with noisy quantum circuits: Across tested noise rates and grid sizes, depolarising errors did not significantly change final fidelity compared with noiseless results.The experiments used error rates 10^-3, 10^-4, and 10^-6 on 2 × 2, 1 × 6, and 2 × 3 grids.
- Optimisation with noisy quantum circuits: Error detection usually produced a small fidelity improvement and less erratic optimiser behavior, but its benefit on larger grids requires more detailed experiments.Detected-error runs were discarded and measurements continued until the intended number of valid measurements was obtained.
- Concluding remarks: A 5 × 5 instance was estimated at 50 qubits, around 25 layers, and approximate fully-connected two-qubit depth 325, under the ∼N-layer heuristic.This estimate includes the depth required to prepare the initial state.
- Concluding remarks: The study’s Hubbard benchmark has simpler structure than typical electronic-structure Hamiltonians, motivating comparable VQE complexity analyses for other systems.The authors identify such analyses as an important direction for future work.
- Concluding remarks: Choosing the optimal classical optimiser remains an important challenge for larger Hubbard instances and other many-body problems.The authors suggest evaluating combinations, modifications, adaptive optimisation, and machine-learning-based methods.
Appendix A: Alternative fermion encodings
The appendix compares fermion-to-qubit encodings and their circuit costs, emphasizing locality, auxiliary-qubit overhead, and parallelisation constraints. Under the stated assumptions, the Ball-Verstraete-Cirac encoding is substantially cheaper than the Bravyi-Kitaev superfast encoding.
- Ball-Verstraete-Cirac encoding: The Ball-Verstraete-Cirac encoding avoids long Jordan-Wigner Z strings by adding auxiliary modes, while preserving local qubit operators.Each fermionic mode is generally associated with an auxiliary mode, and vertical hopping terms use these auxiliaries.
- Ball-Verstraete-Cirac encoding: 17 is the total two-qubit depth for all time-evolution steps under the Ball-Verstraete-Cirac transformation.This assumes arbitrary two-qubit gates have depth 1 and no locality restrictions.
- Bravyi-Kitaev superfast encoding: The Bravyi-Kitaev superfast estimate is inflated by overlapping encoded hopping terms that must be split into separately implemented groups.The authors note that a carefully ordered parity implementation might reduce this overhead.
- Bravyi-Kitaev superfast encoding: 82 is the total two-qubit depth for all time-evolution steps under the Bravyi-Kitaev superfast encoding.The estimate includes six groups each for horizontal and vertical terms plus two onsite-term circuits of depth 5.
- Comparison: 29 would remain the estimated depth even if all commuting horizontal, vertical, and onsite groups could be implemented simultaneously.This hypothetical optimisation would still exceed the Ball-Verstraete-Cirac depth.
Appendix B: Implementation on hardware
The hardware implementation maps the ansatz onto constrained qubit layouts using swap networks and architecture-specific parallelisation. Circuit depth depends strongly on connectivity, with nearest-neighbour layouts cheaper than Google Sycamore for the same ansatz layer.
- Nearest-neighbour architecture: 4nx + 1 is the final circuit depth per ansatz layer on the described nearest-neighbour architecture.Horizontal and vertical Jordan-Wigner-adjacent terms are arranged to be physically adjacent, while onsite interactions add depth 1.
- Connectivity overhead: The hardware layout requires extra swaps because Jordan-Wigner-adjacent qubits are not physically adjacent on Sycamore.The ansatz operators UR, UL, and V are therefore split into multiple layers.
- Google Sycamore architecture: 6nx + 1 is the final circuit depth per layer on Google Sycamore for even nx.For odd nx, the depth becomes 6nx + 2 because vertical hopping terms cannot always be parallelised with UR.
- Google Sycamore architecture: Odd nx increases the Sycamore depth to 6nx + 2 per layer.The odd-width circuit shown for a 5 × 3 grid implements V URUL in depth 7 before the remaining ansatz-layer accounting.
- Gate complexity: 136 gates is the upper bound for the special 2 × 4 system with 2 layers, versus 1000 gates in prior work.The paper reports this as more than a factor-of-seven reduction.
Appendix C: The Number Preserving anstaz
The Number Preserving ansatz restricts evolution to a chosen occupation-number sector, making initial-state selection important. Spread-out computational-basis states can help, but pre-initialisation is inconsistent and certain sectors are inaccessible from computational-basis starts.
- Initial states: All gates in the Number Preserving ansatz conserve fermion number, restricting VQE to the selected occupation-number subspace.A computational-basis initial state can reduce initial-state complexity relative to starting from the non-interacting ground state, but may require more layers.
- Initial states: Spread-out fermion placements generally outperform top-corner placements at lower depths on the 3 × 3, six-fermion example.Only the ground-state initialisation reaches fidelity 0.99; the other tested placements reach about 0.96 at depth 5.
- Initial states: For the 3 × 3 grid, preparing the non-interacting ground state takes depth 8, while each Number Preserving layer takes depth 7.Thus, a computational-basis start does not appear advantageous in this example.
- Pre-initialising ansatz parameters: Pre-initialising parameters can improve one grid and worsen another, with the difference diminishing as depth increases.The procedure improves the 2 × 3 case at depth 2 but worsens the 3 × 3 case in the reported comparisons.
- Occupation-number limitation: A computational-basis start cannot reach the non-interacting ground state when one mode is occupied.All Hamming-weight-1 computational-basis states lie in the Hamiltonian’s null space, which hopping evolution preserves.
Appendix D: Simulation choices
The appendix tests how coupling strength, filling, degeneracy, and measurement sampling affect simulation results. Performance remains strong for moderate U, while stronger coupling and realistic measurements increase the required resources.
- Effect of U: For all tested U ≤4, the final fidelity remains quite high at depths calibrated for U = 2.Fidelity decreases as U increases because the ansatz is initialised from the U = 0 ground state.
- Effect of U: At U = 4, reaching 0.99 fidelity requires only 1 or 2 extra ansatz layers beyond U = 2, whereas U = 8 requires roughly twice as many layers.The required depth generally increases with U.
- Measurement requirements: Moving from U = 2 to U = 8 would likely require 16 times more measurements to maintain the same statistical uncertainty.The paper attributes this to energy terms scaled by U and an uncertainty that increases linearly with U.
- Half-filling: Half-filled instances require depths comparable to those for the full Hamiltonian at the same ansatz settings.The appendix also tracks infidelity, ground-state error, and double-occupancy error for a 2 × 4 half-filled example.
- Measurement requirements: The measured standard error decreases as 1/√m with the number of measurements m.The reported experiments motivate using 10^4 measurements per Hamiltonian term to estimate energy with accuracy about 10^-2.
Appendix E: Preparing the initial state of the non-interacting Hubbard Hamiltonian
The appendix compares several ways to prepare the non-interacting Hubbard ground state and implement its two-dimensional fermionic Fourier transform. The analyses show that parity-correction and swap-network constructions have different depth trade-offs across architectures and grid sizes.
- Naive implementation: The naive FFT separates horizontal and vertical transforms, but vertical parity corrections prevent parallel execution across columns.The horizontal transforms can run in parallel across rows; the vertical component requires nonlocal Z-string corrections.
- Naive implementation: For an n × n lattice, the naive FFT depth scales as Θ(n^3) when the one-dimensional FFT depth is O(n).Here TF(n) denotes the depth of the one-dimensional fermionic Fourier transform.
- Asymptotically efficient implementation: The ancilla-based construction implements the parity-recovery operator Γ through parity-basis conversion, ancilla movement, phase updates, and reversal operations.The construction uses CNOT, SWAP, and CZ gates, with separate depth contributions for its four stages.
- Asymptotically efficient implementation: The ancilla-based method uses O(N) gates and can reduce its parity-conversion depth to O(√ny) on nearest-neighbour or O(log ny) on fully connected architectures.These reductions would slightly lower the grid size at which this approach begins outperforming alternatives.
- Asymptotically efficient implementation: The ancilla-based implementation has unrestricted depth 2ny + 8nx − 2, while nearest-neighbour restrictions add overhead to the relevant stages.For example, the second stage has depth 4nx + ny − 2 on nearest-neighbour hardware but 4nx with arbitrary pair interactions.
- Depth comparison: The optimized FFT depth is slightly below the analytic prediction because stages overlap when combined.The one-dimensional transforms use TF(nx) = nx − 1 and TF(ny) = ny − 1 in the cited implementation.
3. Fermionic swap networks for the FFT
Fermionic swap networks avoid Jordan–Wigner parity corrections by moving qubits until interacting orbitals become adjacent. Their depth scales linearly with the total number of sites, although they are especially advantageous for small grids.
- Fermionic swap construction: FSWAP networks move qubits while preserving the Jordan–Wigner relative phases, eliminating the Z-strings required for vertical hopping terms.The network uses fermionic swaps to make otherwise nonlocal interactions nearest-neighbour interactions.
- Fermionic swap construction: The UL and UR swap pattern brings orbitals from adjacent rows together so that vertical hopping terms can be applied locally.The swaps follow the snake ordering used for the Jordan–Wigner encoding.
- Fermionic swap construction: The network repeats its swap sequence N/8 − 1 times, requiring 3N/2 layers of fermionic swaps overall.This ensures that all neighbouring orbitals become adjacent at least once.
- Depth comparison: The swap-network FFT has O(N) depth, compared with the O(√N) scaling reported for the ancilla-based approach, but is superior for small grid sizes.The comparison concerns circuit depth, and the small-grid advantage is explicitly reported for the swap-network approach.
a. Modified swap network
The modified swap network swaps entire columns in parallel, enabling vertical FFT interactions without Jordan–Wigner parity corrections. Its depth is favorable for small square grids relative to the asymptotically efficient ancilla-based construction.
- Modified swap network: The modified method repeatedly swaps entire columns with parallel FSWAP gates so that all vertical interactions become local in the Jordan–Wigner ordering.It alternates UL and UR rounds across columns.
- Modified swap network: After columns are rearranged, FFT rounds on odd- and even-numbered columns can be applied in parallel while the columns are bounced across the grid.The procedure continues until all ny − 1 FFT rounds have been applied to both column sets.
- Complexity: The modified swap-network depth is TF(nx) + 2nx · TF(ny).TF(n) is the depth of the one-dimensional FFT on n qubits.
- Comparison: The appendix compares a naive implementation, Jiang et al.’s asymptotically efficient method, and two swap-network constructions, including this modification.The comparison covers FFT circuit depth and direct Slater-determinant preparation by Givens rotations.
- Comparison: For square lattices, the modified swap network is more efficient when TF(n) < (20n − 4)/(2n − 1), a condition satisfied for n ≤ 11 when TF(n) = n − 1.The comparison is between the modified swap-network and asymptotically efficient implementations.
- Direct preparation: Direct Slater-determinant preparation with Givens rotations has depth nxny − 1 and is reported as more efficient than the compared FFT approaches.This method exploits known occupation numbers to reduce the number of Givens rotations.