Few-qubit quantum-classical simulation of strongly correlated lattice fermions

Abstract

We study a proof-of-principle example of the recently proposed hybrid quantum-classical simulation of strongly correlated fermion models in the thermodynamic limit. In a ‘two-site’ dynamical mean-field theory (DMFT) approach we reduce the Hubbard model to an effective impurity model subject to self-consistency conditions. The resulting minimal two-site representation of the non-linear hybrid setup involves four qubits implementing the impurity problem, plus an ancilla qubit on which all measurements are performed. We outline a possible implementation with superconducting circuits feasible with near-future technology.

1 Introduction

Using highly controllable quantum devices to study other quantum systems, i.e., quantum simulation [1–4], offers a means to tackle strongly correlated fermion models that are intractable on classical computers. This is vital for understanding complex quantum materials [5] with strong electronic correlations that exhibit a plethora of exciting physical phenomena of immediate technological interest. Examples of such effects include the Mott metal-insulator transition [6, 7], colossal magnetoresistance [8], and high-temperature superconductivity [9, 10].

Classical numerical methods have limited ability to study even significantly simplified toy models of strongly correlated fermions. For instance, exact diagonalization faces exponential scaling with the system size, while quantum Monte Carlo methods [11, 12] are often crippled by the infamous fermionic sign problem [13]. Tensor network methods [14–18] are powerful in one spatial dimension where they track strong correlations accurately. However, in higher dimensional systems, correlations tend to grow more quickly with system size, making these methods computationally challenging.

Another well-established approach to the study of strongly correlated fermionic lattice systems is dynamical mean-field theory (DMFT) [19]. It reduces the complexity of the original problem, e.g., the Hubbard model [20] in the thermodynamic limit, by mapping it onto a simpler impurity problem that is subject to a self-consistency condition relating its properties to those of the original model. Since an impurity problem is local, the mapping corresponds to neglecting spatial fluctuations. In the limit of infinite spatial dimensions this mapping is exact, but for finite dimensions it is an approximation. Nonetheless for lattice geometries with a large coordination number, self-consistently solving the impurity problem can yield an accurate approximate solution to the original Hubbard problem.

The ‘impurity’ itself consists of a single lattice site taken from the original problem, and so inherits on-site interactions from the Hubbard model. This impurity site is then immersed into a time-dependent, self-consistent mean-field with which it can dynamically exchange fermions. The mean-field thus attempts to model the rest of the lattice and by being dynamical can describe retardation phenomena. Overall the impurity problem can be represented by a Hamiltonian in which the interacting impurity site is coupled to a discrete set of non-interacting ‘bath’ sites. The bath sites represent the mean-field and if there is an infinite number of them then the self-consistency condition is guaranteed to be satisfied. However, in practical implementations only a finite number of bath sites are used, which restricts the frequency resolution of the bath so self-consistency condition can only be fulfilled approximately. Nevertheless, many strongly-correlated features, e.g., the Mott transition, are still be captured correctly [19]. For a study of different bath discretisation strategies in DMFT, see [21].

Although DMFT maps a Hubbard model to an impurity model this is still a non-trivial quantum many-body problem to solve because of the interactions at the impurity site. It is usually solved by classical numerical methods, e.g., specialised versions of those used to tackle the original problem, which attempt to keep track of the quantum correlations between impurity and bath sites. Again this limits the number of bath sites that can be treated accurately.

Here, we consider an alternative approach where the impurity problem is solved with a quantum simulator, thus avoiding many issues that are inherent to the classical methods. Quantum simulation of fermionic models has so far been mostly restricted to the analogue paradigm, especially with ultracold atoms in optical lattices [22]. Digital simulation approaches, akin to universal quantum simulators [23], have started to emerge in recent years, for example based on superconducting circuits [24–27]. Different quantum simulation schemes for the Hubbard model have been proposed [28–31]. The number of qubits in these digital simulators is, however, presently rather small. A direct implementation of the Hubbard model would suffer from severe finite size effects. It is nevertheless still possible for a digital quantum simulator with a restricted number of qubits to describe fermionic models directly in the thermodynamic limit when the DMFT approach is adopted.

To demonstrate this method we focus on the minimal incarnation of DMFT, the so-called ‘two-site’ DMFT [32], where the impurity model consists of one impurity site and only one bath site, both with local Hilbert space dimension four, subjected to two specially chosen self-consistency conditions. Since two-site DMFT considers only the smallest possible impurity model, the approach cannot match the accuracy of full DMFT, but it can still give a qualitatively correct description of the infinite-dimensional Hubbard model, and its simplicity makes it a good starting point before advancing to more accurate schemes. For explicit details of two-site DMFT and its features compared to full DMFT we refer to Ref. [32].

The two-site system corresponds to four qubits, two for the impurity site and two for the bath site, while a fifth, ancillary qubit is used for measurements. This number of qubits is readily available in current digital quantum simulator platforms, with IBM having made a five-qubit quantum processor available to the public [33]. A nine-qubit processor has already been demonstrated in superconducting circuits [26, 27, 34]. Trapped-ion technologies also allow for digital quantum simulations with up to six qubits [35, 36]. Being commensurate with current state-of-the-art technology is a further justification for studying this minimal model. Our scheme is readily generalisable to a larger number of qubits allowing for more accurate simulations and potentially offering an exponential speed-up over classical Hamiltonian-based DMFT methods [37]. For example, the number of multiqubit Mølmer-Sørensen gates scales only linearly with the number of bath sites, enabling efficient simulations [28, 29, 37, 38].

The self-consistency conditions are taken care of iteratively in a classical feedback loop, which thus completes the non-linear, hybrid quantum-classical device we introduce. Dynamical mean-field simulations have already been proposed for such hybrid devices [37, 39]. Quantum gates similar to the ones needed in the two-site scheme have been used in demonstrating digital quantum simulation of fermionic models with superconducting circuits [26, 38]. We thus focus on superconducting circuits as a candidate platform, although, e.g., trapped ions [28, 35, 40, 41] could also be considered.

This paper is organised as follows. In Section 2, we further elucidate the framework of DMFT applied to the Hubbard model in infinite dimensions. Section 3 introduces the two-site DMFT scheme in detail. Section 4 discusses the implementation of this two-site scheme with special attention to superconducting circuits. In Section 5, we show the results of our analysis. We end with a summary in Section 6 and give an outline of the single-qubit interferometry measurement scheme in the Appendix.

2 Hubbard model in infinite dimensions and dynamical mean-field theory

A standard model to describe strongly correlated electron systems in thermodynamic equilibrium is the Hubbard Hamiltonian

$$\begin{aligned} \hat{H}=-t\sum_{\langle j, k \rangle\sigma} \bigl( \hat {c}^{\dagger}_{j,\sigma} \hat{c}_{k,\sigma} + \mathrm{H.c.} \bigr)+ U \sum_{j} \hat{n}_{j,\downarrow} \hat{n}_{j,\uparrow}. \end{aligned}$$

(1)

In this model, electrons with spin projections $\sigma=\downarrow, \uparrow$ ‘hop’ between adjacent lattice sites with tunnelling energy t. This process is described in the first term, where $\langle j, k \rangle$ denotes the sum over all nearest-neighbour sites j and k, and $\hat{c}^{\dagger}_{j,\sigma}$ and $\hat{c}_{k,\sigma}$ denote the fermionic creation and annihilation operators, respectively. The electrons interact with on-site Coulomb repulsion $U>0$, described in the latter term by the product of the local number operators $\hat {n}_{j,\downarrow}=\hat{c}^{\dagger}_{j,\downarrow}\hat {c}_{j,\downarrow}$ and $\hat{n}_{j,\uparrow}=\hat{c}^{\dagger }_{j,\uparrow}\hat{c}_{j,\uparrow}$.

Here, we consider the paramagnetic Hubbard model in an infinite-dimensional Bethe lattice in the thermodynamic limit at zero temperature. This setup has very simple self-consistency relations, which makes it an ideal test-bed for a proof-of-principle demonstration of a hybrid quantum-classical scheme.

The DMFT approach [19] to solving this model consists in neglecting spatial fluctuations around a single lattice site and replacing the rest of the many-body lattice in the thermodynamic limit by a time-translation-invariant, self-consistent mean-field $\Delta(\tau-\tau')$ (or $\Delta(\omega)$ in the frequency domain), as illustrated in Figure 1(a). The isolated lattice site can dynamically exchange fermions with the mean-field at time instants $\tau'$ and τ. This allows one to include retardation effects that are important in the presence of strong correlations. In short, the dynamical mean-field approach reduces the complexity of the full Hubbard model to an effective single-site system which is a slightly more benign many-body problem to solve. In infinite dimensions, DMFT becomes exact as the irreducible self-energy of the lattice model becomes strictly local in space, $\Sigma_{{\mathrm{latt}},jk}(\omega)=\delta_{jk} \Sigma_{{\mathrm{latt}},jj}(\omega)$, and its skeleton diagrams agree with those of a single-site, or impurity, model [19].

Figure 1

Dynamical mean-field theory. (a) DMFT neglects spatial fluctuations around a single lattice site j and replaces the rest of the lattice with an effective mean-field $\Delta(\tau-\tau')$ with which the isolated site dynamically exchanges fermions, subject to the self-consistency condition $G^{R}_{\mathrm{imp}}(\omega)= G^{R}_{{\mathrm{latt}},jj}(\omega)$. Here, $G^{R}_{\mathrm{imp}}(\omega)$ is the impurity Green function and $G^{R}_{{\mathrm{latt}},jj}(\omega)$ is the local part of the lattice Green function. (b) In Hamiltonian-based DMFT methods, one considers an impurity model which describes the local part of the Hubbard model directly and represents the mean-field as a set of non-interacting bath sites that are connected to the central, interacting impurity site. (c) The minimal representation of DMFT involves the impurity site, with on-site interaction U and chemical potential μ, coupled via the hybridization energy V to only one bath site. The bath has on-site energy $\epsilon_{c}$ that corresponds to the mean-field $\Delta(\tau-\tau')$ and is subject to two self-consistency conditions.

The solution of the effective single-site, or impurity, problem also yields the solution of the infinite-dimensional Hubbard model due to the self-consistency condition. This leads to the retarded single-particle impurity Green function in the frequency domain being given by

$$\begin{aligned} G^{R}_{\mathrm{imp}}(\omega)=\frac{1}{\omega+\mu-\Delta(\omega)-\Sigma _{\mathrm{imp}}(\omega)}, \end{aligned}$$

(2)

where μ is the chemical potential, and $\Sigma_{\mathrm{imp}}(\omega )$ denotes the impurity self-energy. We set $\hbar=1$ throughout the paper. The impurity Green function describes the response of the many-body system after a localized removal or addition of a particle on the impurity site and is defined in the time domain and at zero temperature as

$$\begin{aligned} iG^{R}_{\mathrm{imp}}(\tau)=\theta(\tau) \bigl\langle \bigl\{ \hat{c}_{\sigma}(\tau ), \hat{c}^{\dagger}_{\sigma}(0) \bigr\} \bigr\rangle , \end{aligned}$$

(3)

where i is the imaginary unit, τ is real time, $\{\cdot,\cdot \}$ denotes the anticommutator, $\theta(\tau)$ is the Heaviside step function, and the average is computed in the ground-state $|GS\rangle$ of the impurity model. The fermionic creation and annihilation operators are given in the Heisenberg picture. In the paramagnetic phase the Green function is spin symmetric and we therefore only need to work out $G^{R}_{\mathrm{imp}}(\omega)$ for one spin configuration.

The initially unknown mean-field $\Delta(\omega)$ has to be chosen such that $G^{R}_{\mathrm{imp}}(\omega)$ matches the local part of the retarded lattice Green function $G^{R}_{{\mathrm{latt}},jj}(\omega)$, i.e.,

$$\begin{aligned} G^{R}_{\mathrm{imp}}(\omega)=G^{R}_{{\mathrm{latt}},jj}( \omega), \end{aligned}$$

(4)

where j is the (arbitrarily chosen) lattice site from which the removal or addition of a particle occurs in the translationally invariant lattice model. The DMFT self-consistency condition Eq. (4) implies

$$\begin{aligned} \Sigma_{\mathrm{imp}}(\omega)=\Sigma_{{\mathrm{latt}},jj}(\omega), \end{aligned}$$

(5)

i.e., the impurity self-energy matches the local self-energy of the Hubbard model in the infinite-dimensional Bethe lattice.

In the general case, the DMFT self-consistency loop is iterated as follows (see also Ref. [19]). (i) First, guess the local self-energy $\Sigma_{{\mathrm{latt}},jj}(\omega)$. (ii) The local lattice Green function can be computed as $G^{R}_{{\mathrm{latt}},jj} (\omega)=\int_{-\infty}^{\infty} d\epsilon \rho _{0}(\epsilon)/ [\omega+\mu-\epsilon-\Sigma_{{\mathrm{latt}},jj}(\omega) ]$, where $\rho_{0}(\epsilon)=\sqrt {4{t^{*}}^{2}-\epsilon^{2}}/2\pi{t^{*}}^{2}$ is the non-interacting density of states of a Bethe lattice. The constant ${t^{*}}$ emerges from the requirement that the Hubbard hopping needs to be scaled as $t \sim t^{*}/\sqrt{z}$ to avoid a diverging kinetic energy per lattice site in the limit of infinite coordination, $z \rightarrow\infty$ [19]. (iii) With Eqs. (4) and (5), we obtain $\Delta(\omega)$ from Eq. (2) and the impurity model is then defined. (iv) Compute the impurity Green function and obtain the impurity self-energy $\Sigma_{\mathrm{imp}}(\omega )$. There are several means to do this [19]. (v) Set $\Sigma^{\mathrm{new}}_{{\mathrm{latt}},jj}(\omega)=\Sigma_{\mathrm{imp}} (\omega)$. (vi) Check if the self-energy has converged. If not, go to step (ii) and repeat.

Once self-consistent, the solution of the impurity problem then gives access to local single-particle properties of the original lattice model. For example, the local lattice spectral function is given by

$$\begin{aligned} {A_{{\mathrm{latt}},jj}(\omega)=-\operatorname{Im}\bigl[G^{R}_{{\mathrm{latt}},jj}( \omega+i\eta )\bigr]/\pi=-\operatorname{Im}\bigl[G^{R}_{{\mathrm{imp}}}( \omega+i\eta)\bigr]/\pi}, \end{aligned}$$

where η is a positive infinitesimal.

In Hamiltonian-based impurity solvers, one parameterizes $\Delta (\omega)$ by a set of bath sites (see Figure 1(b)). For any finite number of bath sites, the self-consistency condition (4) can only be approximately satisfied and in the extreme ‘two-site’ DMFT it turns out to be more suitable to reformulate Eq. (4) in a manner specially focused on this minimal representation [32] (see Section 3). Note that two-site DMFT is only able to provide a qualitatively correct description of the Hubbard model even in infinite dimensions [32].

3 Quantum simulator based on two-site DMFT

In terms of the single-impurity Anderson model (SIAM), the smallest impurity problem involves one fermionic site corresponding to the impurity and only one fermionic site corresponding to the entire mean-field as described in the previous section. Since two qubits are needed to encode the local Hilbert space of a fermionic site, we only require four physical qubits to implement this representation in the lab. The SIAM Hamiltonian for only one bath site reads

$$\begin{aligned} \hat{H}_{\mathrm{SIAM}}=&U\hat{n}_{1\downarrow} \hat{n}_{1\uparrow }-\mu\sum_{\sigma}\hat{n}_{1\sigma}+ \sum_{\sigma} \epsilon_{c} \hat{c}^{\dagger}_{2\sigma}\hat{c}_{2\sigma}+ \sum_{\sigma} V \bigl(\hat{c}^{\dagger}_{1\sigma} \hat{c}_{2\sigma} + \mathrm {H.c.} \bigr). \end{aligned}$$

(6)

Here, U is the Hubbard interaction at the impurity site 1, and μ is the chemical potential that controls the electron filling in the grand canonical ensemble. Furthermore, $\epsilon_{c}$ and V describe the on-site energy of the non-interacting bath site 2 and hybridization between the impurity and the bath site, respectively, and give the mean-field as

$$\begin{aligned} \Delta(\omega)=\frac{V^{2}}{\omega-\epsilon_{c}}. \end{aligned}$$

(7)

See Figure 1(c) for illustration of the two-site SIAM. The parameters $\epsilon_{c}$ and V are initially unknown and they need to be determined iteratively such that two self-consistency conditions are satisfied. For details of the derivation and motivation of these conditions we refer to Ref. [32].

The first condition is that the electron filling $n_{\mathrm{imp}}$ of the impurity site and the filling $n=\langle n_{j\downarrow} \rangle+ \langle n_{j\uparrow} \rangle$ of the lattice model match, i.e.,

$$\begin{aligned} n_{\mathrm{imp}} \equiv n . \end{aligned}$$

(8)

The second self-consistency condition is given by

$$\begin{aligned} V^{2}=\mathcal{Z}M_{2}^{(0)}= \mathcal{Z} \int_{-\infty}^{\infty} d\epsilon \epsilon^{2} \rho_{0}(\epsilon)=\mathcal{Z} {t^{*}}^{2}, \end{aligned}$$

(9)

where quasiparticle weight reads

$$\begin{aligned} \mathcal{Z}= \biggl[1- \frac{d \operatorname{Re}[\Sigma_{\mathrm{imp}}(\omega+i\eta )]}{d\omega} \bigg|_{\omega=0} \biggr]^{-1}. \end{aligned}$$

(10)

In Eq. (9), $M_{2}^{(0)}$ is the second moment of the non-interacting density of states, and the final equality follows from the semicircular density of states of the Bethe lattice.

3.1 Two-site DMFT protocol

The hybrid quantum-classical device implementing two-site DMFT consists of a few-qubit digital quantum simulator in which the impurity Green function is measured and of a classical feedback loop in which the parameters of the two-site SIAM are updated. The two-site DMFT protocol is summarized in Figure 2 and proceeds as follows (see also Ref. [32]).

Figure 2

Non-linear hybrid quantum-classical scheme. A digital quantum simulator works in conjunction with a classical feedback loop to perform a proof-of-principle demonstration of a two-site DMFT calculation.

• 1. First fix U and μ to the desired values in the SIAM and set the unknown parameters $\epsilon_{c}$ and V equal to an initial guess.

• 2. Measure the interacting Green function $iG^{R}_{\mathrm {imp}}(\tau)$. This can be done using, e.g., single-qubit interferometry (see details in the Appendix).

• 3. After Fourier-transforming the impurity Green function, the impurity self-energy is obtained classically from the Dyson equation

  $$\begin{aligned} \Sigma_{\mathrm{imp}}(\omega) = {G}_{\mathrm{imp}}^{R(0)}( \omega )^{-1}-G^{R}_{\mathrm{imp}}(\omega)^{-1}. \end{aligned}$$

  (11)

  Here, the non-interacting impurity Green function is given by

  $$\begin{aligned} {G}_{\mathrm{imp}}^{R(0)}(\omega)^{-1}=\omega+ \mu-\Delta(\omega). \end{aligned}$$

  (12)

  From the derivative of the self-energy one obtains the quasiparticle weight $\mathcal{Z}$ which directly yields the updated hopping parameter V via Eq. (9). The update for $\epsilon_{c}$ is found by minimizing the difference $|n_{\mathrm{imp}}-n|$ [32].

• 4. Steps 2 and 3 need to be repeated until V and $\epsilon _{c}$ are self-consistent, and $n_{\mathrm{imp}}=n$.

The self-consistent Green function $G^{R}_{\mathrm{imp}}(\omega)$ and self-energy $\Sigma_{\mathrm{imp}}(\omega)$ thus obtained are used to calculate approximations to local single-particle properties of the Hubbard model. Note that for larger systems the two-site DMFT steps need to be replaced with the general DMFT self-consistency loop outlined in Section 2.

4 Quantum algorithm for the single-impurity Anderson model with superconducting circuits

Here, we consider the quantum gates of the digital quantum simulator part in Figure 2, with special focus on superconducting circuits as the platform of choice [26, 27, 38].

4.1 Jordan-Wigner transformation of the SIAM

To implement the two-site SIAM with qubits, the fermionic creation and annihilation operators need to be mapped onto tensor products of spin operators which then act on the qubits via quantum gates. In order to obtain as simple quantum gates as possible in Section 4.3 and in the Appendix, we consider an ordering of the qubits where the first two qubits encode the spin ↓ for both fermionic sites while the last two correspond to spin ↑. This is achieved via the Jordan-Wigner transformation given explicitly as

$$\begin{aligned} \hat{c}^{\dagger}_{1\downarrow} =&\hat{\sigma}_{1}^{-} = \frac {1}{2} \bigl(\hat{\sigma}_{1}^{x}-i\hat{ \sigma}_{1}^{y} \bigr) , \end{aligned}$$

(13)

$$\begin{aligned} \hat{c}^{\dagger}_{2\downarrow} =&\hat{\sigma}^{z}_{1} \hat{\sigma }_{2}^{-} = \frac{1}{2}\hat{\sigma}^{z}_{1} \bigl(\hat{\sigma }_{2}^{x}-i\hat{\sigma}_{2}^{y} \bigr), \end{aligned}$$

(14)

$$\begin{aligned} \hat{c}^{\dagger}_{1\uparrow} =&\hat{\sigma}^{z}_{1} \hat{\sigma}^{z}_{2} \hat{\sigma}_{3}^{-} = \frac{1}{2}\hat{\sigma}^{z}_{1} \hat{\sigma }^{z}_{2} \bigl(\hat{\sigma}_{3}^{x}-i \hat{\sigma}_{3}^{y} \bigr), \end{aligned}$$

(15)

$$\begin{aligned} \hat{c}^{\dagger}_{2\uparrow} =&\hat{\sigma}^{z}_{1} \hat{\sigma }^{z}_{2} \hat{\sigma}^{z}_{3} \hat{\sigma}_{4}^{-} = \frac{1}{2}\hat {\sigma}^{z}_{1} \hat{\sigma}^{z}_{2} \hat{\sigma}^{z}_{3} \bigl(\hat{\sigma }_{4}^{x}-i\hat{\sigma}_{4}^{y} \bigr) , \end{aligned}$$

(16)

and $\hat{c}_{j\sigma}= (\hat{c}_{j\sigma}^{\dagger} )^{\dagger}$. Here, $\hat{\sigma}^{x}_{l}$, $\hat{\sigma}^{y}_{l}$, and $\hat{\sigma}^{z}_{l}$ are spin-$\frac{1}{2}$ Pauli operators for qubit l. For larger systems the use of the Jordan-Wigner transformation, e.g., $\hat{c}^{\dagger}_{j\downarrow}= ( \prod_{p<2j-1}\hat{\sigma}^{z}_{p} ) \hat{\sigma}^{-}_{2j-1}$, $\hat {c}^{\dagger}_{j\uparrow}= ( \prod_{p<2j}\hat{\sigma}^{z}_{p} ) \hat{\sigma}^{-}_{2j}$, $\hat{c}_{j\sigma}= (\hat {c}_{j\sigma}^{\dagger} )^{\dagger}$, becomes important, as considered in Refs. [29, 37]. The hybridization terms in the SIAM, which now involve many spins, can be implemented efficiently and scalably with multiqubit Mølmer-Sørensen gates, the number of which scales only linearly with the number of bath sites [28, 29, 37, 38].

With the mappings in Eqs. (13)-(16), the hybridization terms in the SIAM described in Eq. (6) transform into

$$\begin{aligned} V \bigl(\hat{c}^{\dagger}_{1\downarrow}\hat{c}_{2\downarrow} + \mathrm{H.c.} \bigr)=\frac{V}{2} \bigl(\hat{\sigma}_{1}^{x} \hat {\sigma}_{2}^{x} + \hat{\sigma}_{1}^{y} \hat{\sigma}_{2}^{y} \bigr), \end{aligned}$$

(17)

and

$$\begin{aligned} V \bigl(\hat{c}^{\dagger}_{1\uparrow}\hat{c}_{2\uparrow} + \mathrm {H.c.} \bigr)=\frac{V}{2} \bigl(\hat{\sigma}_{3}^{x} \hat{\sigma }_{4}^{x} + \hat{\sigma}_{3}^{y} \hat{\sigma}_{4}^{y} \bigr). \end{aligned}$$

(18)

The number operators become

$$\begin{aligned}& \hat{n}_{1\downarrow} = \frac{1}{2} \bigl(\hat{I} - \hat{\sigma }_{1}^{z} \bigr), \end{aligned}$$

(19)

$$\begin{aligned}& \hat{n}_{2\downarrow} = \frac{1}{2} \bigl(\hat{I} - \hat{\sigma }_{2}^{z} \bigr), \end{aligned}$$

(20)

$$\begin{aligned}& \hat{n}_{1\uparrow} = \frac{1}{2} \bigl(\hat{I} - \hat{\sigma }_{3}^{z} \bigr), \end{aligned}$$

(21)

$$\begin{aligned}& \hat{n}_{2\uparrow} = \frac{1}{2} \bigl(\hat{I} - \hat{\sigma }_{4}^{z} \bigr), \end{aligned}$$

(22)

and thus the interaction term can be written as

$$\begin{aligned} U\hat{n}_{1\downarrow}\hat{n}_{1\uparrow}=\frac{U}{4}\bigl(\hat{\sigma }_{1}^{z} \hat{\sigma}_{3}^{z} - \hat{ \sigma}_{1}^{z} - \hat{\sigma}^{z}_{3} \bigr), \end{aligned}$$

(23)

up to a constant. The total Hamiltonian then reads

$$\begin{aligned} \hat{H}_{\mathrm{SIAM}}={}& \frac{U}{4} \bigl( \hat{ \sigma}^{z}_{1} \hat {\sigma}^{z}_{3} - \hat{\sigma}^{z}_{1} - \hat{\sigma}^{z}_{3} \bigr) + \frac{\mu}{2} \bigl( \hat{\sigma}^{z}_{1} + \hat{\sigma}^{z}_{3} \bigr) - \frac{\epsilon_{c}}{2} \bigl( \hat{ \sigma}^{z}_{2} + \hat{\sigma }^{z}_{4} \bigr) \\ &{}+ \frac{V}{2} \bigl( \hat{\sigma}^{x}_{1}\hat{ \sigma}^{x}_{2} + \hat {\sigma}^{y}_{1} \hat{\sigma}^{y}_{2} + \hat{\sigma}^{x}_{3} \hat{\sigma}^{x}_{4} + \hat{\sigma}^{y}_{3} \hat{\sigma}^{y}_{4} \bigr), \end{aligned}$$

(24)

where we have dropped constant terms.

4.2 Quantum gates in superconducting circuits

We now consider how the Jordan-Wigner transformed SIAM in Eq. (24) can be implemented in an experimental arrangement based on superconducting circuits. We present two alternative approaches. The first one couples the qubits with a transmission line resonator, which leads to the so-called XY gate between the qubits. The second approach is the Controlled-Z _ϕ (CZ _ϕ ) gate, which can be obtained via a capacitive coupling of nearest-neighbour transmon qubits without using a resonator. These CZ _ϕ gates have been implemented with high fidelities of above 99% for a variant of transmon qubits called ‘Xmon’ qubits [42].

XY gates with resonators

The basic Hamiltonian coupling a set of qubits to the resonator has the form of a detuned Jaynes-Cummings model. By adiabatically eliminating the resonator one obtains, when the resonator is in the vacuum state, the well-known XY model for a pair of qubits l and m as

$$ \hat{H}_{XY}=\frac{g_{l} g_{m}}{2\Delta}\bigl(\hat{\sigma}^{x}_{l} \hat{\sigma }^{x}_{m}+\hat{\sigma}^{y}_{l} \hat{\sigma}^{y}_{m}\bigr). $$

(25)

Here, Δ is the detuning between the qubit level-spacing and the resonator mode, and $g_{l}$ is the coupling constant between qubit l and the resonator. The XY gate is universal for quantum computation and simulation in combination with single qubit gates, and is the natural interaction customarily employed in superconducting circuits.

CZ _ϕ gates with capacitive couplings

To perform the CZ _ϕ gate, one qubit is kept at a fixed frequency while the other carries out an adiabatic trajectory near an appropriate resonance of the two-qubit states. By varying the amplitude of this trajectory one can tune the conditional phase ϕ. The unitary for the CZ _ϕ is given by

$$\begin{aligned} {\mathrm{CZ}}_{\phi}= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & e^{i\phi} \end{pmatrix}. \end{aligned}$$

(26)

4.3 Quantum gate decomposition of the time-evolution operator

In order to use quantum gates for time-evolution, we utilize a Trotter decomposition of the time-evolution operator corresponding to $\hat {H}_{\mathrm{SIAM}}$ in Eq. (24). The first order Trotter expansion is given by

$$\begin{aligned} \hat{U}(\tau)={}&e^{-i\hat{H}_{\mathrm{SIAM}} \tau} \approx \bigl(e^{-i \frac{V}{2} (\hat{\sigma}^{x}_{1} \hat{\sigma}^{x}_{2} + \hat{\sigma }^{y}_{1} \hat{\sigma}^{y}_{2}) \frac{\tau}{N}} e^{-i \frac{V}{2} (\hat{\sigma}^{x}_{3} \hat{\sigma}^{x}_{4} + \hat{\sigma}^{y}_{3} \hat {\sigma}^{y}_{4}) \frac{\tau}{N}} e^{-i \frac{U}{4} \hat{\sigma }^{z}_{1} \hat{\sigma}^{z}_{3} \frac{\tau}{N}} \\ &{} \times\ e^{i (\frac{U}{4}-\frac{\mu}{2} ) \hat{\sigma }^{z}_{1} \frac{\tau}{N}} e^{i (\frac{U}{4}-\frac{\mu}{2} ) \hat{\sigma}^{z}_{3} \frac{\tau}{N}} e^{i \frac{\epsilon_{c}}{2} \hat {\sigma}^{z}_{2} \frac{\tau}{N}} e^{i \frac{\epsilon_{c}}{2} \hat {\sigma}^{z}_{4} \frac{\tau}{N}} \bigr)^{N}. \end{aligned}$$

(27)

Here, N is the number of Trotter (i.e., time) steps and $\frac{\tau }{N}$ is the size of the time step. In what follows, we use the two alternative approaches for quantum gates outlined in Section 4.2 to implement Eq. (27).

XY gates

As shown in Section 4.2, the XY gate, given by the expression $XY = {\exp [-i\frac{V}{2}(\hat{\sigma}^{x}_{l} \hat {\sigma}^{x}_{m} + \hat{\sigma}^{y}_{l} \hat{\sigma}^{y}_{m}) \frac {\tau}{n} ]}$, naturally appears when considering the use of a resonator quantum bus [43]. The quantum circuit for a single Trotter step with these gates is shown in Figure 3(a).

Figure 3

Quantum gates for one Trotter step. A single Trotter step is shown for (a) the XY method and (b) the CZ _ϕ method. Here, B is the entangling gate $B = \exp (-i\frac{U}{4}\hat{\sigma}^{z}_{1} \hat{\sigma}^{z}_{3} \frac{\tau}{N} )$, A is a two-qubit gate given by ${A = \exp (-i\frac{V}{2}\hat{\sigma}^{z}_{l} \hat{\sigma}^{z}_{m} \frac{\tau}{N} )}$, acting on qubits l and m, and the quantum gates C and D are single qubit $\hat{\sigma}^{z}$-gates, given by $C = \exp [i (\frac{U}{4}-\frac{\mu}{2} )\hat{\sigma}^{z}_{l} \frac{\tau}{N} ]$ and $D = \exp (i\frac{\epsilon_{c}}{2} \hat{\sigma}^{z}_{l}\frac{\tau}{N} )$, acting on qubit l. Finally, $X_{\phi}$ and $Y_{\phi}$ are ϕ-rotations along the x and y axis, respectively.

CZ _ϕ gates

To be able to utilize the CZ _ϕ gates, we write the time-evolution operator in Eq. (27) in terms of $\hat{\sigma}^{z}_{l} \hat{\sigma}^{z}_{m}$ (ZZ) interactions, taking into account that

$$\begin{aligned} \hat{\sigma}^{x}_{l} \hat{\sigma}^{x}_{m} = \mathcal{R}^{(l)}_{y} \bigl(\tfrac{\pi}{2}\bigr) \hat{ \sigma}^{z}_{l} \mathcal{R}^{(l)}_{y} \bigl(-\tfrac {\pi}{2}\bigr) \mathcal{R}^{(m)}_{y} \bigl( \tfrac{\pi}{2}\bigr) \hat{\sigma }^{z}_{m} \mathcal{R}^{(m)}_{y} \bigl(-\tfrac{\pi}{2}\bigr), \end{aligned}$$

(28)

and

$$\begin{aligned} \hat{\sigma}^{y}_{l} \hat{\sigma}^{y}_{m} = \mathcal{R}^{(l)}_{x} \bigl(-\tfrac{\pi}{2}\bigr) \hat{ \sigma}^{z}_{l} \mathcal{R}^{(l)}_{x} \bigl(\tfrac {\pi}{2}\bigr) \mathcal{R}^{(m)}_{x} \bigl(- \tfrac{\pi}{2}\bigr) \hat{\sigma }^{z}_{m} \mathcal{R}^{(m)}_{x} \bigl(\tfrac{\pi}{2}\bigr), \end{aligned}$$

(29)

where $\mathcal{R}^{(l)}_{\alpha} (\theta) = \exp(-i\frac{\theta }{2} \hat{\sigma}_{l}^{\alpha})$ is the rotation along the α-axis of qubit l. Note that in the computational basis, one can write, e.g.,

$$\begin{aligned} \exp \biggl(-i \frac{\phi}{2} \hat{\sigma}^{z}_{1} \hat{\sigma}^{z}_{2} \biggr)= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & e^{i\phi} & 0 & 0 \\ 0 & 0 & e^{i\phi} & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}, \end{aligned}$$

(30)

where we have neglected global phases. Thus, we have the decomposition

$$\begin{aligned} &\exp \biggl(-i \frac{\phi}{2} \hat{\sigma}^{z}_{1} \hat{\sigma }^{z}_{2} \biggr)=\mathcal{R}^{(1)}_{x}( \pi) {\mathrm{CZ}}_{\phi} \mathcal {R}^{(1)}_{x}(\pi) \mathcal{R}^{(2)}_{x}(\pi) {\mathrm{CZ}}_{\phi} \mathcal{R}^{(2)}_{x}(\pi), \end{aligned}$$

(31)

where the tunable ${\mathrm{CZ}}_{\phi}$-gate is given by Eq. (26).

The time-evolution operator in Eq. (27) in terms of ZZ interactions is given by

$$\begin{aligned} \hat{U}(\tau)={}&e^{-i\hat{H}_{\mathrm{SIAM}} \tau} \approx \bigl(\mathcal{R}^{(1234)}_{y} \bigl(\tfrac{\pi}{2}\bigr) e^{-i \frac{V}{2} \hat {\sigma}^{z}_{1} \hat{\sigma}^{z}_{2} \frac{\tau}{N}} e^{-i \frac {V}{2} \hat{\sigma}^{z}_{3} \hat{\sigma}^{z}_{4} \frac{\tau}{N}} \mathcal{R}^{(1234)}_{y}\bigl(-\tfrac{\pi}{2}\bigr) \\ &{} \times\mathcal{R}^{(1234)}_{x}\bigl(-\tfrac{\pi}{2}\bigr) e^{-i \frac{V}{2} \hat{\sigma}^{z}_{1} \hat{\sigma}^{z}_{2} \frac{\tau}{N}} e^{-i \frac {V}{2} \hat{\sigma}^{z}_{3} \hat{\sigma}^{z}_{4} \frac{\tau}{N}} \mathcal{R}^{(1234)}_{x}\bigl( \tfrac{\pi}{2}\bigr) \\ &{} \times\ e^{-i \frac{U}{4} \hat{\sigma}^{z}_{1} \hat{\sigma}^{z}_{3} \frac{\tau}{N}} e^{i (\frac{U}{4}-\frac{\mu}{2} ) \hat {\sigma}^{z}_{1} \frac{\tau}{N}} e^{i (\frac{U}{4}-\frac{\mu }{2} ) \hat{\sigma}^{z}_{3} \frac{\tau}{N}} e^{i \frac{\epsilon _{c}}{2} \hat{\sigma}^{z}_{2} \frac{\tau}{N}} e^{i \frac{\epsilon _{c}}{2} \hat{\sigma}^{z}_{4} \frac{\tau}{N}} \bigr)^{N}, \end{aligned}$$

(32)

where $\mathcal{R}^{(1234)}_{\alpha}(\phi)=\mathcal{R}^{(1)}_{\alpha } (\phi)\mathcal{R}^{(2)}_{\alpha} (\phi) \mathcal {R}^{(3)}_{\alpha} (\phi)\mathcal{R}^{(4)}_{\alpha}(\phi)$. The sequence of gates for one Trotter step is depicted in Figure 3(b).

A single Trotter step contains 5 ZZ two-qubit gates (corresponding to the A and B gates in Figure 3(b)) between nearest-neighbour qubits, 2 SWAP gates (for the B gate which acts on qubits 1 and 3), and 20 single-qubit rotations. We note that a SWAP-gate amounts to three CZ _ϕ gates, and a ZZ-gate amounts to two CZ _ϕ gates (see Eq. (31)). This number can be optimised further if we consider different orderings for odd and even Trotter steps as in Figure 4, such that subsequent gates may be suppressed. This reorganisation of interactions does not in principle affect the Trotter error. Hence, for a pair of Trotter steps, the number of gates is reduced, and we may only consider 10 ZZ two-qubit gates between nearest-neighbour qubits, 4 SWAP gates, and 32 single-qubit rotations.

Figure 4

Reordering of quantum gates. The ordering of gates shown for (a) an odd Trotter step and (b) an even Trotter step in the CZ _ϕ method. The gates depicted in red can be omitted as they cancel out during a sequence of time steps.

5 Results

We focus on the half-filled case, i.e., $\mu=\frac{U}{2}$ and $\epsilon_{c}=0$, which requires the least amount of quantum gates, since the C and D gates in Section 4 vanish. Note that since the value of $\epsilon_{c}$ is fixed in this case, it need not be updated in the self-consistency loop. We use $t^{*}$, the Hubbard hopping in infinite dimensions, as our unit of energy, hence time τ is measured in units of $1/t^{*}$. Note that τ refers here to the time in the evolution operator $\hat{U}(\tau)$, not to the actual time to run the experiment.

We show in Figure 5 the state fidelities $\mathcal{F}=| \langle\Psi(\tau) | \Psi_{T}(\tau) \rangle|^{2}$, where $|\Psi(\tau )\rangle$ denotes the state obtained with exact time-evolution using the full, non-Trotterized operator $\hat{U}(\tau)=\exp(-i\tau\hat {H}_{\mathrm{SIAM}})$ corresponding to the two-site SIAM in Eq. (6), and $| \Psi_{T}(\tau) \rangle$ is the state evolved using either the XY or CZ _ϕ quantum gates, for various Trotter steps N up to time $\tau=6/t^{*}$. Note that the number of qubits corresponding to the two-site SIAM is fixed, leaving only N as the parameter to be varied for increased accuracy. We use the initial state $|\Psi(\tau=0)\rangle=\hat{c}^{\dagger}_{1\downarrow}|GS\rangle /||\hat{c}^{\dagger}_{1\downarrow}|GS\rangle||$, where $|GS\rangle$ is the ground-state of the two-site SIAM in Eq. (6), which is a relevant state for obtaining the impurity Green function at zero temperature (see Eq. (3)). As expected, using XY gates displays superior fidelities, since CZ _ϕ gates require an extra factorization of the hybridization term (see Section 4). For $N=24$ steps, the state fidelity using XY gates remains over 99% throughout the evolution. In what follows, we use only XY gates for the time-evolution for concreteness.

Figure 5

Time-evolution of the state fidelity. State fidelities $\mathcal{F}=| \langle \Psi(\tau) | \Psi_{T}(\tau) \rangle |^{2}$ using the XY method (blue diamonds) and CZ _ϕ gates (red stars) obtained with (a) 6, (b) 12, (c) 18, and (d) 24 Trotter steps up to time $\tau=6/t^{*}$. We set $U=4t^{*}$ and $V=t^{*}$.

As shown in Section 3, the main object of interest is the retarded impurity Green function. One possibility to measure $iG^{R}_{\mathrm{imp}}(\tau)$ is single-qubit interferometry (see the Appendix for details), which raises the total number of qubits in the experimental arrangement to five. In Figure 6 we plot the impurity Green function obtained from evolving the state with XY gates compared to exact evolution of the two-site SIAM for different N. We see that the Green function from the XY approach starts to follow the curve of the exact Green function better for increasing N. In our subsequent analysis, we use $N=24$ up to $\tau=6/t^{*}$ to study what two-site DMFT physics can be captured with the digital approach.

Figure 6

Impurity Green function in the time domain. The retarded impurity Green function $iG^{R}_{\mathrm{imp}}(\tau)$ obtained with (a) 6, (b) 12, (c) 18, and (d) 24 Trotter steps up to time $\tau=6/t^{*}$ using the XY method (blue diamonds). Comparison is given to the exact Green function (red dashed line). We set $U=4t^{*}$ and $V=t^{*}$.

To obtain the impurity Green function in the frequency domain, we first consider some known and general analytic properties of the retarded Green function in Eq. (3). This Green function can be written as a sum of the particle and hole contributions as

$$\begin{aligned} iG^{R}_{\mathrm{imp}}(\tau)=\theta(\tau)\sum _{j} \bigl( \bigl| \bigl\langle j \bigl| \hat{c}^{\dagger}_{1\sigma} \bigr| GS \bigr\rangle \bigr|^{2} e^{-i\omega_{j} t} + \bigl| \bigl\langle j | \hat{c}_{1\sigma} | GS \bigr\rangle \bigr|^{2} e^{i\omega_{j} t} \bigr), \end{aligned}$$

(33)

where $|j\rangle$ is an eigenstate of $\hat{H}_{\mathrm{SIAM}}$ with eigenenergy $E_{j}$, and $\omega_{j}=E_{j}-E_{GS}$. In two-site DMFT, the interacting Green function is a four-pole function [32], which limits the number of terms in the above summation to four. Moreover, in the presence of particle-hole symmetry, we have $| \langle j | \hat{c}^{\dagger}_{1\sigma} | GS \rangle |^{2}= | \langle j | \hat{c}_{1\sigma} | GS \rangle |^{2}$, and Eq. (33) can be written as

$$\begin{aligned} iG^{R}_{\mathrm{imp}}(\tau)=2 \bigl[ \alpha_{1} \cos(\omega_{1}\tau) + \alpha _{2} \cos(\omega_{2} \tau) \bigr]\theta(\tau), \end{aligned}$$

(34)

where $\alpha_{j}=\vert \langle j | \hat{c}^{\dagger}_{1\sigma} | GS \rangle \vert ^{2}$. Thus, to obtain the impurity Green function in the frequency domain as

$$\begin{aligned} G^{R}_{\mathrm{imp}}(\omega+i\eta)={}&\alpha_{1} \biggl( \frac{1}{\omega +i\eta-\omega_{1}} + \frac{1}{\omega+i\eta+\omega_{1}} \biggr) \\ &{}+\alpha_{2} \biggl( \frac{1}{\omega+i\eta-\omega_{2}} + \frac{1}{\omega+i\eta+\omega_{2}} \biggr), \end{aligned}$$

(35)

we need to extract the unknown residues $\alpha_{j}$ and poles $\omega _{j}$ by fitting an expression of the form in Eq. (34) to the measurement data of $iG^{R}_{\mathrm{imp}}(\tau)$, as shown in Figure 7(a). This method to determine $\alpha_{j}$ and $\omega_{j}$ is far more reliable and requires fewer time steps than numerically Fourier-transforming the $iG^{R}_{\mathrm{imp}}(\tau)$ data. It can also be readily generalised to larger systems by including more terms in the sum in Eq. (33). Figure 7(b) shows the real part of the impurity Green function in the frequency domain, ${\operatorname{Re}} [ G^{R}_{\mathrm{imp}}(\omega+i\eta) ]$, with residues and poles obtained from the fit in Figure 7(a), while in Figure 7(c) we plot the real part of the impurity self-energy, ${\operatorname{Re}} [ \Sigma_{\mathrm{imp}}(\omega+i\eta) ]$, obtained utilizing the Dyson equation (11). We clearly see the four-pole structure of the Green function, while the self-energy has two poles. The results are in excellent agreement with the exact solution of the two-site SIAM, with the poles of the self-energy using fitted $\alpha_{j}$ and $\omega_{j}$ differing from the exact solution by 2%.

Figure 7

The retarded impurity Green function and self-energy in the frequency domain. (a) The residues and poles of the Green function can be obtained from a fit of the form in Eq. (34) (red dashed line) to the $G^{R}_{\mathrm{imp}}(\tau)$ data from the XY method with 24 Trotter steps (blue diamonds). (b) The real part of the impurity Green function, $\operatorname{Re} [ G^{R}_{\mathrm{imp}}(\omega+i\eta) ]$ (blue line), with residues and poles obtained from the fit from (a), compared to the exact Green function (red dashed line). (c) Same as in (b), but for the self-energy $\operatorname{Re} [ \Sigma_{\mathrm{imp}}(\omega+i\eta) ]$. We set $U=4t^{*}$ and $V=t^{*}$. In (b) and (c), we have broadened the peaks with $\eta=0.01$ for clarity.

Once we have obtained the impurity Green function, and thus the impurity self-energy, we proceed according to the two-site DMFT protocol in Section 3 until self-consistency has been reached. In DMFT we are interested in the local lattice spectral function $A_{{\mathrm{latt}},jj}(\omega)$ which, at self-consistency, is given by the impurity spectral function $A_{\mathrm{imp}}(\omega)$. In the paramagnetic phase of the infinite-dimensional Hubbard model, the spectral function has a three peak structure with an upper and a lower Hubbard band, corresponding to empty and doubly occupied sites, respectively, and a quasiparticle peak with integrated spectral weight $\mathcal{Z}$ between the bands [19]. In two-site DMFT, since the self-energy has two poles, this three peak structure can be qualitatively reproduced with the spectral function [32]

$$\begin{aligned} A(\omega)=\rho_{0} \bigl[\omega+\mu- \Sigma_{\mathrm{imp}}(\omega) \bigr], \end{aligned}$$

(36)

where $\rho_{0}$ is the non-interacting density of states of the Bethe lattice. Figure 8 shows the spectral function in Eq. (36) where the impurity self-energy has been obtained both from the XY method and from exact numerics of the two-site SIAM using the interactions $U=5t^{*}$ and $U=8t^{*}$. We notice that for $U=5t^{*}$, the Hubbard bands from the XY method are slightly dislocated and the quasiparticle peak is slightly narrower compared with the exact solution of the two-site SIAM, but the agreement is still very good. The overall shape of the spectral function from the XY method is unchanged compared to the exact case. This underestimation of the width of the quasiparticle peak stems from the fact that the fitting procedure in Figure 7(a) causes the negative of the derivative of the self-energy in the XY method to be a bit larger than the exact value from the two-site SIAM, i.e.,

$$\begin{aligned} -\frac {d{\operatorname{Re}}[\Sigma^{XY}_{\mathrm{imp}}(\omega+i\eta)]}{d\omega} \bigg|_{\omega=0} \gtrsim-\frac{d{\operatorname{Re}}[\Sigma^{\mathrm{exact}}_{\mathrm{imp}} (\omega+i\eta)]}{d\omega} \bigg|_{\omega=0}, \end{aligned}$$

(37)

which leads to $\mathcal{Z}$ in Eq. (10) from the XY method to be slightly smaller than in the exact solution of the two-site SIAM, i.e., $\mathcal{Z}^{XY} \lesssim\mathcal{Z}^{\mathrm{exact}}$. For $U=8t^{*}$, the two spectral functions agree with maximum relative error of 10⁻⁸, since in this case $V=0$ is found to be the self-consistent solution, whence the Trotterized evolution operator in Eq. (27) matches full evolution operator of the two-site SIAM, and thus there is no Trotter error. We observe that in Figure 8 the central quasiparticle peak vanishes, which is characteristic of insulating behaviour. See Ref. [32] for a discussion of the artifacts of the spectral functions in two-site DMFT compared to full DMFT.

Figure 8

Spectral functions in the metallic and insulating phases. Spectral functions obtained with the XY method with 24 Trotter steps (blue line) and exact solution of the two-site SIAM (red dashed line). The parameters of the two-site SIAM are iterated to self-consistency with (a) $U=5t^{*}$ and (b) $U=8t^{*}$.

To study the transition between the two types of spectral functions in Figure 8, we plot in Figure 9 the self-consistent quasiparticle weight $\mathcal{Z}$ obtained from the XY method as a function of the interaction U for different Trotter steps N. We also show $\mathcal{Z}$ from the exact solution of the two-site SIAM for comparison. We see that the digital approach captures the correct trend of the curve, but in the metallic side underestimates to a small degree the values of $\mathcal{Z}$ for interactions close to $U=U_{c}=6t^{*}$, which is the critical interaction for Mott transition in two-site DMFT at half-filling [32]. These results are consistent with the spectral functions in Figure 8. The underestimation of $\mathcal{Z}$ can be diminished by increasing N, as shown in Figure 9. It is noteworthy to mention that two-site DMFT overestimates the quasiparticle weight compared to full DMFT for interactions $U< U_{c}$, as demonstrated in Ref. [32]. Above $U_{c}$, we find $\mathcal{Z}=0$ to be the self-consistent solution, corresponding to the insulating phase.

Figure 9

Quasiparticle weight as a function of interaction U . Self-consistent quasiparticle weight $\mathcal{Z}$ obtained from the XY method with 24 (blue diamonds), 36 (red circles), and 48 Trotter steps (yellow squares), compared to the exact solution of the two-site SIAM (purple stars). Inset: Same plot zoomed into the region around the critical interaction, $U_{c}=6t^{*}$.

6 Summary

We have proposed a quantum algorithm for two-site DMFT to be run on a small digital quantum simulator with a classical feedback loop, allowing the qualitative description of the infinite-dimensional Hubbard model in the thermodynamic limit. We have considered two alternative quantum gate decompositions consistent with state-of-the-art technology in superconducting circuits for the time-evolution operator. We found that an increasing number of Trotter steps improves the fidelity of our digital scheme to qualitatively describe the Mott transition. Our work therefore provides an interesting application for small-scale quantum devices. It also paves the way for more accurate quantum simulations of strongly correlated fermions in various lattice geometries, which are relevant to novel quantum materials, when the general self-consistency condition and larger number of qubits are used.

Appendix: single-qubit interferometry for the impurity Green function

Here, we present a measurement scheme for the retarded impurity Green function.

1.1 A.1 Definitions

The retarded zero temperature impurity Green function in the time domain can be written as

$$\begin{aligned} G^{R}_{\mathrm{imp}}(\tau)=\theta(\tau) \bigl[ G^{>}_{\mathrm{imp}}( \tau )-G^{< }_{\mathrm{imp}}(\tau) \bigr], \end{aligned}$$

(38)

where the ‘greater’ and ‘lesser’ Green functions are given by

$$\begin{aligned} &G^{>}_{\mathrm{imp}}(\tau)=-i\bigl\langle \hat{c}_{1\sigma}(\tau) \hat {c}_{1\sigma}^{\dagger}(0) \bigr\rangle , \end{aligned}$$

(39)

$$\begin{aligned} &G^{< }_{\mathrm{imp}}(\tau)=i\bigl\langle \hat{c}^{\dagger}_{1\sigma}(0) \hat {c}_{1\sigma}(\tau) \bigr\rangle , \end{aligned}$$

(40)

respectively. The average is computed in the ground-state $|GS\rangle$ of the two-site SIAM in Eq. (6). Here, σ can be either ↓ or ↑ since we are considering a spin-symmetric case (i.e., $G^{R}_{\downarrow }=G^{R}_{\uparrow}$), and the ĉ-operators are given in the Heisenberg picture with respect to $\hat{H}_{\mathrm{SIAM}}$, i.e.,

$$\begin{aligned} \hat{c}_{1\sigma}(\tau)=\hat{U}^{\dagger}(\tau)\hat{c}_{1\sigma } \hat{U}(\tau)=e^{i\tau\hat{H}_{\mathrm{SIAM}}} \hat{c}_{1\sigma} e^{-i\tau\hat{H}_{\mathrm{SIAM}}}. \end{aligned}$$

(41)

One possibility to measure the impurity Green function $G^{R}_{\mathrm{imp}}(\tau)$ is to use a single-qubit Ramsey interferometer [44] which was used in Ref. [37] in the more general non-equilibrium case. To this end, we introduce an ancilla qubit in addition to the ‘system’ qubits, raising the total number of qubits needed to implement the two-site DMFT scheme to five.

1.2 A.2 Jordan-Wigner transformation

The greater and lesser components, $G^{>}_{\mathrm{imp}}(\tau)$ and $G^{<}_{\mathrm{imp}}(\tau)$, must be written in terms of spin operators by again mapping the $\hat{c}_{1\sigma}$ and $\hat{c}^{\dagger }_{1\sigma}$ operators onto Pauli operators via the Jordan-Wigner transformation. For concreteness, we focus on the case $\sigma= \downarrow$. We obtain

$$\begin{aligned} G^{>}_{\mathrm{imp}}(\tau)={}&{-}\frac{i}{4} \bigl( \bigl\langle \hat{U}^{\dagger }(\tau)\hat{\sigma}^{x}_{1} \hat{U}(\tau) \hat{\sigma}^{x}_{1} \bigr\rangle + i\bigl\langle \hat{U}^{\dagger}(\tau)\hat{\sigma}^{x}_{1} \hat{U}(\tau) \hat{\sigma}^{y}_{1} \bigr\rangle -i\bigl\langle \hat{U}^{\dagger}(\tau)\hat {\sigma}^{y}_{1} \hat{U}( \tau) \hat{\sigma}^{x}_{1} \bigr\rangle \\ &{}+ \bigl\langle \hat{U}^{\dagger}(\tau)\hat{\sigma}^{y}_{1} \hat{U}(\tau ) \hat{\sigma}^{y}_{1} \bigr\rangle \bigr), \end{aligned}$$

(42)

and

$$\begin{aligned} G^{< }_{\mathrm{imp}}(\tau)={}&\frac{i}{4} \bigl( \bigl\langle \hat{\sigma}^{x}_{1} \hat{U}^{\dagger}(\tau) \hat{ \sigma}^{x}_{1}\hat{U}(\tau)\bigr\rangle - i\bigl\langle \hat{ \sigma}^{x}_{1} \hat{U}^{\dagger}(\tau) \hat{\sigma }^{y}_{1}\hat{U}(\tau) \bigr\rangle +i\bigl\langle \hat{ \sigma}^{y}_{1} \hat {U}^{\dagger}(\tau) \hat{ \sigma}^{x}_{1} \hat{U}(\tau)\bigr\rangle \\ &{} + \bigl\langle \hat{\sigma}^{y}_{1} \hat{U}^{\dagger}( \tau) \hat{\sigma}^{y}_{1}\hat{U}(\tau) \bigr\rangle \bigr). \end{aligned}$$

(43)

1.3 A.3 Measurement protocol

Each of the terms of the form $\langle \hat{U}^{\dagger}(\tau) \hat {\sigma}^{\alpha}_{1} \hat{U}(\tau) \hat{\sigma}^{\beta}_{1} \rangle$, where $\alpha, \beta\in\{x, y\}$, can be measured in the interferometer. This can be seen as follows. We denote the state of the system qubits by $\hat{\rho}_{\mathrm{sys}}=|GS\rangle\langle GS|$, where $|GS\rangle$ is the ground-state of the system. We initialize the ancilla qubit in the state $|0\rangle$, yielding the total density operator $\hat{\rho}_{\mathrm{tot}}=|0\rangle\langle0 | \otimes\hat {\rho}_{\mathrm{sys}}$. The total system then undergoes the following evolution:

• 1. At time $t=0$, a Hadamard gate $\hat{\sigma}_{H}=\frac {1}{\sqrt{2}} (\hat{\sigma}^{z}+\hat{\sigma}^{x} )$ is applied on the ancilla qubit, creating the superposition $|0\rangle _{\mathrm{ancilla}} \rightarrow (|0\rangle_{\mathrm{ancilla}} + |1\rangle _{\mathrm{ancilla}} )/\sqrt{2}$.

• 2. A Controlled-Pauli gate $\hat{\sigma}^{\alpha}_{1}$ is applied on the impurity qubit 1 if the ancilla qubit has state $|0\rangle$.

• 3. The system qubits undergo time evolution according to the unitary $\hat{U}(\tau)$ which is decomposed into quantum gates.

• 4. Another controlled Pauli gate $\hat{\sigma}^{\beta}_{1}$ is applied on the impurity qubit 1 if the ancilla qubit has state $|1\rangle$.

• 5. Another Hadamard gate is applied on the ancilla qubit.

Denoting the total unitary in steps 2-4 by T̂, the state of the ancilla qubit after this evolution is given by

$$\begin{aligned} \hat{\rho}_{\mathrm{ancilla}}={}&\operatorname {Tr}_{\mathrm{sys}} \bigl[ \hat{ \sigma}_{H} \hat {T} \hat{\sigma}_{H} \hat{ \rho}_{\mathrm{tot}} \hat{\sigma}_{H} \hat {T}^{\dagger} \hat{ \sigma}_{H} \bigr] \\ ={}&\frac{1+\operatorname{Re}[F(\tau)]}{2}|0\rangle \langle0 | - i\frac{ \operatorname{Im}[F(\tau)]}{2}|0 \rangle \langle1 | + i\frac{ \operatorname{Im}[F(\tau)]}{2}|1 \rangle \langle0 | \\ &{} + \frac {1-\operatorname{Re}[F(\tau)]}{2}|1 \rangle \langle1 |, \end{aligned}$$

(44)

where $F(\tau)=\operatorname {Tr}_{\mathrm{sys}} [ \hat{T}_{1}^{\dagger }(\tau) \hat{T}_{0}(\tau) \hat{\rho}_{\mathrm{sys}} ]$. We have denoted the controlled unitaries as $\hat{T}_{1}(\tau)=\hat {\sigma}^{\alpha}_{1} \hat{U}(\tau)$ and $\hat{T}_{0}(\tau)= \hat {U}(\tau)\hat{\sigma}^{\beta}_{1}$. Note that since the same $\hat {U}(\tau)$ appears in both unitaries, only the Pauli gates $\hat {\sigma}^{\alpha/ \beta}_{1}$ need to be controlled, as described above. Note also that $F(\tau)=\langle\hat{U}^{\dagger}(\tau) \hat {\sigma}^{\alpha}_{1} \hat{U}(\tau) \hat{\sigma}^{\beta}_{1} \rangle$. We can rewrite the state of the ancilla qubit as

$$\begin{aligned} \hat{\rho}_{\mathrm{ancilla}}=\frac{1}{2} \bigl(\hat{I} + \operatorname{Re} \bigl[F(\tau)\bigr]\hat{\sigma}_{z} +\operatorname{Im}\bigl[F(\tau)\bigr] \hat{\sigma}_{y} \bigr), \end{aligned}$$

(45)

whence $\operatorname {Tr}_{\mathrm{ancilla}} [\hat{\rho}_{\mathrm{ancilla}} \hat {\sigma}^{z} ]=\operatorname{Re}[F(\tau)]$, and $\operatorname {Tr}_{\mathrm{ancilla}} [\hat{\rho}_{\mathrm{ancilla}} \hat {\sigma}^{y} ]=\operatorname{Im}[F(\tau)]$. Thus, repeated measurements of the $\hat{\sigma}^{z}$ and $\hat{\sigma }^{y}$ components of the ancilla qubit yield the real and imaginary parts of the term $\langle \hat{U}^{\dagger}(\tau) \hat{\sigma}^{\alpha }_{1} \hat{U}(\tau) \hat{\sigma}^{\beta}_{1} \rangle$. See Figure 10 for the quantum network of the scheme.

Figure 10

Quantum network to measure the $\pmb{\langle GS | \hat{U}^{\dagger}(\tau) \hat{\sigma}^{\alpha}_{1} \hat{U}(\tau) \hat{\sigma}^{\beta}_{1} |GS\rangle}$ contribution to the Green function $\pmb{G^{R}_{\mathrm{imp}}(\tau)}$ . The time-evolution operator $\hat{U}(\tau)$ is composed of a set of quantum gates according to the main text.