Fast, accurate, high-resolution simulation of large-scale Fermi-Hubbard models on a digital quantum processor

Abstract

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 using up to 90 Trotter steps, and quantitatively extract velocity ratios $v_c/v_s$ which match classical simulations across a range of model parameters. We then extend experiments to $L=60$ (120 qubits) and long 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; increasing the TDVP bond dimension through $χ= 4096$ expands the range of evolution times within which agreement has RMSE $\sim 1\%$ before the approaches diverge. Owing to the large scale of the simulation and the use of efficient overhead-free error-suppression techniques, for simulated evolution times at the limit of quantum/classical agreement ($t\gtrsim 5$ in natural hopping units), the wall-clock runtime of the quantum processor is up to $3000\times$ faster than an optimized TDVP simulation using $χ= 4096$. 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.

Figures (26)

• Figure 1: Application-aware compilation for fermionic simulation. (a) The fermion-to-spin mapping represents a length-$L$ chain of sites with $L=8$ here (top row, purple circles) as a system of $2L$ qubits, one for each site/spin combination (bottom row, red and blue circles). We use a pair-interleaved ordering $\downarrow \uparrow \uparrow \downarrow \downarrow \uparrow \uparrow \downarrow \cdots$. The nearest-neighbor hopping interactions (indicated by the $\leftrightarrow$ arrows in the top row) map to "short-hop" and "long-hop" terms (bottom arrows). (b) The Trotterized time-evolution circuit with layers color-coded for simplicity. Red indicates state preparation. Purple indicates the single-qubit $R_Z$ for onsite and chemical potential. Blue indicates short hopping terms, onsite interactions are yellow, and $\mathrm{fSWAP}$ for spin exchange are green. Subsequent Trotter layers are mirrored to improve the Trotter error scaling of the overall circuit. (c) As a companion to (b), this panel tracks the degree of freedom encoded by a given qubit. The $\mathrm{fSWAP}$ layer exchanges adjacent spin-up and spin-down qubits. Alternating red and blue shading indicates qubit pairs used to encode a single fermion. (d) Embedding the qubitized spin chain into a heavy-hex coupling graph on an IBM quantum device. Qubits are colored according to which spin they encode---spin-up (red) or spin-down (blue). Qubits in use have blue shading around the chain. The snaking layout is chosen to avoid qubits with low gate fidelity (indicated by orange color coding) relative to the median.
• Figure 2: Spin-charge separation exhibited in the evolution of a central vacancy in a Néel initial state for $L=31$ over a range of repulsive couplings in each column ($U > 0$). All simulations run to a total time $t = 9\, t_h^{-1}$. Simulations in panels (a,d) use 45 Trotter steps, (b,e) use 60 Trotter steps, and (c,f) use 90 Trotter steps. (a--c) Heatmaps of the charge tracer correlator ${C^{\rm c}_i(t) := \langle n_{i,\uparrow}(t) + n_{i,\downarrow}(t)\rangle - \langle n_{i,\uparrow}(0) + n_{i,\downarrow}(0)\rangle}$. Time flows top-to-bottom, and the horizontal axis indicates site location on the one-dimensional chain. The wavefront boundary, as identified by the protocol detailed in Supplementary Material, Sec. \ref{['sec:wavefront']}, is shown as a dashed line, and the velocity of the wavefront is reported in each panel. (d--f) Heatmaps of the spin tracer observable, ${C^{\rm s}_i(t) := 4(\langle S_i^z(t) S_{i_*}^z(t)\rangle - \langle S_i^z(t) \rangle \langle S_{i_*}^z(t)\rangle})$, where $i_*$ is the site of the initial vacancy. Axes oriented as above. The same wavefront extraction has been applied as above. (g--i) Extracted wavefront velocities as a function of $U/t_h \in \{0, 2, 4, 6, 8, 10, 12, 14\}$. All TDVP simulations run on a lattice of size $L=23$ to a total time $t = 8.9\, t_h^{-1}$, using a time step size of $0.1\, t_h^{-1}$ and max bond dimension of $\chi=1024$. For the digital quantum simulations we use the following conditions: For $U/t_h \in \{0, 2, 4\}$ Trotter step size $\Delta t = 0.2\, t_h^{-1}$; for $U/t_h \in \{6, 8, 10\}$$\Delta t = 0.15\, t_h^{-1}$; for $U \in \{12, 14\}t^{-1}_h$$\Delta t = 0.1\, t_h^{-1}$. The error bars are computed using 2000 bootstrap samples and correspond to one standard deviation.
• Figure 3: Digital quantum simulation of the Fermi-Hubbard model for two initial states of an $L=60$ chain with interaction strength $U/t_h = -2$: a Néel state (left) and a Néel state with a central vacancy defect (right). (a) Heatmap showing the evolution of the per-site spin-up occupation $\langle n_{i,\uparrow} \rangle$. To enhance the visibility of small late-time occupations, the colormap applies the transformation $\sigma(\langle n_{i,\uparrow} \rangle)$, where $\sigma(x) = [1 + \tanh(k(x-1/2))/\tanh(k/2)]/2$. (b)--(d) Time evolution of the per-site spin occupations, $\langle n_{i,\uparrow} \rangle$ and $\langle n_{i,\downarrow} \rangle$, for selected sites. The site shown in each panel is indicated by the purple marker in the one-dimensional chain schematic shown in the inset. (e) Heatmap analogous to (a) for an initial state in which the central site ($i=29$) is a vacancy. The black line indicates the causal light cone. (f)--(h) Time evolution plots analogous to (b)--(d) for the same set of sites. Error bars are smaller than the markers and have been omitted for clarity.
• Figure 4: Quantum simulator outputs benchmarked against TDVP for a range of bond dimensions, $\chi \in \{64, 128, 256, 512, 1024, 2048, 4096\}$. (a) Time evolution of the occupation expectation value of a representative spin-orbital, $\langle n_{46,\uparrow}\rangle(t)$ for site $i=46$. Quantum hardware results are shown as purple circles; TDVP results are shown as solid lines, with color indicating bond dimension $\chi$. The light rectangular box indicates the region of zoomed-in dynamics treated in panel (b). (b) Zoom of (a), showing that larger-$\chi$ TDVP simulations track the quantum data to later times before diverging. (c) RMSE calculated between the quantum and TDVP occupation values across all 120 spin orbitals, shown for all tested values of $\chi$ at the discrete times for which quantum measurements exist (markers). The dashed horizontal line and green shading indicate the estimated noise floor, $\mathrm{RMSE}=0.01$, and vertical dotted lines mark the last step where each curve represents $<2\%$ error. Red shading at right for $t\gtrsim 5.2 t^{-1}_{h}$ indicates the zone in which correctness of quantum and classical simulations is indeterminate. (d) Wall-clock runtime required to simulate evolution to time $t$ using quantum hardware (purple lines and markers) and TDVP (blue lines and markers, shaded by $\chi$).
• Figure S1: Schematic quantum circuit implementing a single Trotter step $U_{\mathrm{step}}$ of Eq. \ref{['eq:trotter_step']}, together with the operation $\mathcal{P}$ relabeling fermionic modes at the end, shown for $L=3$ sites. The two layers of $R_Z$ rotations realize the one-qubit term $H_{1\mathrm{Q}}$ of Eq. \ref{['eq:mu_and_U_terms']}, the two layers of $R_{XX}$ and $R_{YY}$ rotations realize the "short hops" in the Hamiltonian $H_\mathrm{S}$ of Eq. \ref{['eq:H_S']}, and the layer of $R_{ZZ}$ rotations implements the two-qubit onsite interaction term $H_{U,2\mathrm{Q}}$ of Eq. \ref{['eq:mu_and_U_terms']}. The "long hops" $H_\mathrm{L}$ of Eq. \ref{['eq:H_L']} are converted into the "short hops" using the layer of fSWAPs, $\mathcal{F}$, and the final virtual permutation of the fermionic modes $\mathcal{P}$, denoted by the mismatch between the input and the output qubit wire labels.