Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution

TL;DR

This paper introduces resource-efficient quantum algorithms for finding eigenstates and thermal states, avoiding deep circuits and complex optimization, with demonstrated implementations on quantum simulators and hardware.

Contribution

It presents quantum imaginary time evolution and quantum Lanczos algorithms as practical, low-resource methods for quantum state preparation, improving over existing techniques.

Findings

Algorithms require less space and time per iteration than classical methods.

Implementation demonstrated on Rigetti quantum virtual machine and Aspen-1 hardware.

Potential for efficient quantum simulation of thermal states and excited states.

Abstract

The accurate computation of Hamiltonian ground, excited, and thermal states on quantum computers stands to impact many problems in the physical and computer sciences, from quantum simulation to machine learning. Given the challenges posed in constructing large-scale quantum computers, these tasks should be carried out in a resource-efficient way. In this regard, existing techniques based on phase estimation or variational algorithms display potential disadvantages; phase estimation requires deep circuits with ancillae, that are hard to execute reliably without error correction, while variational algorithms, while flexible with respect to circuit depth, entail additional high-dimensional classical optimization. Here, we introduce the quantum imaginary time evolution and quantum Lanczos algorithms, which are analogues of classical algorithms for finding ground and excited states. Compared to their classical counterparts, they require exponentially less space and time per iteration, and can be implemented without deep circuits and ancillae, or high-dimensional optimization. We furthermore discuss quantum imaginary time evolution as a subroutine to generate Gibbs averages through an analog of minimally entangled typical thermal states. Finally, we demonstrate the potential of these algorithms via an implementation using exact classical emulation as well as through prototype circuits on the Rigetti quantum virtual machine and Aspen-1 quantum processing unit.

Tables

Table 1. Table 1: QPUs: 1-qubit QITE and QLanczos.

Trotter stepsize	nSamples	$\delta$	s	$ϵ$
0.2	100000	0.01	0.75	${10}^{-2}$

Table 2. Table 2: QPUs: 2-qubit QITE and QLanczos.

Trotter stepsize	nSamples	$\delta$	s	$ϵ$
0.5	100000	0.1	0.75	${10}^{-2}$

Table 3. Table 3: QPUs: 1-qubit METTS.

$\beta$	Trotter stepsize	nSamples	nMETTs	$\delta$
1.5	0.15	1500	70	0.01
2.0	0.20	1500	70	0.01
3.0	0.30	1500	70	0.01
4.0	0.40	1500	70	0.01

Table 4. Table 4: QVM: 2-qubit QITE and QLanczos.

Trotter stepsize	nSamples	$\delta$	s	$ϵ$
0.5	100000	0.1	0.75	${10}^{-2}$

Table 5. Table 5: QVM: 1-qubit METTS.

$\beta$	Trotter stepsize	nSamples	nMETTs	$\delta$
1.0	0.10	1500	70	0.01
1.5	0.15	1500	70	0.01
2.0	0.20	1500	70	0.01
3.0	0.30	1500	70	0.01
4.0	0.40	1500	70	0.01

Table 6. Table 6: QVM: 2-qubit METTS.

$\beta$	Trotter stepsize	nSamples	nMETTs	$\delta$
1.0	0.10	30000	100	0.1
1.5	0.15	30000	100	0.1
2.0	0.20	30000	100	0.1
3.0	0.30	30000	100	0.1
4.0	0.40	30000	100	0.1

Table 7. Table 7: VQE simulation parameters for (a) 1D Heisenberg with applied field and (b) 1D AFM transverse field Ising. Conv. refers to the convergence criterion used. We note for the last case, the VQE optimization could not reach within 1 % percent \% of the ground state energy, so we set the convergence criterion to 2 % percent \%

model	n-site	conv.	$\alpha$	$\gamma$	$d$	$N$	$P_{total}$
a	4	1$\%$	0.602	0.101	8	800	25,600
a	6	1$\%$	0.602	0.101	20	8400	403,200
b	4	1$\%$	0.602	0.101	12	800	12,800
b	6	2$\%$	0.602	0.101	12	2890	69,360

Table 8. Table 8: Hyperparameters sweep for 6-site 1D Heisenberg model using a circuit depth of 20 for a total of 9000 optimization steps. The step at which the calculation converged is recorded under column T 𝑇 T . ’-’ indicates that VQE failed to converge.

$\alpha$	$\gamma$	$T$
0.400	0.066	-
0.400	0.101	-
0.400	0.133	-
0.602	0.066	-
0.602	0.101	8400
0.602	0.133	8800
0.800	0.066	-
0.800	0.101	-
0.800	0.133	-

Table 9. Table 9: QITE simulation parameters for (a) 1D Heisenberg with applied field and (b) 1D AFM transverse field Ising. Conv. indicates the convergence criterion used. We used 2 % percent \% for the final calculation to facilitate comparison with VQE which failed to converge to within 1 % percent \% .

model	n-site	conv.	$\Delta \tau$	$D$	$T$	K	$P_{total}$
a	4	1$\%$	0.1	4	7	4	12,544
a	6	1$\%$	0.1	4	17	6	47,872
b	4	1$\%$	0.2	4	7	4	12,544
b	6	2$\%$	0.2	4	8	6	22,528

Table 10. Table 10: Total Pauli string expectation values in VQE and QITE for (a) 1D Heisenberg with applied field and (b) 1D AFM transverse field Ising. The total number of Pauli strings to be measured in QITE can be further reduced by using only Pauli strings with only an odd number of Y ^ ^ 𝑌 \hat{Y} operators due to the real nature of the Hamiltonian. We show this reduced number in brackets.

model	n-site	VQE	QITE
a	4	25,600	12,544(5,880)
a	6	403,200	47,872(22,440)
b	4	12,800	12,544(5,880)
b	6	69,360	22,528(10,560)

Equations

$$e^{-\beta \hat{H}}=(e^{-\Delta \tau \hat{h}[1]}{e}^{-\Delta \tau \hat{h}[2]}\dots {)}^n+\mathcal{O}\left(\Delta \tau \right);n=\frac{\beta }{\Delta \tau }$$

$$|{\Psi }^{\prime }⟩=e^{-\Delta \tau \hat{h}[m]}|\Psi ⟩.$$

$$c=1-2\Delta \tau ⟨\Psi |\hat{h}[m]|\Psi ⟩+\mathcal{O}({\Delta \tau }^2).$$

$$\hat{A}[m]=\sum _{i_1\dots i_k}a[m{]}_{i_1\dots i_k}{\hat{\sigma }}_{i_1}\dots {\hat{\sigma }}_{i_k}\equiv \sum _{I}a[m{]}_I{\hat{\sigma }}_I.$$

$$S_{I,I^{\prime }}=⟨\Psi |{\hat{\sigma }}_I^{†}{\hat{\sigma }}_{I^{\prime }}^{}|\Psi ⟩,b_I=\frac{-i}{\sqrt{c}}⟨\Psi |{\hat{\sigma }}_I^{†}\hat{h}[m]|\Psi ⟩.$$

$$|{\bar{\Psi }}^{\prime }⟩\equiv c^{-1/2}{e}^{-\Delta \tau \hat{h}[l]}|\Psi ⟩=e^{-i\Delta \tau \hat{A}[l]}|\Psi ⟩.$$

$$\hat{A}[l]=\sum _{i_1{i}_2\dots i_D}a[l{]}_{i_1{i}_2\dots i_D}{\hat{\sigma }}_{i_1}{\hat{\sigma }}_{i_2}\dots {\hat{\sigma }}_{i_D}=\sum _{I}a[l{]}_I{\hat{\sigma }}_I,$$

$$|{\Delta }_0⟩=\frac{|{\bar{\Psi }}^{\prime }⟩-|\Psi ⟩}{\Delta \tau },|\Delta ⟩=-i\hat{A}[l]|\Psi ⟩.$$

$$f(a[l])=f_0+\sum _I{b}_{I}a[l{]}_I+\sum _{IJ}a[l{]}_I{S}_{IJ}a[l{]}_J$$

$$f_0$$

$$S_{IJ}$$

$$b_I$$

$$\left(\mathbf{S}+{\mathbf{S}}^T\right)\mathbf{a}[l]=-\mathbf{b}$$

$$\begin{array}{rl}{b}_I& =-2\text{mbox}Im\left[⟨\Psi |{\hat{\sigma }}_I^{†}|{\Delta }_0⟩\right],\end{array}$$

$$y(D=1)=1$$

$$y(D=2)=6$$

$$\hat{I}\hat{Y},\hat{X}\hat{Y},\hat{Z}\hat{Y},\hat{Y}\hat{I},\hat{Y}\hat{X},\hat{Y}\hat{Z}.$$

$$y(D=3)=28$$

$$\begin{array}[]{lll}\hat{I}\hat{I}\hat{Y},&\dots,&\hat{I}\hat{Y}\hat{Z}\\ \hat{X}\hat{I}\hat{Y},&\dots,&\hat{X}\hat{Y}\hat{Z}\\ \hat{Y}\hat{I}\hat{I},&\dots,&\hat{Y}\hat{Z}\hat{Z}\\ \hat{Z}\hat{I}\hat{Y},&\dots,&\hat{Z}\hat{Y}\hat{Z}\\ \end{array}$$

$$y(D+1)=3y(D)+(4^D-y(D)).$$

$$y(D)=2^D\frac{2^D-1}{2}.$$

$$H=\sum _{l=1}^m\hat{h}[l],$$

$$\frac{e^{-\beta \hat{H}}|{\Psi }_0⟩}{\Vert e^{-\beta \hat{H}}|{\Psi }_0⟩\Vert }\simeq \frac{{\left(e^{-\Delta \tau \hat{h}[1]}\dots e^{-\Delta \tau \hat{h}[m]}\right)}^n|{\Psi }_0⟩}{\Vert {\left(e^{-\Delta \tau \hat{h}[1]}\dots e^{-\Delta \tau \hat{h}[m]}\right)}^n|{\Psi }_0⟩\Vert }.$$

$$\begin{array}{rl}{C}_s(\hat{A},\hat{B})& =⟨{\Psi }_s|\hat{A}\otimes \hat{B}|{\Psi }_s⟩-⟨{\Psi }_s|\hat{A}|{\Psi }_s⟩⟨{\Psi }_s|\hat{B}|{\Psi }_s⟩\\ & \le \Vert \hat{A}\Vert \Vert \hat{B}\Vert e^{-\text{dist}(A,B)/C}.\end{array}$$

$$N_q=k(2C{)}^d{\ln }^d\left(2\sqrt{2}nm{\epsilon }^{-1}\right)$$

$$\Vert |{\Psi }_{mn}⟩-|{\Phi }_{mn}⟩\Vert \le \epsilon .$$

$$\Vert |{\Psi }_s⟩-|{\Phi }_s⟩\Vert =\Vert |{\Psi }_s⟩-{\hat{U}}_s|{\Phi }_{s-1}⟩\Vert$$

$$\le$$

$$|{\Psi }_s⟩=\frac{e^{-\Delta \tau \hat{h}[i_s]}|{\Psi }_{s-1}⟩}{\Vert e^{-\Delta \tau \hat{h}[i_s]}|{\Psi }_{s-1}⟩\Vert },$$

$$\begin{array}{rl}& {\Vert {\text{tr}}_{\backslash R_v}(|{\Psi }_s⟩⟨{\Psi }_s|)-{\text{tr}}_{\backslash R_v}(|{\Psi }_{s-1}⟩⟨{\Psi }_{s-1}|)\Vert }_1\\ & \le \Vert e^{\Delta \tau \hat{h}[i_s]}{\Vert }^{-1}{e}^{-\frac{v}{C}}\le 2e^{-\frac{v}{C}},\end{array}$$

$$\Vert e^{-\Delta \tau \hat{h}[i_s]}\Vert \ge \Vert I-\Delta \tau \hat{h}[i_s]\Vert \ge 1-\Delta \tau \ge 1/2.$$

$$\Vert |\eta {⟩}_{AB}-(I\otimes \hat{V})|\nu {⟩}_{AB}\Vert \le 2\sqrt{\delta }.$$

Full text

An important application for a quantum computer is to compute the ground-state $\Psi$ of a Hamiltonian $\hat{H}$ Feynman (1982); Abrams and Lloyd (1997). This arises in simulations, for example, of the electronic structure of molecules and materials, Lloyd (1996); Aspuru-Guzik et al. (2005); Kandala et al. (2017a, 2019) as well as in more general optimization problems. While efficient ground-state determination cannot be guaranteed for all Hamiltonians, as this is a QMA-hard problem Kempe et al. (2006), several heuristic quantum algorithms have been proposed, including adiabatic state preparation with quantum phase estimation Farhi et al. (2000); Kitaev (1995) (QPE) and quantum-classical variational algorithms, such as the quantum approximate optimization algorithm Farhi et al. (2014); Otterbach et al. (2017); Moll et al. (2018) and variational quantum eigensolver Peruzzo et al. (2014); McClean et al. (2016); Grimsley et al. (2019). Despite many advances, these algorithms also have potential disadvantages, especially in the context of near-term quantum computing architectures with limited quantum resources. For example, phase estimation produces a nearly exact eigenstate, but appears impractical without error correction, while variational algorithms, though somewhat robust to coherent errors, are limited in accuracy by a fixed Ansatz, and involve high-dimensional noisy classical optimizations McClean et al. (2018).

In classical simulations, different strategies are employed to numerically determine nearly exact ground-states. One popular approach is imaginary-time evolution, which expresses the ground-state as the long-time limit of the imaginary-time Schrödinger equation $-\partial_{\beta}|\Phi(\beta)\rangle=\hat{H}|\Phi(\beta)\rangle$, $|\Psi\rangle=\lim_{\beta\to\infty}\frac{|\Phi(\beta)\rangle}{\|\Phi(\beta)\|}$ (for $\langle\Phi(0)|\Psi\rangle\neq 0$). Unlike variational algorithms with a fixed Ansatz, imaginary-time evolution always converges to the ground-state, as distinguished from imaginary-time Ansatz optimization McArdle et al. (2018). Another family of approaches are variants of the iterative Lanczos method Lanczos (1950). The Lanczos iteration constructs the Hamiltonian matrix $\mathbf{H}$ in a Krylov subspace $\{|\Phi\rangle,\hat{H}|\Phi\rangle,\hat{H}^{2}|\Phi\rangle\ldots\}$; diagonalizing $\mathbf{H}$ yields a variational estimate of the ground-state which tends to $|\Psi\rangle$ for a large number of iterations. For an $N$-qubit Hamiltonian, the classical complexity of imaginary time evolution and Lanczos algorithm scales as $\sim\exp{(\mathcal{O}(N))}$ in space and time. Exponential space comes from storing $\Phi(\beta)$ or the Lanczos vector, while exponential time comes from the cost of Hamiltonian multiplication $\hat{H}|\Phi\rangle$, as well as, in principle, though not in practice, the $N$-dependence of the number of propagation steps or Lanczos iterations. Thus it is natural to consider quantum versions of these algorithms that can overcome the exponential bottlenecks.

Here we describe the quantum imaginary time evolution (QITE), the quantum Lanczos (QLanczos) and the quantum analog of the minimally entangled typical thermal states (QMETTS) algorithm, to determine ground-states, ground and excited states and thermal states on a quantum computer. Under the assumption of finite correlation length, these methods rigorously use exponentially reduced space and time per propagation step or iteration, compared to their direct classical counterparts. Even when such assumptions do not hold, the inexact versions of the QITE and QLanczos algorithms remain valid heuristics that can be applied within a limited computational budget, and offer advantages over existing ground-state quantum algorithms, as they do not use deep circuits and converge to their solutions without non-linear optimization. A crucial common component is the efficient implementation of the non-Hermitian operation of an imaginary-time step $e^{-{\Delta\tau}\hat{H}}$ (for small ${\Delta\tau}$) assuming a finite correlation length in the state. Non-Hermitian operations are not natural on a quantum computer and are usually achieved using ancillae and postselection, but we describe how to implement imaginary time evolution on a given state without these resources. The lack of ancillae and complex circuits make our algorithms potentially suitable for near-term quantum architectures. We demonstrate the algorithms on spin and fermionic Hamiltonians using exact classical emulation, and demonstrate proof-of-concept implementations on the Rigetti quantum virtual machine (QVM) and Aspen-1 quantum processing units (QPUs).

Quantum Imaginary Time Evolution. Define a geometric $k$-local Hamiltonian $\hat{H}=\sum_{m}\hat{h}[m]$ (where each term $\hat{h}[m]$ acts on at most $k$ neighboring qubits on an underlying graph) and a Trotter decomposition of the corresponding imaginary-time evolution,

[TABLE]

applied to a state $|\Psi\rangle$. After a single Trotter step, we have

[TABLE]

The basic idea is that the normalized state $|\bar{\Psi}^{\prime}\rangle=|\Psi^{\prime}\rangle/\|\Psi^{\prime}\|$ is generated from $|\Psi\rangle$ by a unitary operator $e^{-i{\Delta\tau}\hat{A}[m]}$ acting on a neighborhood of the qubits acted on by $\hat{h}[m]$, where $\hat{A}[m]$ can be determined from tomography of $|\Psi\rangle$ in this neighborhood up to controllable errors. This is illustrated by the simple example where $|\Psi\rangle$ is a product state. The squared norm $c=\langle\Psi|e^{-2{\Delta\tau}\hat{h}[m]}|\Psi\rangle$ can be calculated from the expectation value of $\hat{h}[m]$, requiring measurements over $k$ qubits,

[TABLE]

Because $|\Psi\rangle$ is a product state, $|\Psi^{\prime}\rangle$ is obtained applying the unitary operator $e^{-i{\Delta\tau}\hat{A}[m]}$ also on $k$ qubits. $\hat{A}[m]$ can be expanded in terms of an operator basis, e.g. the Pauli basis $\{\hat{\sigma}_{i}\}$ on $k$ qubits,

[TABLE]

Up to $\mathcal{O}({\Delta\tau})$, the coefficients $a[m]_{I}$ are defined by the linear system $\mathbf{S}\mathbf{a}[m]=\mathbf{b}$ where the elements of $\mathbf{S}$ and $\mathbf{b}$ are expectation values over $k$ qubits,

[TABLE]

In general, $\mathbf{S}$ has a null space; to ensure $\mathbf{a}[m]$ is real, we minimize $\|\bar{\Psi}^{\prime}-(1-i{\Delta\tau}\hat{A}[m])\Psi\|^{2}$ w.r.t. real variations in $\mathbf{a}[m]$ (see SI). Because the solution is determined from a linear problem, there are no local minima.

In this simple case, the normalized result of the imaginary time evolution step could be represented by a unitary update over $k$ qubits, because $|\Psi\rangle$ had correlation length zero. After the initial step, this is no longer the case. However, for a more general $|\Psi\rangle$ with finite correlations over at most $C$ qubits (i.e. correlations between observables separated by distance $L$ are bounded by $\exp(-L/C)$), $|\bar{\Psi}^{\prime}\rangle$ can be generated by a unitary acting on a domain of width at most $O(C)$ qubits surrounding the qubits acted on by $\hat{h}[m]$. This follows from Uhlmann’s theorem Uhlmann (1976), which states that two pure states with marginals close to each other must be related by a unitary transformation on the purifying subsystem (see SI). The unitary $e^{-i{\Delta\tau}\hat{A}[m]}$ can then be determined by measurements and solving the least squares problem in this domain (Fig. 1). For example, for a nearest-neighbor local Hamiltonian on a $d$-dimension cubic lattice, the domain size $D$ is bounded by $O(C^{d})$. In many physical systems, we expect the maximum correlation length throughout the Trotter steps to increase with $\beta$ and saturate for $C_{\text{max}}\ll N$ Hastings and Koma (2006). Fig. 1 shows the mutual information between qubits $i$ and $j$ as a function of imaginary time in the 1D and 2D ferromagnetic transverse field Ising models computed by tensor network simulation (see SI), demonstrating a monotonic increase and clear saturation.

The above replacement of imaginary time evolution steps by unitary updates can be extended to more general Hamiltonians, such as ones with long-range interactions and fermionic Hamiltonians. For fermions, in particular, the locality of the corresponding qubit Hamiltonian depends on the qubit mapping. In principle, a geometric $k$-local fermionic Hamiltonian can be mapped to a geometric local qubit Hamiltonian Bravyi and Kitaev (2002); Verstraete and Cirac (2005), allowing above techniques to be applied directly. Alternatively, we conjecture that by constructing Eq. (4) with a local fermionic basis, the unitary update can be constructed over a domain size $D\sim O(C^{d})$, $C$ being the fermionic correlation length (see SI).

Cost of QITE. The number of measurements and classical storage at a given time step (starting propagation from a product state) is bounded by $\exp(O(C^{d}))$ (with $C$ the correlation length at that time step), since each unitary at that step acts on at most $O(C^{d})$ sites; classical solution of the least squares problem has a similar scaling $\exp(O(C^{d}))$, as does the synthesis and application as a quantum circuit (composed of two-qubit gates) of the unitary $e^{-i{\Delta\tau}\hat{A}[m]}$. Thus, space and time requirements are bounded by exponentials in $C^{d}$, but are polynomial in $N$ when one is interested in a local approximation of the state (or quasi-polynomial for a global approximation); the polynomial in $N$ comes from the number of terms in $H$ (see SI for details).

The exponential dependence on $C^{d}$ can be greatly reduced in many cases, for example if $\hat{A}[m]$ has a locality structure, e.g. if it is (approximately) a $p$-local Hamiltonian (i.e. all $a[m]_{i_{1}\ldots i_{k}}$ in Eq. (4) are zero except for those where at most $p$ of the $\hat{\sigma}_{i}$ operators differ from the identity) then the cost of tomography becomes only $C^{O(dp)}$, while the cost of finding and implementing the unitary is $O(pC^{d}T_{e})$, $T_{e}$ being the cost of computing one entry of $\hat{A}[m]$ Berry et al. (2015). If we assume further that $\hat{A}[m]$ is geometric local, the cost of tomography is reduced further to $O(pC^{d})$. However, it is important to note that even if $C$ is too large to construct the unitaries exactly, we can still run the algorithm as a heuristic, truncating the unitary updates to domain sizes that fit the computational budget. This gives the inexact QITE algorithm, described and studied in detail below.

Compared to a direct classical implementation of imaginary time evolution, the cost of a QITE time-step (for bounded correlation length $C$) is linear in $N$ in space and polynomial in $N$ in time, thus giving an exponential reduction in space and time. Note that a finite correlation length $C_{0}$ in the ground-state does not generally imply an efficient classical strategy. In the SI, we analyze multiple classical heuristics under the assumption of finite ground-state correlations, including: truncating the problem size at the ground-state correlation length $C_{0}$, classical simulation in the Heisenberg representation, and tensor network calculations Vidal (2004a); Schollwöck (2011); Schuch et al. (2007); Haferkamp et al. (2018).

Inexact QITE. Given limited resources, for example on near-term devices, we can choose to measure and construct the unitary over a domain $D$ smaller than induced by correlations, to fit the computational budget. For example, if $D=1$, this gives a mean-field approximation of the imaginary time evolution, and larger $D$ gives successively better approximations to the ground-state. Importantly, while the unitary is no longer an exact representation of the imaginary time evolution, there is no issue of a local minimum in its construction, although the energy is no longer guaranteed to decrease at every step. In this case, one can apply inexact imaginary time evolution until the energy stops decreasing; the energy will still be a variational upper bound. One can also use the quantum Lanczos algorithm, described later.

Illustrative QITE calculations. To illustrate the QITE algorithm, we have carried out exact classical emulations (assuming perfect expectation values and gates) for several Hamiltonians (see SI): short-range 1D Heisenberg (with and without a field); 1D AFM transverse-field Ising; long-range 1D Heisenberg with spin-spin coupling $J_{ij}=(|i-j|+1)^{-1}$; 1D Hubbard at half-filling; a 6-qubit MAXCUT Farhi et al. (2014); Otterbach et al. (2017); Moll et al. (2018) instance, and a minimal basis 2-qubit dihydrogen molecular Hamiltonian O’Malley et al. (2016). To assess the feasibility of implementation on near-term quantum devices, we have carried out noisy classical emulation (sampling expectation values and with an error model) using the Rigetti quantum virtual machine (QVM) and a physical simulation using the Rigetti Aspen-1 QPUs, for a single-qubit field model $(\hat{X}+\hat{Z})/\sqrt{2}$ Lamm and Lawrence (2018) and a 1D AFM transverse-field Ising model. We also carried out measurement resource estimates for QITE on the short-range 1D Heisenberg (with field) model studied in Ref. Kandala et al. (2017a) with VQE, and the 1D AFM transverse-field Ising model; we compared with resource estimates using the publicly available VQE implementation in IBM’s Qiskit. We carried out QITE using different fixed domain sizes $D$ for the unitary or fermionic unitary (see SI for descriptions of simulations and models).

Fig. 2a-2f and 3 show the energy obtained by QITE as a function of $\beta$ and $D$ for the various models. As we increase $D$, the asymptotic ($\beta\to\infty$) energies rapidly converge to the exact ground-state. For small $D$, the inexact QITE tracks the exact QITE for a time until the correlation length exceeds $D$. Afterwards, it may go down or up. The non-monotonic behavior is strongest for small domains; in the MAXCUT example, the smallest domain $D=2$ gives an oscillating energy; the first point at which the energy stops decreasing is a reasonable estimate of the ground-state energy. In all models, increasing $D$ past a maximum value (less than $N$) no longer affects the asymptotic energy, showing that the correlations have saturated (this is true even in the MAXCUT instance).

Figs. 2g, 2h show an estimate from classical emulation of the number of Pauli string expectation values to be measured in the QITE algorithm as well as the hardware-efficient VQE ansatz (using the optimization protocol in Ref. Kandala et al. (2017a)) to obtain an energy accuracy of 1% in the 1D Heisenberg model with field $J=B=1$ (Fig. 2g) and 1% or 2% in the 1D AFM transverse-field Ising model (Fig. 2h, the looser threshold was chosen to enable convergence of VQE). QITE is competitive with VQE for the 4-site model and requires significantly fewer measurements in the 6-site model. While the number of measurements could potentially be reduced in VQE by different optimizers and Ansätze, the data suggests that QITE is a promising alternative to VQE on near-term devices.

Figs. 2e and 2f show the results of running the QITE algorithm on Rigetti’s QVM and Aspen-1 QPUs for 1- and 2- qubits, respectively. The error bars are due to gate, readout, incoherent and cross-talk errors. Sufficient samples were used to ensure that sampling error is negligible. Encouragingly for near-term simulations, despite these errors it is possible to converge to a ground-state energy close to the exact energy for the 1-qubit case. This result reflects a robustness that is sometimes informally observed in imaginary time evolution algorithms in which the ground state energy is approached even if the imaginary time step is not perfectly implemented. In the 2-qubit case, although the QITE energy converges, there is a systematic shift which is reproduced on the QVM using available noise parameters for readout, decoherence and depolarizing noise Rig . Remaining discrepancies between the emulator and hardware are likely attributable to cross-talk between parallel gates not included in the noise model (see SI). However, reducing decoherence and depolarizing errors in the QVM or using different sets of qubits with improved noise characteristics (see SI) all lead to improved convergence to the exact ground-state energy.

Quantum Lanczos algorithm. Given the QITE subroutine, we now consider how to formulate a quantum Lanczos algorithm, which is an especially economical realization of a quantum subspace method McClean et al. (2017); Colless et al. (2018). An important practical motivation is that the Lanczos algorithm typically converges much more quickly than imaginary time evolution, and often in physical simulations only tens of iterations are needed to converge to good precision. In addition, Lanczos provides a natural way to compute excited states. Consider the sequence of imaginary time vectors $|\Phi_{l}\rangle=e^{-l{\Delta\tau}\hat{H}}|\Phi\rangle$, $l=0,1,\ldots n$, where $c_{l}=\|\Phi_{l}\|$. In QLanczos, we consider the vectors after even numbers of time steps $|\Phi_{0}\rangle,|\Phi_{2}\rangle\ldots$ to form a basis for the ground-state. (SI describes the equivalent treatment in terms of normalized imaginary time vectors). These vectors define an overlap matrix whose elements can be computed entirely from norms, $S_{ll^{\prime}}=\langle\Phi_{l}|\Phi_{l^{\prime}}\rangle=c^{2}_{(l+l^{\prime})/2}$, where $c_{(l+l^{\prime})/2}$ is the norm of another integer time step vector, and the overlap matrix elements for $n/2$ vectors can be accumulated for free after $n$ steps of time evolution. The Hamiltonian matrix elements satisfy the identity $H_{ll^{\prime}}=\langle\Phi_{l}|\hat{H}|\Phi_{l^{\prime}}\rangle=\langle\Phi_{(l+l^{\prime})/2}|\hat{H}|\Phi_{(l+l^{\prime})/2}\rangle$. Although the Hamiltonian has $\sim n^{2}$ matrix elements in the basis of the $\Phi_{l}$ states, there are only $\sim n$ unique elements, and importantly, each is a simple expectation value of the energy during the imaginary time evolution. This economy of matrix elements is a property shared with the classical Lanczos algorithm. Whereas the classical Lanczos iteration builds a Krylov space in powers of $\hat{H}$, QLanczos builds a Krylov space in powers of $e^{-2{\Delta\tau}\hat{H}}$; in the limit of small ${\Delta\tau}$ these Krylov spaces are identical. Diagonalization of the QLanczos Hamiltonian matrix is guaranteed to give a ground-state energy lower than that of the last imaginary time vector $\Phi_{n}$ (while higher roots approximate excited states).

With a limited computational budget, we can use inexact QITE to generate $\Phi_{l}$, $\Phi_{l}^{\prime}$. However, in this case the above expressions for $S_{ll^{\prime}}$ and $H_{ll^{\prime}}$ in terms of expectation values no longer exactly hold, which can create numerical issues (e.g. the overlap may no longer be positive). To handle this, as well as errors due to noise and sampling in real experiments, the QLanczos algorithm needs to be stabilized by ensuring that successive vectors are not nearly linearly dependent (see SI).

We demonstrate the QLanczos algorithm using classical emulation on the 1D Heisenberg Hamiltonian, as used for the QITE algorithm in Fig. 2 (see SI). Using exact QITE (large domains) to generate matrix elements, exact quantum Lanczos converges much more rapidly than imaginary time evolution. Convergence of inexact QITE (small domains), however, can both be faster and reach lower energies than inexact quantum Lanczos. We also assess the feasibility of QLanczos in presence of noise, using emulated noise on the Rigetti QVM as well as on the Rigetti Aspen-1 QPUs. In Fig. 2, we see that QLanczos also provides more rapid convergence than QITE with both noisy classical emulation as well as on the physical device for 1 and 2 qubits.

Quantum thermal averages. The QITE subroutine can be used in a range of other algorithms. For example, we discuss how to compute thermal averages \mathrm{Tr}\big{[}\hat{O}e^{-\beta\hat{H}}\big{]}/\mathrm{Tr}\big{[}e^{-\beta\hat{H}}\big{]} using imaginary time evolution. Several procedures have been proposed for quantum thermal averaging, ranging from generating the finite-temperature state explicitly by equilibration with a bath Terhal and DiVincenzo (2000), to a quantum analog of Metropolis sampling Temme et al. (2011a) that relies on phase estimation, as well as methods based on ancilla based Hamiltonian simulation with post-selection Chowdhury and Somma (2017) and approaches based on recovery maps Brandão and Kastoryano (2019). However, given a method for imaginary time evolution, one can generate thermal averages of observables without any ancillae or deep circuits. This can be done by adapting to the quantum setting the classical minimally entangled typical thermal state (METTS) algorithm White (2009a); Stoudenmire and White (2010a), which generates a Markov chain from which the thermal average can be sampled. The QMETTS algorithm can be carried out as follows (i) start from a product state, carry out imaginary-time evolution (using QITE) up to time $\beta$ (ii) measure the expectation value of $\hat{O}$ to produce its thermal average (iii) measure a product operator such as $\hat{Z}_{1}\hat{Z}_{2}\ldots\hat{Z}_{N}$, to collapse back onto a random product state (iv) repeat (i). Note that in step (iii) one can measure in any product basis, and randomizing the product basis can be used to reduce the autocorrelation time and avoid ergodicity problems in sampling. In Fig. 4 we show the results of quantum METTS (using exact classical emulation) for the thermal average $\langle\hat{H}\rangle$ as a function of temperature $\beta$, for the 6-site Heisenberg model for several temperatures and domain sizes; sufficiently large $D$ converges to the exact thermal average at each $\beta$; error bars reflect only finite QMETTS samples. We also show an implementation of quantum METTS on the Aspen-1 QPU and QVM with a 1-qubit field model (Fig. 4b), and using the QVM for a 2-qubit AFM transverse field Ising model (Fig. 4c).

Conclusions In summary, the quantum analogs of imaginary-time evolution, Lanczos and METTS algorithms we have presented enable a new class of eigenstate and thermal state quantum simulations, that can be carried out without ancillae or deep circuits and that, for bounded correlation length, achieve exponential reductions in space and time per iteration relative to known classical counterparts. Encouragingly, these algorithms appear useful in conjunction with near-term quantum architectures, and serve to demonstrate the power of quantum elevations of classical simulation techniques, in the continuing search for quantum supremacy.

Acknowledgments. MM, GKC, FGSLB, ATKT, AJM were supported by the US NSF via RAISE-TAQS CCF 1839204. MJO’R was supported by an NSF graduate fellowship via grant No. DEG-1745301; the tensor network algorithms were developed with the support of the US DOD via MURI FA9550-18-1-0095. EY was supported by a Google fellowship. CS was supported by the US DOE via DE-SC0019374. GKC is a Simons Investigator in Physics and a member of the Simons Collaboration on the Many-Electron Problem. The Rigetti computations were made possible by a generous grant through Rigetti Quantum Cloud services supported by the CQIA-Rigetti Partnership Program. We thank GH Low, JR McClean, R Babbush for discussions, and the Rigetti team for help with the QVM and QPU simulations.

Author contributions and data availability.

MM, CS, GKC designed the algorithms. FGSLB established the mathematical proofs and error estimates. EY and MJO’R performed classical tensor network simulations. MM, CS, ATKT carried out classical exact emulations. ATKT and AJM designed and carried out the Rigetti QVM and QPU experiments. All authors contributed to the discussion of results and writing of the manuscript. The code used to generate the data presented in this study can be publicly accessed on GitHub at $\mathsf{https://github.com/mariomotta/QITE.git}$

Appendix A SUPPLEMENTAL INFORMATION

A.1 Representing imaginary-time evolution by unitary maps

In this section, we discuss how to emulate imaginary time evolution by measurement-assisted unitary circuits acting on suitable domains.

As discussed in the main text, we map the scaled non-unitary action of $e^{-\Delta\tau\hat{h}[l]}$ on a state $\Psi$ to that of a unitary $e^{-i\Delta\tau\hat{A}[l]}$, i.e.

[TABLE]

where $c=\langle\Psi|e^{-2\Delta\tau\hat{h}[l]}|\Psi\rangle$. $\hat{h}[l]$ acts on $k$ qubits; $\hat{A}$ is Hermitian and acts on a domain of $D$ qubits around the support of $\hat{h}[l]$, and is expanded as a sum of Pauli strings acting on the $D$ qubits,

[TABLE]

where $I$ denotes the index $i_{1}i_{2}\ldots i_{D}$. Define

[TABLE]

Our goal is to minimize the difference $||\Delta_{0}-\Delta||$. If the unitary $e^{-i\Delta\tau\hat{A}[l]}$ is defined over a sufficiently large domain $D$ (related to the correlation length of $|\Psi\rangle$, see Section A.3) then this error minimizes at $\sim 0$, for small $\Delta\tau$. Minimizing for real $a[l]$ corresponds to minimizing the quadratic function $f(a[l])$

[TABLE]

where

[TABLE]

whose minimum obtains at the solution of the linear equation

[TABLE]

In general, $\mathbf{S}+\mathbf{S}^{T}$ may have a non-zero null-space. Thus, we solve Eq. (13) either by applying the generalized inverse of $\mathbf{S}+\mathbf{S}^{T}$ or by an iterative algorithm such as conjugate gradient.

Note that the above results for geometric local Hamiltonians can be extended to Hamiltonians with long-range terms. The primary difference is a modification of the domain over which the unitary acts. For example, for a spin Hamiltonian with long-range pairwise terms, the scaled action of $e^{-\Delta\tau\hat{h}[l]}$ (if $\hat{h}[l]$ acts on qubits $i$ and $j$) can be emulated, to within accuracy $\varepsilon$, by a unitary constructed in the neighborhoods of $i$ and $j$, with the domain size given by the result in Eq. (34), see also discussion below).

For fermionic Hamiltonians, we replace the Pauli operators in Eq. (7) by fermionic field operators, i.e. $\sigma\in\{1,f,f^{\dagger},f^{\dagger}f\}$, and conjecture that the analogous result holds for the domains for the fermionic operators as for spin operators. For a number conserving Hamiltonian, such as the fermionic Hubbard Hamiltonian treated in Fig. 3 in the main text, we retain only those terms with an equal number of creation and annihilation operators, to conserve particle number.

A.2 Real Hamiltonians and states

As the cost to construct the quantities ${\bf{S}}$, ${\bf{b}}$ and to solve the linear system (13) increases exponentially with $D$, it is natural to seek out simplifications that can be made rigorously. For example, as described above, to emulate the non-unitary action of a number-conserving fermionic Hamiltonian, we can consider $\hat{A}[l]$ to contain only fermionic operator strings that also conserve particle number. In the main text, we also considered the case where $\hat{A}[l]$ is itself approximately $p$-local, which removes the exponential dependence on $D$.

Another common scenario concerns Hamiltonians and states that only have real matrix elements and coefficients in the $Z$ computational basis. Then, since $\langle\Psi|$ and $|\Delta_{0}\rangle$ are real in the computational basis and

[TABLE]

$b_{I}\equiv 0$ unless the matrix elements of $\hat{\sigma}_{I}^{\dagger}$ have non-zero imaginary part. When the Pauli basis is used, this means that $b_{I}\equiv 0$ unless $\hat{\sigma}_{I}^{\dagger}$ contains an odd number of $\hat{Y}$ operators.

The number of such operators for a domain size $D$, $y(D)$, is $2^{D}\frac{2^{D}-1}{2}\simeq\frac{4^{D}}{2}$. This can be shown by induction over $D$. For $D=1$, one has

[TABLE]

Pauli strings with an odd number of $\hat{Y}$’s (i.e. just $\hat{Y}$). For $D=2$, one has

[TABLE]

such strings, namely

[TABLE]

For $D=3$, the number grows to

[TABLE]

with strings

[TABLE]

and so on. Pauli strings of length $D+1$ containing an odd number of $\hat{Y}$ operators are obtained either by attaching a $\hat{Y}$ operator to a length-$D$ string containing an even number of $\hat{Y}$ operators, or by attaching a $\hat{I}$, $\hat{X}$, $\hat{Z}$ operator to a length-$D$ string containing an odd number of $\hat{Y}$ operators. Therefore, $O(D)$ obeys the recursion relation

[TABLE]

This recursion relation is solved, with the initial condition $y(D=1)=1$, by

[TABLE]

This result means that both $\mathbf{b}$ and $\mathbf{S}+{\mathbf{S}}^{T}$ can be assembled from $y(D)$ Pauli string expectation values, roughly half the number of measurements needed if one did not assume real Hamiltonians and states. Further, the dimension of $\mathbf{b}$ and $\mathbf{S}+{\mathbf{S}}^{T}$ is $y(D)$ and $y(D)\times y(D)$ respectively. Asymptotically, this reduces the cost of solving the linear system Eq. (13) by a factor of $1/8$, assuming dense matrix techniques.

A.3 Rigorous run time bounds

Here we present a more detailed analysis of the running time of the algorithm. Consider a $k$-local Hamiltonian

[TABLE]

acting on a $d$-dimensional lattice with $\|h_{i}\|\leq 1$, where $\|*\|$ is the operator norm.

Note that, if a quantum chemistry system is studied with an orthonormal basis of spatially localized states, such states can be approximately positioned on a lattice, and the results of this section apply on length scales larger than the size of the employed basis functions. In imaginary-time evolution one typically applies Trotter formulae to approximate

[TABLE]

for an initial state $|\Psi_{0}\rangle$ (which we assume to be a product state). This approximation leads to an error which can be made as small as one wishes by increasing the number of time steps $n$. Let $|\Psi_{s}\rangle$ be the state (after renormalization) obtained by applying $s$ terms $e^{-{\Delta\tau}\hat{h}[i]}$ from \big{(}e^{-{\Delta\tau}\hat{h}[1]}\ldots e^{-{\Delta\tau}\hat{h}[m]}\big{)}^{n}; with this notation $|\Psi_{mn}\rangle$ is the state given by Eq. (23). In the QITE algorithm, instead of applying each of the operators $e^{-{\Delta\tau}\hat{h}[i]}$ to $|\Psi_{0}\rangle$ (and renormalizing the state), one applies local unitaries $\hat{U}_{s}$ which should approximate the action of the original operator. Let $|\Phi_{s}\rangle$ be the state after $s$ unitaries have been applied.

Let $C$ be an upper bound on the correlation length of $|\Psi_{s}\rangle$ for every $s$: we assume that for every $s$, and every pair of observables $\hat{A}$ and $\hat{B}$ acting on domains separated by $\text{dist}(A,B)$ sites,

[TABLE]

Theorem 1**.**

For every $\varepsilon>0$, there are unitaries $\hat{U}_{s}$ each acting on

[TABLE]

qubits, such that

[TABLE]

Proof.

We have

[TABLE]

To bound the first term we use our assumption that the correlation length of $|\Psi_{s-1}\rangle$ is smaller than $C$. Consider a region $R_{v}$ of all sites that are a distance at most $v$ (in the Manhattan distance on the lattice) of the sites in which $h_{i_{s}}$ acts. Let $\text{tr}_{\backslash R_{v}}(|\Psi_{s}\rangle\langle\Psi_{s}|)$ be the reduced state on $R_{v}$, obtained by partial tracing over the complement of $R_{v}$ in the lattice. Since

[TABLE]

it follows from Eq. (24) and Lemma 9 of Brandão and Horodecki (2015) that

[TABLE]

where we used that for $n\geq 2\beta$,

[TABLE]

Above, $\|*\|_{1}$ is the trace norm. The key result in our analysis is Uhlmann’s theorem (see e.g. Lemmas 11 and 12 of Brandão and Horodecki (2015)). It states that two pure states with nearby marginals must be related by a unitary on the purifying system.

In more detail, if $|\eta\rangle_{AB}$ and $|\nu\rangle_{AB}$ are two states with partial traces $|\eta\rangle_{A}$ and $|\nu\rangle_{A}$ over the complement of $A$, such that $\||\eta\rangle_{A}-|\nu\rangle_{A}\|_{1}\leq\delta$, then there exists a unitary $\hat{V}$ acting on $B$ such that

[TABLE]

Applying Uhlmann’s theorem to $|\Psi_{s}\rangle$ and $|\Psi_{s-1}\rangle$, with $B=R_{v}$, and using Eq. (29), we find that there exists a unitary $\hat{U}_{s}$ acting on $R_{v}$ s.t.

[TABLE]

which by Eq. (A.3) implies

[TABLE]

Choosing $\nu=2C\ln(2\sqrt{2}nm\varepsilon^{-1})$ as the width of the support of the approximating unitaries, the error term above is $\varepsilon$. The support of the local unitaries is $k\nu^{d}$ qubits (as this is an upper bound on the number of qubits in $R_{\nu}$). Therefore each unitary $\hat{U}_{s}$ acts on at most

[TABLE]

qubits. ∎

Finding $\hat{U}_{s}$: In the algorithm we claim that we can find the unitaries $\hat{U}_{s}$ by solving a least-square problem. This is indeed the case if we can write them as $\hat{U}_{s}=e^{i\frac{\hat{A}[s]}{n}}$ with $\hat{A}[s]$ a Hamiltonian of constant norm. Then for sufficiently large $n$, $\hat{U}_{s}=I+i\frac{\hat{A}[s]}{n}+O\left(\frac{1}{n^{2}}\right)$ and we can find $\hat{A}[s]$ by performing tomography of the reduced state over the region where $\hat{U}_{s}$ acts and solving the linear problem given in the main text. Because we apply Uhlmann’s Theorem to

[TABLE]

using $e^{-{\Delta\tau}\hat{h}[i_{s}]}=I-{\Delta\tau}\hat{h}[i_{s}]+O\left(\frac{1}{n^{2}}\right)$ and following the proof of the Uhlmann’s Theorem, we find that the unitary can indeed be taken to be close to the identity, i.e. $\hat{U}_{s}$ can be written as $e^{i\frac{\hat{A}[s]}{n}}$.

Total Run Time: Theorem 1 gives an upper bound on the maximum support of the unitaries needed for a Trotter update, while tomography of local reduced density matrices gives a way to find the unitaries. The cost for tomography is quadratic in the dimension of the region, so it scales as $\exp(O(N_{q}))$. This is also the cost to solve classically the linear system which gives the associated Hamiltonian $\hat{A}[s]$ and of finding a circuit decomposition of $\hat{U}_{s}=e^{i\frac{\hat{A}[s]}{n}}$ in terms of two-qubit gates. As this is repeated $mn$ times, for each of the $mn$ terms of the Trotter decomposition, the total running time (of both quantum and classical parts) is

[TABLE]

This is exponential in $C^{d}$, with $C$ the correlation length, and quasi-polynomial in $n$ (the number of Trotter steps) and $m$ (the number of local terms in the Hamiltonian. Note that typically $m=O(N)$, with $N$ the number of sites). While this an exponential improvement over the $\exp(O(N))$ scaling classically, the quasi-polynomial dependence on $m$ can still be prohibitive in practice. Below we show how to improve on that.

Local Approximation: If one is only interested in a local approximation of the state (meaning that all the local marginals of $|\Phi_{nm}\rangle$ are close to the ones of $e^{-\beta\hat{H}}|\Psi_{0}\rangle$, but not necessarily the global states), then the support of the unitaries becomes independent of the number of terms of the Hamiltonian $m$ (while for global approximation we have a polylogarithmic dependence on $m$):

Theorem 2**.**

For every $\varepsilon>0$, there are unitaries $\hat{U}_{s}$ each acting on

[TABLE]

qubits such that, for every connected region $S$ of size at most $|S|$,

[TABLE]

Proof.

Consider the unitaries $\hat{U}_{s}$ obtained in the proof of Theorem 1 satisfying Eq. (32).

Consider the replacement of the local term of the Trotter expansion by the unitary $\hat{U}_{s}$ for all local terms which are more than $2C\log(1/\delta)$ sites away from the region $S$. Because the correlation length is always smaller than $C$, we find by Lemma 9 of Brandão and Horodecki (2015) that the total error $\eta$ on the reduced density matrix in region $S$ can be bounded as

[TABLE]

For the local terms which are at most a distance $2C\log(1/\delta)$ from the region $S$, in turn, the total error is bounded by the sum of each individual term, giving

[TABLE]

Choosing

[TABLE]

gives the result. ∎

Non-local Terms: Suppose the Hamiltonian has a term $\hat{h}[q]$ acting on qubits which are not nearby, e.g. on two sites $i$ and $j$. Then $e^{-{\Delta\tau}\hat{h}[q]}$ can still be replaced by a unitary, which only acts on sites $i$ and $j$ and qubits in the neighborhoods of the two sites. This is the case if we assume that the state has a finite correlation length and the proof is again an application of Uhlmann’s theorem (we follow the same argument from the proof of Theorem 1 but define $R_{v}$ in that case as the union of the neighborhoods of $i$ and $j$). Note however that the assumption of a finite correlation length might be less natural for models with long range interactions.

Scaling with temperature and increase of correlation length: Our discussion has been based on the assumption that the correlation length $C$ is small on all intermediate states. Here we discuss the range of validity of the assumption.

Let us begin with an example where the correlation length can increase very quickly with number of local terms applied (this was communicated to us by Guang Hao Low). Consider a projection on two qubits $\hat{P}_{i,i+1}=|0,0\rangle\langle 0,0|_{i,i+1}+|1,1\rangle\langle 1,1|_{i,i+1}$. Then

[TABLE]

with $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{(}2)$, is the GHZ state $(|0\ldots 0\rangle+|1\ldots 1\rangle)/\sqrt{2}$, which has correlation length $C=n$. While the projector $\hat{P}_{i,i+1}$ cannot appear as a local term $e^{-{\Delta\tau}\hat{h}[i]}$ in the Trotter decomposition, this example show that we cannot expect a speed-of-sound bound on the spread of correlations for a circuit with non-unitary gates; indeed the example shows a depth two circuit can already create long range correlations.

However, we expect that generically the correlations do grow ballistically. Consider the state

[TABLE]

after $n$ rounds have been applied. Let us assume the Hamiltonian acts on a line, is translation invariant and has nearest-neighbor interactions. Then the state is a matrix product state of bond dimension at most $2^{n}$. For matrix product states we can bound the correlations as follows (see e.g. Lemma 22 of Brandão and Horodecki (2015))

[TABLE]

where we define the gap of the matrix-product-state as $\Delta:=1-\lambda$, with $\lambda$ the second largest eigenvalue of the transfer matrix of the matrix product state (normalized so that the largest eigenvalue is one). In the GHZ example above, the gap $\Delta=0$ and that is the reason for the fast build up of correlations. Typically we expect the gap to be independent of $n$ or decrease mildly as $1/\text{poly}(n)$.

From the above, we can replace a non-unitary local Trotter term applied to $|\psi_{n}\rangle$ by an unitary acting on $O(n/\Delta)$ qubits. Taking $n=O(\beta)$ to reach temperature $\beta$ in the imaginary time evolution, the support of the unitaries would scale as $O(\beta/\Delta)$. Assuming $\Delta$ is a constant, we find a linear increase in temperature.

We also expect the linear growth of correlations/unitary support with inverse temperature also to hold generically in two dimensions, although there the analysis is more subtle as rigorous results for the expected behavior of the transfer operator (which becomes a one-dimensional tensor product operator) and its gap are not available.

A.4 Spreading of correlations

In the main text, we argued that the correlation volume $V$ of the state $e^{-\beta H}|\Psi\rangle$ is bounded for many physical Hamiltonians and saturates at the ground-state with $V\ll N$ where $N$ is the system size. To numerically measure correlations, we use the mutual information between two sites, defined as

[TABLE]

where $S(i)$ is the von Neumann entropy of the density matrix of site $i$ ($\rho(i)$) and similarly for $S(j)$, and $S(i,j)$ is the von Neumann entropy of the two-site density matrix for sites $i$ and $j$ ($\rho(i,j)$).

To compute the mutual information in Fig. 1 in the main text, we used matrix product state (MPS) and finite projected entangled pair state (PEPS) imaginary time evolution for the spin-$1/2$ 1D and 2D FM transverse field Ising model (TFI)

[TABLE]

where the sum over $\langle i,j\rangle$ pairs are over nearest neighbors. We use the parameter $h=1.25$ for the 1-D calculation and $h=3.5$ for the 2-D calculations as the ground-state is gapped in both cases. It is known that the ground-state correlation length is finite.

MPS. We performed MPS imaginary time evolution (ITE) on a 1-D spin chin with $L=50$ sites with open boundary conditions. We start from an initial state that is a random product state, and perform ITE using time evolution block decimation (TEBD) Vidal (2004b); Schollwöck (2005) with a first order Trotter decomposition. In this algorithm, the Hamiltonian is separated into terms operating on even and odd bonds. The operators acting on a single bond are exponentiated exactly. One time step is given by time evolution of odd and even bonds sequentially, giving rise to a Trotter error on the order of the time step $\Delta\tau$. In our calculation, a time step of $\Delta\tau=0.001$ was used.

We carry out ITE simulations with maximum bond dimension of $D=80$, but truncate singular values less than 1.0e-8 of the maximum singular value. In the main text, the ITE results are compared against the ground state obtained via the density matrix renormalization group (DMRG)). This should be equivalent to comparing to a long-time ITE ground state. The long-time ITE ($\beta=38.352$) ground state reached an energy per site of -1.455071, while the DMRG ground-state energy per site is -1.455076. The relative error of the nearest neighbor correlations is on the order of $10^{-4}$ to $10^{-3}$, and about $10^{-2}$ for correlations between the middle site and the end sites (a distance of 25 sites). The error in fidelity between the two ground states was about $5\times 10^{-4}$.

PEPS. We carried out finite PEPS Nishino and Okunishi (1996); Verstraete and Cirac (2004); Verstraete et al. (2006); Orús (2014) imaginary time evolution for the two-dimensional transverse field Ising model on a lattice size of $21\times 31$. The size was chosen to be large enough to see the spread of mutual information in the bulk without significant effects from the boundary. The mutual information was calculated along the long (horizontal) axis in the center of the lattice. The standard Trotterized imaginary time evolution scheme for PEPS Verstraete et al. (2008) was used with a time step $\Delta\tau=0.001$, up to imaginary time $\beta=6.0$, starting from a random product state. To reduce computational cost from the large lattice size, the PEPS was defined in a translationally invariant manner with only 2 independent tensors Jordan et al. (2008) updated via the so-called “simple update” procedure Jiang et al. (2008). The simple update has been shown to be sufficiently accurate for capturing correlation functions (and thus $I(i,j)$) for ground states with relatively short correlation lengths (compared to criticality) Lubasch et al. (2014a, b). We chose a magnetic field value $h=3.5$ which is detuned from the critical field ($h\approx 3.044$) but still maintains a correlation length long enough to see interesting behavior.

Accuracy: Even though the simple update procedure was used for the tensor update, we still needed to contract the $21\times 31$ PEPS at at every imaginary time step $\beta$ for a range of correlation functions, amounting to a large number of contractions. To control the computational cost, we limited our bond dimension to $D=5$ and used an optimized contraction scheme Xie et al. (2017), with maximum allowed bond dimension of $\chi=60$ during the contraction. Based on converged PEPS ground state correlation functions with a larger bond dimension of $D=8$, our $D=5$ PEPS yields $I(i,i+r)$ (where $r$ denotes horizontal separation) at large $\beta$ with a relative error of $\approx 1\%$ for $r=1-4$, $5\%$ or less for $r=5-8$, and $10\%$ or greater for $r>8$. At smaller values of $\beta$ ($<0.5$) the errors up to $r=8$ are much smaller because the bond dimension of 5 is able to completely support the smaller correlations (see Fig. 1, main text). While error analysis on the 2D Heisenberg model Lubasch et al. (2014a) suggests that errors with respect to $D=\infty$ may be larger, such analysis also confirms that a $D=5$ PEPS captures the qualitative behavior of correlation in the range $r=5-10$ (and beyond). Aside from the bond dimension error, the precision of the calculations is governed by $\chi$ and the lattice size. Using the $21\times 31$ lattice and $\chi=60$, we were able to converge entries of single-site density matrices $\rho(i)$ to a precision of $\pm 10^{-6}$ (two site density matrices $\rho(i,j)$ had higher precision). For $\beta=0.001-0.012$, the smallest eigenvalue of $\rho(i)$ fell below this precision threshold, leading to significant noise in $I(i,j)$. Thus, these values of $\beta$ are omitted from Fig. 1 (main text) and the smallest reported values of $I$ are $10^{-6}$, although with more precision we expect $I\to 0$ as $r\to\infty$.

Finally, the energy and fidelity errors were computed with respect to the PEPS ground state of the same bond dimension at $\beta=10.0$ (10000 time steps). The convergence of the quantities shown in Fig. 1 (main text) thus isolates the convergence of the imaginary time evolution, and does not include effects of other errors that may result from deficiencies in the wavefunction Ansatz.

A.5 Comparison to classical algorithms

In the main text, we noted that QITE provided an exponential speedup per iteration over the direct classical implementation of imaginary time evolution algorithm, given bounded correlation length $C$ during the evolution. We now compare to some other possible classical algorithms.

We first note that a finite correlation length $C_{0}$ in the ground-state does not itself imply an efficient classical strategy. For example, a simple heuristic is to solve the problem locally, e.g. to truncate the problem size at correlation length $C_{0}$ of the ground-state and solve by exact diagonalization, which can be done in time $\exp(O(C_{0}d))$ in $d$ spatial dimensions. But this will not generally converge to the correct ground-state in a frustrated Hamiltonian, as this would efficiently solve NP-hard classical satisfiability problems even though these have $C_{0}=0$; physical examples include glassy models.

Similarly, as QITE defines a quantum circuit for the imaginary time evolution, we might attempt to use it for a faster classical simulation. If we are only interested in local observables, we can apply the circuit in the Heisenberg picture in a classical emulation. However, this gives an extra exponential dependence on the number of previous time-steps: after the unitaries associated to $(e^{-\Delta\tau\hat{h}[1]}e^{-\Delta\tau\hat{h}[2]}\dots)^{l}$ have been applied, the cost of applying the next unitary scales as $\exp(O(lD))$, with $D$ the domain size of the unitaries, instead of $\exp(O(D))$ in QITE.

Alternatively, if $|\Psi\rangle$ is represented by a tensor network in a classical simulation, then $e^{-\Delta\tau\hat{h}[l]}|\Psi\rangle$ can be represented as a classical tensor network with increased bond dimension Vidal (2004a); Schollwöck (2011). However, the bond dimension will scale as $\exp(O(lD))$. Further, apart from the extra exponential dependence on $l$, another potential drawback in the tensor network approach is that we cannot guarantee contracting the resulting classical tensor network for an observable is efficient; it is a #P-hard problem in the worst case in 2D (and even in the average case for Gaussian distributed tensors) Vidal (2004a); Schollwöck (2011); Schuch et al. (2007); Haferkamp et al. (2018).

A.6 Simulation models

We here define, and give some background on, the models used in the QITE and QLanczos simulations.

A.6.1 1 qubit field model

[TABLE]

This Hamiltonian has previously been used as a model for quantum simulations on physical devices in Ref. Lamm and Lawrence (2018).We used $\alpha=\frac{1}{\sqrt{2}}$ and $\beta=\frac{1}{\sqrt{2}}$. In simulations with this Hamiltonian, the qubit is assumed to be initialized in the $Z$ basis.

A.6.2 1D Heisenberg and transverse field Ising model

The 1D short-range Heisenberg Hamiltonian is defined as

[TABLE]

the 1D short-range Heisenberg Hamiltonian in the presence of a field as

[TABLE]

the 1D long-range Heisenberg Hamiltonian as

[TABLE]

and the 1D AFM transverse-field Ising Hamiltonian as

[TABLE]

A.6.3 1D Hubbard model

The 1D Hubbard Hamiltonian is defined as

[TABLE]

where $\hat{n}_{i\sigma}=a^{\dagger}_{i\sigma}a_{i\sigma}$, $\sigma\in\{\uparrow,\downarrow\}$, and $\langle\cdot\rangle$ denotes summation over nearest-neighbors, here with open-boundary conditions. We label the $n$ lattice sites with an index $i=0\dots n-1$, and the $2n-1$ basis functions as $|\varphi_{0}\rangle=|0\uparrow\rangle$, $|\varphi_{1}\rangle=|0\downarrow\rangle$, $|\varphi_{2}\rangle=|1\uparrow\rangle$, $|\varphi_{3}\rangle=|1\downarrow\rangle$ $\dots$. Under Jordan-Wigner transformation, recalling that

[TABLE]

with $p=0\dots 2n-2$ and $q<p$, the Hamiltonian takes the form

[TABLE]

A.6.4 H2 molecule minimal basis model

We use the hydrogen molecule minimal basis model at the STO-6G level of theory. This is a common minimal model of hydrogen chains Hachmann et al. (2006); Motta et al. (2017) and has previously been studied in quantum simulations, for example in O’Malley et al. (2016). Given a molecular geometry (H-H distance $R$) we perform a restricted Hartree-Fock calculation and express the second-quantized Hamiltonian in the orthonormal basis of RHF molecular orbitals as Szabo and Ostlund (1996)

[TABLE]

where $a^{\dagger}$, $a$ are fermionic creation and annihilation operators for the molecular orbitals. The Hamiltonian (55) is then encoded by a Bravyi-Kitaev transformation into the 2-qubit operator

[TABLE]

with coefficients $g_{i}$ given in Table I of O’Malley et al. (2016).

A.6.5 MAXCUT Hamiltonian

The MAXCUT Hamiltonian encodes the solution of the MAXCUT problem. Given a graph $\Gamma=(V,E)$, where $V$ is a set of vertices and $E\subseteq V\times V$ is a set of links between vertices in $V$, a cut of $\Gamma$ is a subset $S\subseteq V$ of $V$. The MAXCUT problem consists in finding a cut $S$ that maximizes the number of edges between $S$ and $S^{c}$ (the complement of $S$). We denote the number of links in a given cut $S$ as $C(S)$.

In Figure 3 of the main text, we consider a graph $\Gamma$ with vertices and links

[TABLE]

respectively. It is easy to verify that $S=\{0,2,4\}$, $\{0,1,2\}$, $\{3,4\}$ and their complements $S^{c}$ are solutions of the MAXCUT problem, with weight $C_{max}=5$.

The MAXCUT problem can be formulated as a Hamiltonian ground-state problem, by (i) associating a qubit to every vertex in $V$, (ii) associating to every partition $S=$ an element of the computational basis (here assumed to be in the $z$ direction) of the form $|z_{0}\dots z_{n-1}\rangle$, where $z_{i}=1$ if $i\in S$ and $z_{i}=0$ if $i\in S^{c}$, and finding the minimal (most negative) eigenvalue of the $2$-local Hamiltonian

[TABLE]

The spectrum of $\hat{C}$ is a subset of numbers $C\in\{0,1\dots|E|\}$.

In the present work, we initialize the qubits in the state $|\Phi\rangle=|+\rangle^{\otimes n}$, where $|+\rangle=\frac{|0\rangle+|1\rangle}{\sqrt{2}}$, and evolve $\Phi$ in imaginary time. Measuring the evolved state at time $\beta$ $|\Phi(\beta)\rangle$ will collapse it onto an element $|z_{0}\dots z_{n-1}\rangle$ of the computational basis, which is also an eigenfunction of $\hat{C}$ with eigenvalue $C$. In Figure 3 in the main text, we illustrate the probability $P(|C|=C_{max})$ that such measurements yield a MAXCUT solution. Note that, even in the presence of oscillations (with the smallest domain size $D=2$) this probability remains above $60\%$.

A.7 Numerical simulation details

A.7.1 QITE stabilization

Sampling noise in the expectation values of the Pauli operators can affect the solution to Eq. (13) that sometimes lead to numerical instabilities. We regularize $\mathbf{S}+\mathbf{S}^{T}$ against such statistical errors by adding a small $\delta$ to its diagonal. To generate the data presented in Figures 2 and 4 of the main text, we used $\delta=0.01$ for 1-qubit calculations and $\delta=0.1$ for 2-qubit calculations.

A.7.2 QLanczos stabilization

In quantum Lanczos, we generate a set of wavefunctions for different imaginary-time projections of an initial state $|\Psi\rangle$, using QITE as a subroutine. The normalized states are

[TABLE]

where $n_{l}$ is the normalization constant. For the exact imaginary-time evolution and $l$, $l^{\prime}$ both even (or odd) the matrix elements

[TABLE]

can be computed in terms of expectation values (i.e. experimentally accessible quantities) only. Indeed, defining $2r=l+l^{\prime}$, we have

[TABLE]

and similarly

[TABLE]

The quantities $n_{r}$ can be evaluated recursively, since

[TABLE]

For inexact time evolution, the quantities $n_{r}$ and $\langle\Phi_{r}|\hat{H}|\Phi_{r}\rangle$ can still be used to approximate $S_{l,l^{\prime}}$, $H_{l,l^{\prime}}$.

Given these matrices, we then solve the generalized eigenvalue equation $\mathbf{H}\mathbf{x}=E\mathbf{S}\mathbf{x}$ to find an approximation to the ground-state $|\Phi^{\prime}\rangle=\sum_{l}x_{l}|\Phi_{l}\rangle$ for the ground state of $\hat{H}$. This eigenvalue equation can be numerically ill-conditioned, as $S$ can contain small and negative eigenvalues for several reasons (i) as $m$ increases the vectors $|\Phi_{l}\rangle$ become linearly dependent; (ii) simulations have finite precision and noise; (iii) $S$, $H$ are computed approximately when inexact time evolution is performed.

To regularize the problem, out of the set of time-evolved states we extract a better-behaved sequence as follows (i) start from $|\Phi_{\text{last}}\rangle=|\Phi_{0}\rangle$ (ii) add the next $|\Phi_{l}\rangle$ in the set of time-evolved states s.t. $|\langle\Phi_{l}|\Phi_{\text{last}}\rangle|<s$, where $s$ is a regularization parameter $0<s<1$ (iii) repeat, setting the $|\Phi_{\text{last}}\rangle=\Phi_{l}$ (obtained from (ii)), until the desired number of vectors is reached. We then solve the generalized eigenvalue equation $\tilde{\mathbf{H}}\mathbf{x}=E\tilde{\mathbf{S}}\mathbf{x}$ spanned by this regularized sequence, removing any eigenvalues of $\tilde{\mathbf{S}}$ less than a threshold $\epsilon$. The exact emulated QLanczos calculations reported in the main text were stabilized with this algorithm (the source of error here is primarily (iii)) using stabilization parameter $s=0.95$ and $\epsilon=10^{-14}$. The stabilization parameters used in the QVM and QPU QLanczos calculations were $s=0.75$ and $\epsilon=10^{-2}$ (the main source of error in the simulations was (ii)). Note that the stabilization procedure is unlikely to fix all possible numerical instabilities, but was sufficient for all models and calculations performed in this work.

A.8 METTS algorithm

The METTS (minimally entangled typical thermal state) algorithm Stoudenmire and White (2010b); White (2009b) is a sampling method to calculate thermal properties based on imaginary time evolution. Consider the thermal average of an observable $\hat{O}$

[TABLE]

where $\{|i\rangle\}$ is an orthonormal basis set, and $Z$ is the partition function. Defining $|\phi_{i}\rangle=P_{i}^{-1/2}e^{-\frac{\beta}{2}\hat{H}}|i\rangle$ with $P_{i}=\langle i|e^{-\beta\hat{H}}|i\rangle$, we obtain

[TABLE]

The summation in Eq. (65) can be estimated by sampling $|\phi_{i}\rangle$ with probability $P_{i}/Z$, and summing the sampled $\langle\phi_{i}|\hat{O}|\phi_{i}\rangle$.

In standard Metropolis sampling for thermal states, one starts from $|\phi_{i}\rangle$ and obtains the next state $|\phi_{j}\rangle$ from randomly proposing and accepting based an acceptance probability. However, rejecting and resetting in the quantum analog of Metropolis Temme et al. (2011b) is complicated to implement on a quantum computer, requiring deep circuits. The METTS algorithm provides an alternative way to sample $|\phi_{i}\rangle$ distributed with probability $P_{i}/Z$ without this complicated procedure. The algorithm is as follows:

Choose a classical product state (PS) $|i\rangle$. 2. 2.

Compute $|\phi_{i}\rangle=P_{i}^{-1/2}e^{-\frac{\beta}{2}\hat{H}}|i\rangle$ and calculate observables of interest. 3. 3.

Collapse the state $|\phi_{i}\rangle$ to a new PS $|i^{\prime}\rangle$ with probability $p(i\rightarrow i^{\prime})=|\langle i^{\prime}|\phi_{i}\rangle|^{2}$ and repeat Step 2.

In the above algorithm, $|\phi_{i}\rangle$ is named a minimally entangled typical thermal state (METTS). One can easily show that the set of METTS sampled following the above procedure has the correct Gibbs distribution Stoudenmire and White (2010b). Generally, $\{|i\rangle\}$ can be any orthonormal basis. For convenience when implementing METTS on a quantum computer, $\{|i\rangle\}$ are chosen to be product states. On a quantum emulator or a quantum computer, the METTS algorithm is carried out as following:

Prepare a product state $|i\rangle$. 2. 2.

Imaginary time evolve $|i\rangle$ with the QITE algorithm to $|\phi_{i}\rangle=P_{i}^{-1/2}e^{-\frac{\beta}{2}\hat{H}}|i\rangle$, and measure the desired observables. 3. 3.

Collapse $|\phi_{i}\rangle$ to another product state by measurement.

In practice, to avoid long statistical correlations between samples, we used the strategy of collapsing METTS onto alternating basis sets Stoudenmire and White (2010b). For instance, for the odd METTS steps, $|\phi_{i}\rangle$ is collapsed onto the $X$-basis (assuming a $Z$ computational basis, tensor products of $|+\rangle$ and $|-\rangle$), and for the even METTS steps, $|\phi_{i}\rangle$ is collapsed onto the $Z$-basis (tensor products of $|0\rangle$ and $|1\rangle$). The statistical error is then estimated by block analysis Flyvbjerg and Petersen (1989).

A.9 Implementation on emulator and quantum processor

We used pyQuil, an open source Python library, to express quantum circuits that interface with both Rigetti’s quantum virtual machine (QVM) and the Aspen-1 quantum processing units (QPUs).

pyQuil provides a way to include noise models in the QVM simulations. Readout error can be included in a high-level API provided in the package and is characterized by $p_{00}$ (the probability of reading $|0\rangle$ given that the qubit is in state $|0\rangle$) and $p_{11}$ (the probability of reading $|1\rangle$ given that the qubit is in state $|1\rangle$). Readout errors can be mitigated by estimating the relevant probabilities and correcting the estimated expectation values. We do so by using a high level API present in pyQuil. A general noise model can also be applied to a gate in the circuit by applying the appropriate Kraus maps. Included in the package is a high level API that applies the same decoherence error attributed to energy relaxation and dephasing to every gate in the circuit. This error channel is characterized by the relaxation time $T_{1}$ and coherence time $T_{2}$. We also include in our emulation our own high-level API that applies the same depolarizing noise channel to every single gate by using the appropriate Kraus maps. The depolarizing noise is characterized by $p_{1}$, the depolarizing probability for single-qubit gates and $p_{2}$, the depolarizing probability for two-qubit gates. We do not include all sources of error in our emulation. We applied the same depolarizing and dephasing channels to each gate operation for all qubits, when in reality, they can vary from qubit to qubit. In addition, noise due to cross-talk between qubits cannot be modeled using the QVM and is another source of discrepancy between the QVM and QPU results.

We investigate the influence of noise on the 2-qubit results obtained via the QVM using different noise parameters;

• •

Noise model 1:

[TABLE]

• •

Noise model 2:

[TABLE]

• •

Noise model 3:

[TABLE]

Noise model 1 reflects realistic parameters that characterize the Aspen-1 QPUs we run our calculations on; $p_{00}$, $p_{11}$, $T_{1}$, and $T_{2}$ are reported values whereas $p_{1}$ and $p_{2}$ are values typically used to benchmark error mitigation algorithms Temme et al. (2017). We repeated 10 calculations for each noise model and note there is practically no variation from run to run. Fig. 5(a) shows that reducing the readout error does not greatly affect the converged ground state energy after readout error mitigation has been performed. However, reducing the other sources of error does improve the converged energy. Note that sufficient measurement samples are used such that the sampling variance is smaller than that due to noise.

We also ran 2-qubit simulations on different pairs of qubits on Aspen-1, with Q1 consisting of qubits 14, 15 and Q2 consisting of qubits 0,1. These two pairs are reported to have different noise characteristics,

• •

Q1:

[TABLE]

• •

Q2:

[TABLE]

Based on this, we expect simulations on Q2 to be worse. Note that in contrast to our QVM calculations, the results from the actual devices varied from run to run. Thus, we present the mean and standard deviation for 10 different runs on each pair. (Similarly, sufficient samples are taken when running the QVM such that the sampling variance is smaller than that due to noise). Fig. 5(b) indeed demonstrates that Q2 provides a less faithful implementation of the quantum algorithm.

A.9.1 Parameters used in QVM and QPUs simulations

In this section, we include the parameters used in our QPU and QVM simulations. Note that all noisy QVM simulations (unless stated otherwise in the text) were performed with noise parameters from noise model 1. We also indicate the number of samples used during measurements for each Pauli operator.

A.10 Comparison of QITE and VQE

To address the feasibility of running QITE for larger systems on near-term devices, we compared the total number of Pauli string measurements needed for both VQE and QITE to obtain the ground state of two different spin models; (a) a 1D Heisenberg chain in a magnetic field with the parameters $J=B=1$; the 4-site instance of this model was studied in Ref. Kandala et al. (2017b), and (b) 1D AFM transverse-field Ising model ($J=h=1/\sqrt{2}$). Specifically, we estimated how many expectation values of Pauli strings would need to be measured to obtain the ground state of a 4-site and 6-site instance. The state under the evolution of QITE or in VQE was said to be converged to the ground state if its energy was within 1$\%$ of the exact ground state energy for the Heisenberg model, and 1$\%$ or 2$\%$ for the Ising model (the relaxed criterion for the Ising model was chosen so that the VQE optimization could complete in a reasonable number of steps). The following section describes how we counted the total number of Pauli string measurements in VQE and QITE.

A.10.1 Counting Pauli strings in VQE

To perform the VQE calculations, we used the hardware-efficient variational Ansatz as described in  Kandala et al. (2017b). This consists of first applying rotation unitaries represented by $U^{q,i}(\boldsymbol{\theta})=Rz_{\theta^{q,i}_{1}}Rx_{\theta^{q,i}_{2}}Rz_{\theta^{q,i}_{3}}$ to all qubits before applying layers of a certain depth $d$; each layer begins by applying $CZ$ gates between nearest-neighbors followed by applying $U^{q,i}(\boldsymbol{\theta})=Rz_{\theta^{q,i}_{1}}Rx_{\theta^{q,i}_{2}}Rz_{\theta^{q,i}_{3}}$ to all qubits again. Details of the circuit can be seen in Fig. 6.

As in  Kandala et al. (2017b), we also used the simultaneous perturbation stochastic approximation (SPSA) algorithm as the optimization protocol. The SPSA algorithm is commonly used because (i) it performs well in the presence of stochastic fluctuations and (ii) it requires only evaluating the objective function twice to update the variational parameters regardless of the number of parameters involved. The performance of the optimizer depends on the hyperparameters $\alpha$ and $\gamma$ as described in  Kandala et al. (2017b) and we found that their reported values of $\alpha=0.602$ and $\gamma=0.101$ also gave the best results for us. Numerical evidence of this is provided later on.

Evaluating the objective function involves estimating the expectation value of the Pauli strings that appear in the Hamiltonian. To prevent sampling errors from influencing the comparison, we evaluated the expectation values exactly. We conducted the VQE calculations using Qiskit, a quantum emulator Python package provided by IBM. The package provides both the SPSA algorithm and a variational Ansatz which we modified to reproduce the exact Ansatz used in Kandala et al. (2017b).

To count the number of Pauli strings needed for convergence, we ran VQE using different layer depths and determined the number of iterations $N$ needed for the algorithm to converge to a state with an energy within a certain percentage (1$\%$ or 2$\%$) of the exact ground state energy. Examples of converged VQE calculations for the 1D Heisenberg model are given in Fig. 7(a) and (b). In each iteration, the objective function was evaluated twice, and the evaluation of the objective function required measuring the expectation value of the $M$ Pauli strings that appear in the Hamiltonian. Therefore, the total number of Pauli strings $P_{total}$ is given as

[TABLE]

We note that the results from one VQE trajectory can differ slightly from the next. For our VQE calculations, we always performed 10 trajectories and analyzed our data using the average trajectory. We summarize the VQE parameters that we found gave the lowest number of total Pauli measurements to converge to the ground state in Table  7. For the 6-site 1D AFM transverse field Ising model, VQE could not converge to within 1$\%$ and we instead used the 2$\%$ convergence criterion. We found that for $\alpha=0.602,\gamma=0.101$, our simulation results for the 6-site 1D Heisenberg model indicates that using a circuit depth of 20 requires the least number of total Pauli measurements. We also ran some tests to determine what values of $\alpha$ and $\gamma$ gave the best result for the 6-site 1D Heisenberg model; the VQE calculation for the 6-site model conducted using $\alpha=0.602$, $\gamma=0.101$ and a circuit depth of 20 converged within 8400 optimization steps. We ran VQE calculations for different $\alpha$ and $\gamma$ using the same circuit depth of 20 and a total of 9000 optimization steps. The data in table  8 clearly shows that $\alpha=0.602$ and $\gamma=0.101$ produced the best result for us.

A.10.2 Counting Pauli strings in QITE

To implement QITE, we used the second-order Trotter decomposition given by

[TABLE]

to carry out the real time evolution. We initialized our state as: (a) $|0101\ldots\rangle$ for the 1D Heisenberg model and (b) maximally-mixed state for the 1D AFM transverse-field Ising model. We converged to the ground state using a time step of ${\Delta\tau}=0.1$ and a domain size $D$ of 4, as seen in Figs. 7(c) and (d). To count the number of Pauli strings, we note that a domain size of 4 implies that to evaluate $e^{-{\Delta\tau}/2\hat{h}[i]}$ involves measuring $4^{4}=256$ Pauli strings (without using the real-valued nature of the Hamiltonian). Therefore, with a total number of Trotter steps $T$, the total number of Pauli strings $P_{total}$ is given as

[TABLE]

We summarize the parameters that we used for QITE to obtain the ground state in table  9. We had no trouble converging our ground state to arbitrary accuracy using QITE but we used the same convergence criterion as for VQE to facilitate comparison.

A.10.3 VQE and QITE

Data from Table  10 suggests that QITE is competitive with VQE with respect to the number of Pauli string measurements. In fact, for the 6 qubit system, the number of measurements needed in QITE was significantly less than in VQE, due largely to the SPSA iterations needed to reach convergence when optimizing the VQE energy. While it is likely that the VQE costs could be lowered by using a better optimizer, or a better VQE Ansatz, we also note that the counts for QITE can also be reduced by using the methods outlined in the main text and earlier sections that discussed how one can economize measurements in QITE. The widespread current implementation of VQE and the observed performance of QITE suggest that it will be practical to implement the QITE protocol for intermediate system sizes on near-term devices.