Simulating a ring-like Hubbard system with a quantum computer

We develop a workflow to use current quantum computing hardware for solving quantum many-body problems, using the example of the fermionic Hubbard model. Concretely, we study a four-site Hubbard ring that exhibits a transition from a product state to an intrinsically interacting ground state as hopping amplitudes are changed. We locate this transition and solve for the ground-state energy with high quantitative accuracy using a variational quantum algorithm executed on an IBM quantum computer. Our results are enabled by a variational ansatz that takes full advantage of the maximal set of commuting ${\mathbb{Z}}_2$ symmetries of the problem and a Lanczos-inspired error mitigation algorithm. They are a benchmark on the way to exploiting near term quantum simulators for quantum many-body problems.

Figure 1

The four-site system is shown left in (a) alongside the energy levels for the nearest-neighbour and next-nearest-neighbour hoppings in (b). In (a) the (next-)nearest-neighbour hopping is visualized with an orange (blue) arrow and labeled with $t$ ($t^{\prime }$). In (b), the fermionic states of the four-site ring are labeled by $(\lambda ,\sigma )$. For half filling with only next-nearest-neighbour hopping $t^{\prime }\ne 0,t=0,E/t^{\prime }=2{\lambda }^2$ the ground state is nondegenerate as shown in blue, while it is degenerate for nearest-neighbour hopping only $t\ne 0,t^{\prime }=0,E/t=\lambda +{\lambda }^{*}$, where the occupied states are orange and half-occupied states are gray.

Figure 2

Example of a sequential adaptive $R_y{R}_z$ circuit with $n_{\mathrm{CZ}}=3$ controlled $Z$ gates (CZ) adapted to the five-qubit backend ibmq_ourense. Solving the four-site Hubbard ring requires only four qubits, hence, the eight initial rotations $R_y$, $R_z$ only act on the first four qubits leaving the fifth qubit $q_4$ idle. Three CZ gates are chosen from the set of all natively implemented controlled gates on the device. Each is followed by $R_y{R}_z$ rotations on the related qubits yielding in total 20 variational parameters $\mathit{\Theta }$, one angle $\mathrm{\Theta }$ for each rotation $R_y$ or $R_z$.

Figure 3

The results for the ground state energy estimates at $U=t/2$ are shown for the three symmetry subspaces: In orange for $B_1$, in blue for $A_1$, and in green for $E$. The lines denote the exact diagonalization ground-state energy. The error bars visualise the weighted least square LEMA result for each CZ sequence (see Appendix pp1-s5), with $n_{\mathrm{CZ}}=3$. The crosses mark the ground-state estimate $E_{\text{no-n}}^c$ using the non-noisy simulation for the state optimized on the quantum computer. The dots represent the minimum-weighted least-square LEMA result of the four CZ sequences. For comparison, the minimum-weighted least-square result of $E_{\mathrm{opt}}^c$ is shown with circles with error bars too small to be shown. For better readability, the data for each of the $n_c=4$ sequences is shown right to the minimum result they belong to.

Figure 4

The procedure to undo the tapering method is shown in this figure for the symmetry ${\mathit{\sigma }}_{\mathcal{B}}$ as an example. First, an ancilla qubit $q_{k_{\mathcal{B}}}$ in the state $\left|\right(1-s_{\mathcal{B}})/2\rangle$ is added to the qubits $q_0-q_3$, which were used for the calculations before. Then, for each $Z_k$ gate in the tensor product representation of the symmetry ${\mathit{\sigma }}_{\mathcal{B}}=Z_{q_{k_{\mathcal{B}}}}{Z}_{k(e,\downarrow )}$, a controlled $X$ gate is applied on $q_{k_{\mathcal{B}}}$ with $q_k$ as control. In this case, there is only one $Z_k$ gate, namely $Z_{k(e,\downarrow )}$, yielding one controlled $X$ gate as shown in the figure.

Figure 5

Result of the measurement of the momentum on the quantum computer ibmq_ourense. We show the probability of measuring $\lambda =-1$ for the ground state for all results shown in Fig. 3 in the symmetry subspaces with $s_{C_2}=1$, where $\lambda \in \{1,-1\}$. The value $|\langle E_0|\mathrm{\Psi }\left({\mathit{\Theta }}_{\mathrm{opt}}\right)\rangle {|}^2$ denotes the overlap of the ground state with the learned state of the VQE in a non-noisy calculation. This shows a correlation of the measurement outcome with the overlap, as it approaches 1 for $B_1$ and 0 for $A_1$ for an overlap of 1, as expected.

Figure 6

The results of the non-noisy simulations for different circuit depths are shown. The generation of the data is done as explained in the main text. To obtain stable ground-state energy estimates for each considered $t^{\prime }/t$ and symmetry subspaces 20 VQE runs are performed and the lowest value is taken as a final result. These are averaged over $t^{\prime }/t$ and the symmetry subspaces to obtain the data shown. With a vertical line the minimal excitation gap within the symmetry subspaces of all points $t^{\prime }/t$ considered is shown.

Figure 7

The results of the VQE optimization for the six-site ring at $U=1.5t$ are shown. The VQE optimization was performed for the two symmetry subspaces involved in the phase transition. In the top panel, the obtained energies are shown for the symmetry subspace $A_1$ (blue) and ${}^3{A}_2$ (orange) along with the spectrum for the lowest nine states. The bottom panel depicts the accuracy of the VQE results with respect to the exact diagonalization result in the respective symmetry subspace. The black lines represent the energy difference of the lowest excited states with respect to the ground state.

Figure 8

We show the results for the three discussed symmetry subspaces, $A_1$ (blue), $B_1$ (orange), and $E$ (green). For each symmetry subspace and each $t^{\prime }/t$ we generate four entangling gate sequences, which are shown for different values at the $x$ axis for graphical reasons. The data within each panel belong to one $t^{\prime }/t$. For each entangling gate sequences, we perform five measurements, which are shown in opaque colors. With the vertical lines, the result for the weighted least square is shown, as discussed in the main text. The small horizontal lines indicate the result, which we obtain by using the parameters $\mathit{\Theta }$ found with the quantum computer and perform a non-noisy simulation. The two long horizontal lines indicate the energy scaling. The lower one corresponds to the exact diagonalization value for the ground-state energy $E_0$ in the respective subspace. The upper of the two long horizontal lines indicates $E_0+0.1t$.

Figure 9

The optimization process for the data obtained on the quantum computer is shown. For each data point, given by $t^{\prime }/t$ and its eigenvalue to the mirror symmetry $s_{\mathcal{M}}$, the optimization process as well as the ground-state energy $E_0$ in gray obtained via exact diagonalization is shown. The rows correspond to results for the same symmetry subspace, in each column the results for a given $t^{\prime }/t$ are shown. For each data point we choose four CZ sequences, so that in each panel four optimization processes are shown.