Method to efficiently simulate the thermodynamic properties of the Fermi-Hubbard model on a quantum computer

Many phenomena of strongly correlated materials are encapsulated in the Fermi-Hubbard model whose thermodynamic properties can be computed from its grand-canonical potential. In general, there is no closed-form expression of the grand-canonical potential for lattices of more than one spatial dimension, but solutions can be numerically approximated using cluster methods. To model long-range effects such as order parameters, a powerful method to compute the cluster's Green's function consists of finding its self-energy through a variational principle. This allows the possibility of studying various phase transitions at finite temperature in the Fermi-Hubbard model. However, a classical cluster solver quickly hits an exponential wall in the memory (or computation time) required to store the computation variables. Here it is shown theoretically that the cluster solver can be mapped to a subroutine on a quantum computer whose quantum memory usage scales linearly with the number of orbitals in the simulated cluster and the number of measurements scales quadratically. A quantum computer with a few tens of qubits could therefore simulate the thermodynamic properties of complex fermionic lattices inaccessible to classical supercomputers.

Figure 1

Diagrammatic representation of the relation between the single-particle Green's function $\mathbf{G}$, the bare Green's function of the noninteracting lattice ${\mathbf{G}}_{\mathbf{0}\mathbf{t}}$, and the self-energy $\mathrm{\Sigma }$.

Figure 2

The Luttinger-Ward functional $\mathrm{\Phi }$ is the sum of all the two-body skeleton diagrams. The functional derivative with respect to $\mathbf{G}$ gives all the diagrams for the computation of $\mathrm{\Sigma }$. In the case where $U=0$, then $\mathrm{\Phi }\left[\mathbf{G}\right]=0$.

Figure 3

The essence of the VCA method is to remove the one-body links (denoted $t$) between small clusters (contained in $\mathbf{V}$) from the lattice $\mathrm{\Gamma }$ and consider only the reference lattice ${\mathrm{\Gamma }}^{\prime }$ whose Hamiltonian ${\mathcal{H}}^{\prime }$ is block diagonal in the Wannier basis and easier to solve than the complete problem $\mathcal{H}$. The reference system generates a manifold of trial self-energies ${\mathrm{\Sigma }}^{\prime }$ parametrized by single-particle parameters ${\mathbf{t}}^{\prime }$. The self-energy functional can be evaluated exactly on this manifold as the interaction part of the Hamiltonian (the $U\mathrm{s}$) is left unchanged. The solution become asymptotically exact as the clusters are made to include more sites.

Figure 4

Reduced Brillouin zone in the reciprocal lattice. The quasimomentum vector $\mathbf{k}$ belongs to the reciprocal lattice of $\mathrm{\Gamma }$, while the $\mathrm{K}$ component belongs to the reciprocal lattice of a single cluster. The $\tilde{\mathbf{k}}$ vector belongs to the reciprocal superlattice (hopping between clusters).

Figure 5

Circuit to simulate the time-dependent correlation function (66) of the cluster Hamiltonian (15). The first part, meant to generate a Gibbs state, is taken from [40]. Register $R$ is used in the modified phase-estimation scheme to prepare a rectangular state between the bath and the system contained in register $Q$. When the bath is traced out the system channel is left in a Gibbs state from which the different correlation functions can be read from the one-qubit register $P$. The size of register $Q$ depends on the number of orbitals in the simulated cluster (typically $n=2L_c$) and the bath size (which can be some constant factor larger than the system register). Register $R$ is used as a digital component and $q$ is therefore the size required for the desired floating point accuracy on reading $s_{*}$. Note that the numbers in the controlled gates of register $R$ denote the index of the qubit which is acting as the control.

Figure 6

Circuit to measure the correlation function (66) from an input Gibbs state. Register $S$ initially contains a given Gibbs state at inverse temperature $\beta$ and register $P$ is a single qubit initialized in the zero state. $P$ is put in a state superposition by applying a Hadamard gate $\mathbb{H}$ and then used to apply the controlled evolution sequence $O_{\mu \nu }\left(\tau \right)\equiv U_S^{†}\left(\tau \right){\sigma }_{\nu }{U}_S\left(\tau \right){\sigma }_{\mu }$ with $U_S\left(\tau \right)=e^{-i{\mathcal{H}}^{\prime }\tau }$ to the system channel. Finally, the state superposition is reversed by a last Hadamard gate and the measurement in repeated to obtain the probability $P\left(\mathcal{M}\right)$, which returns information on the cluster Green's function (30).

Figure 7

Measured probabilities for different $X_{\mu }$ and $Y_{\mu }$ at different times. The time axis $\tau$ is in units $t^{-1}$ of the hopping energy. In this case the cluster parameters are $L_c=2,t=1,U={\mathrm{\Delta }}^{\prime }={\mu }^{\prime }=0$, and $T=0.1$.

Figure 8

Nonzero correlation functions computed from the results of Fig. 7 with the time in units of $t^{-1}$. The function was regularized with an $e^{-\eta \tau }$ term to remove the fast oscillations of the Fourier transform arising from the finiteness of the time domain; $\eta =\frac{\pi }{50}$ was used in this case. The cluster parameters are given by $L_c=2,t=1,U={\mathrm{\Delta }}^{\prime }={\mu }^{\prime }=0$, and $T=0.1$.

Figure 9

Real and imaginary parts of the frequency-dependent Green's functions arising from the correlation functions measured in Fig. 8. The frequency axis $\omega$ is in units of the hopping energy $t$. The cluster parameters are $L_c=2,t=1,U={\mathrm{\Delta }}^{\prime }={\mu }^{\prime }=0$, and $T=0.1$.

Figure 10

Potthoff functional $\mathrm{\Omega }$ for different variational parameters ${\mu }^{\prime }$ and ${\mathrm{\Delta }}^{\prime }$ of a cluster of size $L_c=2$ with parameters $t=1,U=0,\mu =0$, and $T=1$. The cross marks the saddle point at $\left(\begin{array}{c}{\mu }_{*}^{\prime }\\ {\mathrm{\Delta }}_{*}^{\prime }\end{array}\right)=\left(\begin{array}{c}0.0046\\ 0\end{array}\right)$. All the variational parameters are in units of the hopping energy $t$.

Figure 11

Electron momentum-frequency distribution $A\left(k,\omega \right)$ for a lattice with parameters $t=1,U=0,\mu =0$, and $T=1$. The cluster used had a $L_c=2$ site and the saddle point is the same as in Fig. 10. The dashed line is at the chemical potential and frequencies are in units of the hopping energy $t$.

Figure 12

Electron momentum distribution $N\left(k\right)$ for different chemical potentials $\mu$ and temperatures $T$ with $U=0$. The solid lines are the results from the numerical simulation of the quantum algorithm using time steps of size $d\tau =0.02$ up to ${\tau }_{\max }=200$, while the dashed lines come from an imaginary frequency summation. The parameters $T$ and $\mu$ are in units of the hopping energy $t$.

Figure 13

Electron momentum-frequency distribution $A\left(k,\omega \right)$ for a lattice with parameters $t=1,U=4,\mu =-2$, and $T=0.1$. The cluster used had $L_c=2$ site and the saddle point is at $\left(\begin{array}{c}{\mu }_{*}^{\prime }\\ {\mathrm{\Delta }}_{*}^{\prime }\end{array}\right)=\left(\begin{array}{c}-2\\ 0\end{array}\right)$. The dashed line is at the chemical potential. The frequency axis is in units of $t$.

Figure 14

Detailed circuit to prepare an approximate Gibbs state ${\rho }_{QC}\approx {\rho }_{\mathrm{Gibbs}}$ following [40]. The simulated inverse temperature $\beta$ is related to the measurement of $s_{*}$ by Eq. (B10). The initial states of $R$ and $Q$ are taken to be the zero state ${\left|0\right.〉}^{\otimes \left(q+m+n\right)}$; then the Hadamard gate ${\mathbb{H}}^{\otimes q}$ is applied on $R$ and $Q$ is transformed (nonunitarily) to the fully mixed state $\frac{1}{2^{m+n}}{\mathbb{I}}^{\otimes \left(m+n\right)}$. Then $q$ controlled-$U$ operations are applied, where the notation $U_{\tau }=U^{2^{\tau }}$ and $U=e^{-i\frac{{\mathcal{H}}_0}{{‖{\mathcal{H}}_0‖}_{\infty }}}$, with ${\mathcal{H}}_0={\mathcal{H}}^{\prime }+{\mathcal{H}}_B$. An inverse quantum Fourier transform is applied on register $R$ and the string $s_{*}$ is read from the first $q$ qubits. Register $S$ is then left in a simulated Gibbs state ${\rho }_{QC}^S$.