Programmable digital quantum simulation of 2D Fermi-Hubbard dynamics using 72 superconducting qubits

TL;DR

This work demonstrates 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.

Abstract

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\times 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.

Figures (76)

• Figure 1: Selection of initial states considered and their representation on hardware. We use various initial states to probe different physical phenomena. A Néel-ordered background with one or two holons is used to explore magnetic polaron formation and holon attraction, respectively. The holon stripe is used to probe percolation. Two holons on a singlet background is used to investigate spin-charge separation. Finally, an ensemble of states with varying numbers of holons at random positions is used to explore thermalisation as a function of doping. Each spin sector of the model is separately mapped to a Jordan-Wigner representation in which qubits are ordered along a line. Time dynamics is simulated via fermionic swap networks along this line. The two representations are then arranged on the hardware in an interlocking pattern. We evolve up to time 1.2 in units of inverse hopping strength, and measure the time-evolved quantum state in the computational basis. We then compute expectation values of observables, and analyse signals based directly on the measured bit-strings. We perform approximate classical simulations to benchmark the results.
• Figure 2: Fermi-Hubbard dynamics phenomenology explored at the $5\times 5$ scale.a) Deviation from mean holon density (squares in the background), local spin (circles in the foreground) and nearest-neighbour spin correlations (links). We observe that a localisation of the excess holon density (purple squares) correlates with the strongest development of nearest-neighbour correlations as the interaction and flux increase. Note that for non-zero flux and on-site interaction strength $U$, a local spin moment is trapped in the middle, where a lack of neighbour correlations is observed. b) Illustration of paths used to compute the percolation fraction. c) Probability as a function of time for a stripe of holons to connect two opposite sides of the lattice, computed by: post-selecting on shots with less than four doublons and with negative spin-spin correlation across sites, and applying error-mitigation. While mean-field approaches Zaanen2001 assume that the main dynamics of a stripe is to oscillate around its initial position, we observe that the stripe quickly breaks apart, even for the largest interaction strengths considered ($U=8$). The largest tensor network simulation that we considered ($\chi=2048$) agrees with this result. d) Estimating a Wilson line between a holon (open circle) and a doublon (site with both up and down arrows): this is computed by finding shots with a holon and a doublon at Manhattan distance $M$, and computing $\sum_{m,n\in\gamma} S_n^z S^z_m$ for consecutive sites along a path $\gamma$. We sum over all possible paths of the same length. Here we show three paths for $M=4$ as dotted, dashed and bold lines. e) Time-snapshots of the Wilson line between a holon and a doublon as a function of their separation, generated by the spin background. We observe that the potential decreases as a function of time and increases as a function of distance up to a separation $\sim 6$, with the experimental results showing the strongest suppression at zero interaction and flux, as expected for the non-interacting system. The classical simulations seem to converge to a different value. The decrease of the Wilson line for larger distances appears to be a finite-size effect, as going to even larger system sizes we see a movement of the peak to larger separations (see \ref{['sec:spin_charge_u1']}).
• Figure 3: IPR results. Full, charge, and spin versions of the marginal inverse participation ratio (IPR) for an $m=4$-site region of the $5\times5$ system, constructed from different local basis states. The "full" marginal IPR uses the bitstrings directly. The charge IPR is obtained by grouping states with definite charge, while the spin IPR uses a basis grouped by definite spin labels (see \ref{['sec:ipr']} for a full discussion). Results are shown for initial states with one and with six randomly placed holes on top of a Néel background, each averaged over five configurations. The marginal IPR quantifies the delocalisation of the wavefunction in the corresponding local basis: it is equal to $1$ for a state localised on a single basis configuration, and it becomes exponentially small in the marginal size (here $m=4$) for a delocalised state. For an on-site uniform distribution, the corresponding baseline values are $(1/4)^m$ for the full IPR and $(3/8)^m$ for the charge and spin IPRs. We use the convergence to the estimate on the uniform distribution as a measure of thermalisation. In all cases the marginal IPR decays in time, and this decay becomes slower as $U$ increases. The spin IPR approaches its baseline value more rapidly than the charge IPR, indicating faster local equilibration of the spin sector compared to the charge sector.
• Figure 4: Dynamics of the staggered magnetisation at increasing scale. The staggered magnetisation is $M_s = \frac{1}{L_xL_y}\sum_i (-1)^i S^z_i$ where $S^z_j$ is the spin in the $z$ direction at site $j$, and $L_x$, $L_y$ are the system dimensions, ranging from $4\times 4$ to $6 \times 6$, with and without a magnetic flux $\phi$ passing through the plane of the model. Different interaction strengths $U$ are indicated by colour. Error-mitigated hardware data is given by points, with solid curves corresponding to smoothed fits. For the $4\times 4$ case, exact state vector simulations are possible at all values of $U$, indicated by the solid lines in that panel. For larger lattice sizes, exact classical simulations computed by FLO are shown with dashed lines for the non-interacting case ($U=0$), and Majorana propagation simulations are shown with dashed lines. At $U=0$ the Majorana propagation simulations are expected to be accurate simulations of the ideal quantum circuit output, so differences to the exact time dynamics are indicative of Trotter error. Similar plots illustrating performance with scaling can be found in \ref{['sec:further_hw_plots']}.
• Figure 5: Flux is represented on the model by attaching Peierls phases $\phi_{ij}$ to lattice edges $(i,j)$ and dressing the Fermi-Hubbard hopping terms as in \ref{['eq:FH_Ham']}, the phases are directional with $\phi_{ij}=-\phi_{ji}$. When a charged particle moves in a closed loop the state picks up a phase equal to the sum of phases for each edge crossed, as illustrated in the left hand figure. The total phase $\phi$ is proportional to the magnetic flux through the loop, analogous to the Aharonov-Bohm effect with the Peierls phases taking the role of the vector potential. Accordingly, a choice of phases $\phi_{ij}$ specifies a flux profile up to some gauge freedom, on the right we show the phases we use to induce a $\pi$ flux through every plaquette (blank edges indicate a trivial phase). The directionality of the $\pi$ phases for each edge is arbitrary as the effect is the same either way.