Fermionic dynamics on a trapped-ion quantum computer beyond exact classical simulation

TL;DR

An efficient quantum simulation algorithm is implemented on Quantinuum's System Model H2 trapped-ion quantum computer for the time dynamics of a 56-qubit system that is too complex for exact classical simulation, and behaviour is found that differs from predictions made by classical tensor network methods.

Abstract

Simulation of the time-dynamics of fermionic many-body systems has long been predicted to be one of the key applications of quantum computers. Such simulations -- for which classical methods are often inaccurate -- are critical to advancing our knowledge and understanding of quantum chemistry and materials, underpinning a wide range of fields, from biochemistry to clean-energy technologies and chemical synthesis. However, the performance of all previous digital quantum simulations of fermions has been matched by classical methods, and it has thus far remained unclear whether near-term, intermediate-scale quantum hardware could offer any computational advantage in this area. Here, we implement an efficient quantum simulation algorithm on Quantinuum's System Model H2 trapped-ion quantum computer for the time dynamics of a 56-qubit system that is too complex for exact classical simulation. We focus on the periodic spinful 2D Fermi-Hubbard model and present evidence of spin-charge separation, where the elementary electron's charge and spin decouple. In the limited cases where ground truth is available through exact classical simulation, we find that it agrees with the results we obtain from the quantum device. Employing long-range Wilson operators to study deconfinement of the effective gauge field between spinons and the effective potential between charge carriers, we find behaviour that differs from predictions made by classical tensor network methods. Our results herald the use of quantum computing for simulating strongly correlated electronic systems beyond the capacity of classical computing.

Figures (40)

• Figure 1: Fermi-Hubbard lattice, initial state and map of a spin sector on the quantum circuit.Top, The Fermi-Hubbard model instance studied in this work is a double-periodic system of size $|\mathcal{L}|=L_x\times L_y=7\times 4$ (56 qubits), with $\Phi=\pi$ flux on the long direction. Around each highlighted plaquette $\phi_{ij}=\pi/4$, with $i$, $j$ ordered such that they form a directed loop (example on the right) . Each site (blue circles) can accommodate spin-up and spin-down electrons. Middle, The initial state is a dimerised configuration where we place a maximally entangled state, shown as a blue link. Each link corresponds to the $S_z^{\rm total}=0$ triplet state. The empty and full red circles correspond to a holon and a doublon, respectively. Bottom, Structure of the quantum circuit for a single Trotter step on qubits (blue circles) representing a single spin-sector. The grey line shows the structure of the Jordan-Wigner encoding. "Boundary hopping" means a hopping term across the vertical boundary, these are not caught by the swap network and must be implemented with a separate circuit, see \ref{['sec:circuits']}.
• Figure 2: Evolution of local charge, and spin correlations. Local charge density $\langle n_{i}(t)\rangle$ (squares as sites) and spin (connected) correlation function $C^{zz}_{ij}(t)$ between nearest-neighbours (represented by links) as a function of time for a,$U=0$ and b,$U=4$ for the error mitigated (TFLO + GPR) experimental data. In the charge sector, we see diffusion from the initial holon-doublon configuration towards the uniform state. In the free case (a), the charges develop a charge-density profile oscillating in the y direction, while in the interacting case (b), the charge profile is more disordered. In the spin sector, the initial triplet configuration takes $\sim t=0.7$ to melt, leaving behind a residual antiferromagnetic correlation, greater for the interacting case (see also \ref{['fig:doublons_and_triplets']}).
• Figure 3: Evolution of global doublon charge and magnetic correlations.a, Total number of doublons and b, triplet density as a function of time. While the initial state has a fixed energy density, it is not an equilibrium state. In particular, the number of doublons evolves nontrivially with time (panel a). By mitigating the raw data with TFLO + GPR, we recover a signal that achieves better agreement with the ground truth, here represented by Majorana propagation. For $U = 0$, we verified that the latter overlaps exactly with the FLO simulation of the Trotterised circuit (not shown). In both panels, the dotted line represents the exact, untrotterised FLO simulation of the time evolution for $U = 0$. The TFLO + GPR curve has been obtained by symmetrising over doublons and holons (since $N_{\mathrm{doublons}}(t) = N_{\mathrm{holons}}(t)$) via $N^{\mathrm{sym}}_{\mathrm{doublons}} =(N_{\mathrm{doublons}} + N_{\mathrm{holons}})/2$. Since $N_{\mathrm{doublons}}$ and $N_{\mathrm{holons}}$ are not independent, in computing the error bars of $N^{\mathrm{sym}}_{\mathrm{doublons}}$, we assumed perfectly correlated errors (i.e., $\sigma^{\mathrm{sym}}_{\mathrm{doublon}} = \sigma_{\mathrm{doublon}}+\sigma_{\mathrm{holon}}$). Note that for $N_{\rm doublons}$, the results from Majorana propagation are equally distant from the TDVP simulations as from the experimental data, but the experiment is capturing the late increase in doublon number. From panel b, we can see that while the initial melting of the antiferromagnetic order is slower for the interacting system, the spin ordering is lost essentially at $t\gtrsim 0.5$. In $n_{\rm triplets}$, the agreement between the mitigated experimental signal and the results from Majorana propagation is within the error bars, while the results from tensor network TDVP simulations converge to a smaller negative value, indicating more antiferromagnetic order in the tensor network simulations than expected for this state. The insets isolate the late-time dynamics of $n_{\rm triplets}$ for ease of comparison. In both panels, error bars indicate one standard deviation from the mean.
• Figure 4: Wilson open lines. a, In (lattice) gauge theory, the Wilson loops that serve as order parameters of the confinement/deconfinement transition are usually taken in space-time, here sketched as a green rectangle spanning a region of size $\ell$ in space and (Euclidean) time $T$. b, In our analysis, we fix the time slice and analyse the expectation of Wilson lines as a function of perimeter or area for different real times. c, Examples of some of the closed Wilson loops considered in the computation of $\mathcal{W}_d(\mathcal{C})$ in \ref{['fig:wilson_time_averaged']}a below for area $A = 10$ and varying perimeters $p$. d, Different paths in the computation of holon-doublon Wilson lines, located in opposite corners and represented by a white and red circle, respectively. (Green) This path has a full Néel order of the spinons in between. (Yellow) where the Néel order is broken. (Blue) Path where holons are present.
• Figure 5: Time averaged Wilson loops and lines. Time averaged Wilson loops as defined in \ref{['eq:Wilson_loops']} with: a, fixed area, varying perimeter; b, fixed perimeter, varying area. Here, we observe the scaling of space-like Wilson loops with the perimeter. This operator is the order parameter of the $U(1)$ gauge field mediating the interactions between spinons. The nontrivial scaling with perimeter signals deconfinement. c, Time average of the expectation of the (absolute value of the) open Wilson line $V_{hd}(M)$ as a function of the Manhattan distances $M$.The expectation value of this operator measures the effective interaction potential between a doublon and a holon. Increasing interaction with distance signals confinement. Both experimental and TDVP data have been smoothed using GPR. In all panels, square points (Raw) indicate the raw device signal, triangular points (TMPS) indicate device signal corrected by the training with MPS (TMPS) procedure described in \ref{['app:subsec:tflo']}, and the solid lines (GPR) are the curves produced by GPR applied to the triangular (a, b) or square (c) points, in both time and perimeter/area directions. TMPS data have been obtained by training with TDVP data up to $t=0.3$. Dashed lines are produced by sampling $10,000$ shots from the MPS obtained by TDVP with $\chi = 2048$.