Quantum computation of thermal averages in the presence of a sign problem

Giuseppe Clemente, Marco Cardinali, Claudio Bonati, Enrico Calore, Leonardo Cosmai, Massimo D’Elia, Alessandro Gabbana, Davide Rossini, Fabio Sebastiano Schifano, Raffaele Tripiccione, and Davide Vadacchino (QuBiPF Collaboration)∗ Dipartimento di Fisica dell’Università di Pisa and INFN Sezione di Pisa, Largo Pontecorvo 3, I-56127 Pisa, Italy. Università degli Studi di Ferrara and INFN Sezione di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy INFN Sezione di Bari, I-70126 Bari, Italy


I. INTRODUCTION
The computation of thermal averages for a quantum mechanical system is a common challenge to many different fields of modern physics, in diverse areas such as condensed matter or fundamental interactions. In most cases, the issue is connected to the determination of the phase diagram of the system: this is the case, for instance, of Quantum Chromodynamics (QCD), the theory which describes strong interactions within the standard model of particle physics. In this context, the QCD phase diagram is related to fundamental questions regarding the properties of the Universe in its early stages, or the internal structure and composition of compact astrophysical objects, such as neutron stars [1,2].
In the absence of an exact solution, the strongly interacting nature of the system generally makes analytical tools, such as perturbation theory, unfeasible, and demands for a numerical approach. Monte-Carlo (MC) algorithms represent a possible solution to the problem: in a nutshell, the computation of quantum thermal averages is brought to the form of a classical computation. Namely, after properly choosing the basis of the Hilbert space (the so-called computational basis), one performs a classical computation by means of classical MC methods, when possible. The path-integral formulation of quantum mechanics provides a suitable framework where this can be done systematically. The partition function of the quantum thermal system is rewritten in terms of an Euclidean path integral, which resembles the partition function of a classical statistical system in one more dimension (Euclidean time); Euclidean paths, or field configurations, are assigned a Boltzmann-like weight, corresponding to the exponential of minus the action. The path integral, in particular thermal averages, are then computed using standard tools, such as MC importance sampling.
Unfortunately, the prescription depicted above fails dramatically when the weight turns out to be nonpositive defined or, even worse, complex. This constitutes the infamous sign problem, which prevents progress in many fields, such as the accurate determination of the QCD phase diagram at finite baryon density [3][4][5], and can be shown to be a NP-hard task [6]. The problem is usually related to the path-integral formulation itself, or, in other words, to the choice of the computational basis. Indeed, in the canonical formulation of the partition function, all the Hamiltonian eigenstates have a positive weight, as it stems from the hermiticity of the Hamiltonian itself, or of other possible conserved charges coupled to chemical potentials. However, determining the spectrum and the eigenstates of the Hamiltonian usually implies a similar or even greater level of computational difficulty.
The sign problem is often taken as one of the compelling motivations for the development of Quantum Computing (QC). Moving forward from the original proposal of Feynman [7], QC has been nowadays exploiting the possibility to operate algorithms based on quantum logic gates, aimed at either solving certain tasks more efficiently than classically, or in general reproducing the real-time evolution of a system on a digital, or analog, platform [8,9]. Cornerstone examples are the Lloyd's proposal for the quantum simulation of the Schrödinger equation [10], the Shor's algorithm for integer factorization [11], and the Grover's algorithm for searching from an unstructured database [12]. Our current investigation revolves around the question on how one can make use of a quantum computer in order to reproduce the quantum thermal averages and then explore the phase diagram, i.e., on how the sign problem can be solved in practice.
The reconstruction of quantum thermal averages via a quantum computer has been object of a rather intense theoretical activity in the past few years [13][14][15][16][17][18][19][20][21][22][23][24][25]. Here we focus on a specific implementation, proposed in Ref. [16], of the so-called Quantum Metropolis Sampling (QMS), where the main idea is to exploit quantum parallelism in order to sample directly in the eigenbasis of the Hamiltonian. The purpose of this study is to implement the algorithm for a prototype system where one finds a sign problem in the path-integral formulation, showing that this approach completely solves the problem and investigating the possible systematics of the method. In particular, we consider a system of three quantum spins coupled to each other via antiferromagnetic interactions [26,27], a model that has already been used as a QC testbed in several works, and in particular in quantum simulators using trapped ions [28,29] and nuclear magnetic resonance [30]. While it is almost trivial to diagonalize the system Hamiltonian and perform a classical sign problem-free computation, the way in which the quantum computer solves it is general and applies, at least in principle, to quantum systems of arbitrary complexity. Indeed, as we will better explain in the following, everything relies on the possibility of performing the quantum phase estimation (QPE) algorithm for the system at hand: it is at this point that the quantum advantage becomes manifest, since the algorithm is able to do that with no need for an explicit diagonalization of the Hamiltonian. This paper is organized as follows. In Sec. II we give a brief description of the model, showing how a discrete path integral (Trotter-Suzuki decomposition [31,32]) performed in the standard computational basis leads to a sign problem. Section III provides a short summary of the QMS algorithm adopted in this study. Section IV contains a summary of the results obtained with our model, while the possible systematics related to the algorithm are discussed in Sec. V. We conclude with Sec. VI, summarizing our results and discussing possible extensions to more complex systems.
The quantum algorithm illustrated in this study has been implemented on a quantum simulator (Simulator for Universal Quantum Algorithms, SUQA) developed on purpose by one of the authors (GC), which is inspired to the well known Qiskit simulator 1 [33].

II. DESCRIPTION OF THE MODEL
In this section we introduce the system to be used as a playground for testing the QMS algorithm. It is made up of three quantum spin-1/2 variables, with frustrated pair interactions. Its Hamiltonian reads 1 We used our own simulator to add some ad-hoc features regarding the implementation of hybrid quantum-classical operations.
where 1 stands for the 2 × 2 identity operator, σ j (with j = x, y, z) denotes the usual spin-1/2 Pauli matrices, and J is a positive parameter (in order to make the couplings antiferromagnetic). Each of the 1 or σ x operators acts on the spin corresponding to its position in the tensor product, and the full Hilbert space of the system is spanned by eight basis states. A possible choice, which makes the problem trivial, is to select as a basis the eight products of eigenstates (for each spin) of the σ x operator: this is also a basis of eigenstates of the Hamiltonian, which has only two energy levels, the lower one with energy E 0 = −J and degeneracy 6 (two spins aligned in either direction and opposite to the third spin), the higher one with energy E 1 = 3J and degeneracy 2 (three spins aligned in either direction).
We are however interested in working with the standard computational basis, made up of eigenstates of the σ z operators for all spins, which is the one where the sign problem appears. In the following we will label, for simplicity, each of the eight states of this basis by its decimal correspondent, i.e. starting from |0 ≡ |000 to |7 ≡ |111 . In this basis, the partition function of the system reads where β is the inverse temperature (we adopt natural units, = k B = 1). The path-integral approach consists in rewriting the expression above by partitioning the exponential into the product of N smaller "evolution steps" (Trotter-Suzuki decomposition [31,32]), and inserting N − 1 projectors over a complete set of states between each couple of evolution steps. Doing so, we arrive to the well known expression where α N +1 ≡ α 1 and ({α i }) is the space of all possible configurations, which we can interpret as periodic paths labelled by a discrete Euclidean time variable i. If the product of matrix elements turns out to be positive for each path, then the sum over all possible Euclidean paths appearing in Eq. (3) can be sampled by a classical MC algorithm. However, it is easy to show that this is not the case for this system, i.e. that it is always possible to find paths for which the product of matrix elements turns out to be negative. In order to show this fact, let us rewrite the elementary evolution step in the following form: where ε ≡ βJ/N [details about the derivation of Eq. (4) are reported in Appendix A]. Since the Hamiltonian is composed of a symmetric combination of pair-products of σ x operators, which flip the σ z -component of the corresponding spin, it should be clear that the only non-zero non-diagonal matrix elements of exp(−βH/N ) are those between states differing by exactly two spin flips. This fact enables to split the computational basis into two disjoint sets: Non-diagonal elements within the same set are all equal and negative, and in particular their common value is −e ε (1 − e −4ε )/4, while non-diagonal elements across the two sets are exactly zero. In contrast, diagonal elements are all equal to e ε (3 + e −4ε )/4.
To summarize, the discrete path-integral sum in Eq. (3) simplifies by decoupling into two distinct and simpler sums, where paths run either within Set 1 or within Set 2. However, within each sum, paths presenting an odd number 2 of flips between different states carry a negative weight. The contribution of these paths cannot be rewritten in terms of a positive Boltzmann weight, i.e. in the form ] is a real Euclidean action. This invalidates the possibility of using classical methods which approximate the sum by means of a MC sampling of its terms.
In the QMS algorithm described in the next section, the focus is brought again to the basis of eigenstates of the Hamiltonian, which are sampled according to their thermal weight. However this is not done by exploiting the exact solution of the problem (that, in this case, would be trivial), but in a way which can be easily generalized to more complex cases. As we shall see below, the key ingredient to achieve this task will be the QPE algorithm [34,35].

III. THE QUANTUM METROPOLIS SAMPLING ALGORITHM
This section is a brief, yet self-contained, overview of the main steps of the QMS algorithm, to be operated on a quantum computer. Further details can be found in the original publication, Ref. [16].
The key ingredient of QC, which potentially enables an exponential speedup over any known classical algorithm analogue, is quantum parallelism. Namely, the possibility, for a set of N two-level quantum systems (the socalled qubits, quantum analogues of classical bits), to be in a superposition of states, thus simultaneously acting on an exponential number 2 N of classical configurations. Such an ensemble of qubits can be suitably manipulated through a sequence of invertible single-qubit and twoqubit controlled-NOT (CNOT) quantum logic gates (all mathematically representable as unitary operators); this constitutes a universal set of quantum gates [8,9].
The QMS algorithm is a procedure which enables to perform a quantum MC sampling of thermal averages, using a sequence of invertible quantum logic gates on a register of qubits and classical measures on appropriate subportions of it. More specifically, QMS makes use of QPE in order to record the eigenvalues of a given quantum Hamiltonian without explicitly knowing the detailed structure of its eigenvectors. Before entering the overview of QMS, we discuss the working principle of QPE, lying at the heart of it.

A. Quantum phase estimation
The QPE algorithm applies to the following general problem. Suppose a given unitary operator U has an eigenvector |u with eigenvalue e iθ (0 ≤ θ ≤ 2π). Also assume that one is able to prepare the state |u and to perform controlled-U 2 j gates (j ∈ N). The latter operation needs a set of two quantum registers, the first one (1) containing the n qubits necessary to store |u ; the second one (2) containing an ancillary qubit, such that The purpose of QPE is to obtain the best l-bit estimate of θ, for a given integer value of l [34,35]. To achieve this goal, one discretizes the possible values of θ choosing a suitable encoding in binary notation with r bits [up to O(2 −(r+1) ) accuracy]. A sequence of the above controlled-U 2 j gates needs to be applied (varying j = 0, 1, 2, . . .) on r different ancillary qubits, here referred to as the control register. Since the quantum register that stores |u is prepared in an eigenstate of the operators U, U 2 , U 4 , · · · , the state of this register never changes, while the phase factors e iθ , e 2iθ , e 4iθ , . . . are propagated backwards in the control register. Finally, a quantum Fourier transform on the control register can be shown to provide an approximate binary encoding of θ, where each digit is represented by the state of a specific qubit state of the register. In particular, the algorithm outputs the best l-bit estimate of θ with probability > 1−ǫ, if the control register contains r = l+O(log(1/ǫ)) qubits and the obtained result is rounded off to its most significant l bits [34]. For further details on the QPE algorithm, see, e.g., Ref. [9].
In the specific implementation of QMS, the QPE procedure is applied to the unitary operator U = e −iHt , which, for a large class of physically relevant Hamiltonians, can be efficiently simulated on a quantum computer [10], after an appropriate discretization of the real-time evolution (Trotterization) [31,32]. Since an eigenvector |ψ α of H with eigenvalue E α is also an eigenvector of U (t) with eigenvalue e −iEαt , it is thus possible to encode, in the control register, the eigenvalues of H without explicitly knowing the detailed structure of its eigenvectors [35].

B. Implementation of Quantum Metropolis Sampling
We consider a generic quantum system that can be represented (or approximated) by a finite number n of qubits and is described by the Hamiltonian H. In principle, the eigenvalues of H could take any real value, but for the sake of simplicity let us assume that they can be exactly represented by a certain number r of qubits. This is always possible for commensurate energy spacings, by discretizing their possible values and using the binary representation of the integers 0 through 2 r − 1. As emerging from the discussion above in Sec. III A, and will be clarified later in Sec. V B, for systems with energy levels that are not multiple of a certain energy unit this introduces a systematic error, which however can be made arbitrarily small by increasing r (at the cost of increasing the complexity of the algorithm).
The idea of the QMS algorithm is to adapt to the QC framework the classical Markov Chain MC method [36], by generating a sequence of eigenstates of H distributed according to the Boltzmann statistical weight. While some steps of the original classical algorithm can be extended to the quantum case straightforwardly, the accept-reject step requires more care since, due to the no-cloning theorem [37,38], it is not possible to have a backup copy of the original configuration to be used if the proposed update is rejected.
To overcome this difficulty, the QMS algorithm uses four quantum registers: a register representing the system wavefunction |ψ with n qubits, two r-qubits registers holding information about the energy before (|E old ) and after (|E new ) a proposed update, and a single-qubit register with information about the acceptance (|acc ). All these registers can be arranged as in the following, when omitted, the subscripts corresponding to the various registers are intended as in (7). The QMS algorithm then proceeds as follows: Step 0: At the beginning, the input state (7) has to be initialized to zero for all registers, except the one holding the system state, which must be set to an eigenstate |ψ k of the Hamiltonian: |0, 0, 0, ψ k . Any eigenstate will work, so the simulation can be initialized, e.g., by making the wave-function collapse by means of an energy measure.
Step 1: A QPE procedure Φ (old) is performed between register 1 and register 2: E k being the energy of the eigenstate |ψ k . Thus, the global output state after this step is |0, 0, E k , ψ k .
Step 2: The choice of the trial update is implemented by randomly selecting a unitary operator C from a certain set C ⊆ U(2 n ), and applying it to the first register 3 , where x (C) k,p are the matrix elements of C. This operation is followed by a QPE Φ (new) between register 1 and register 3: In summary, this brings the output state of step 1 into a new (formal) superposition of eigenstates Step 3: This step is the quantum analogue of the accept-reject procedure: a one-qubit operator W (E p , E k ) is applied to the acceptance register 4, conditioned (by controlled gates) on both the old and new energy registers (2 and 3), such that Now, a classical measure is done on the acceptance register, which can only take the values 1 or 0. If the output of the measure is 1, one can proceed by performing a classical measure of energy on register 3, making the system qubits collapse onto a given energy eigenstate |1, E p , E k , ψ p . From this state, one resets 4 all the non-system qubits (registers 2, 3 and 4) to zero and then starts again the algorithm from step 1. Things are slightly more complex if the readout of the acceptance register is 0: in this case the trial state has to be rejected, and therefore the update step needs to be reverted, using the procedure described in step 4.
Step 4: First apply the inverse of the operator which represents the sequence of unitary operators applied in steps 2 and 3. Notice that register 3 will be set to zero by U † , while register 2 is kept untouched from step 1, giving us a reference value for the initial energy. A QPE Φ (new) followed by a classical measure on register 3 can return two possible results: if E new = E k , one has succeeded in reverting back the system state into an energy eigenstate with energy 5 E k ; it is thus possible to reset all non-system qubits and start again from step 1.
If instead E new = E k , one has to apply [Φ (new) ] −1 , followed by U and by a measure in the acceptance register (in order to reach a configuration analogous to the one at the end of step 3), and repeat step 4 until E new is equal to E k . From a practical point of view it is reasonable to set a maximum number of iterations for this procedure, after which the whole update is aborted and one has to restart from step 0. We finally emphasize that the possibility of using quantum gates controlled by classical measure is a key ingredient for the implementation of the QMS algorithm. This feature is used at the end of step 3 when, depending on the result of the classical measure on the acceptance register, one accepts (go to step 1) or rejects (go to step 4) the proposed configuration. Moreover, in case of reject, another classical control is needed to check whether E new is equal to E old or not.

C. Measuring observables
Using the QMS algorithm, we can generate a sequence of energy eigenstates whose energy is distributed according to the Boltzmann weight and, in particular, we can easily compute thermal averages of operators commuting with the Hamiltonian H. One is also interested, however, in computing thermal averages of a generic observable O, which does not commute with H. The measure of such an observable on the system register would however spoil the stochastic exactness of the QMS algorithm, since the state register would collapse away from the given energy eigenstate.
This problem is obviously absent in the classical algorithm; a possible way out is to measure the observable O every n th steps of the QMS algorithm, with n th large enough to let the system rethermalize between consecutive measures. A different strategy could be to start from step 0 of the QMS algorithm, perform n th updates to reach thermalization, measure the observable O and finally abort the simulation, repeating this procedure several times. In this way, we have however less or no variability in the initial configuration at step 0 (depending on the initialization procedure adopted), and by performing n th rethermalization updates after every measure we expect to have a more efficient sampling of the energies.
Note that, in both cases, (i.e. using rethermalization or different independent simulations) we have to fix the number n th , and this introduces a systematic error in the algorithm. By selecting n th large enough, this systematic can be made negligible with respect to the statistical errors that are present in every MC simulation, however this means that we cannot improve arbitrarily the accuracy by simply increasing the statistics. A detailed discussion of this point will be presented in Sec. V. Finally, let us explicitly remark that, if we are interested in computing the quantum thermal average of several operators, rethermalization updates are required after the measure of each of these observables, the only exception being the case of observables commuting with each other, which can be measured (in any order) with no need for in-between rethermalization.

IV. NUMERICAL RESULTS
In this section we present the results obtained by applying the QMS algorithm to the frustrated triangle system described by the Hamiltonian in Eq. (1).
States of the frustrated triangle can be represented using n = 3 qubits, while the energy of the system can be exactly encoded by using only one qubit (r = 1). Indeed, since the energy takes only two values (−J and 3J), by shifting and rescaling H we can define the new Hamilto-nianH = (1 + H/J)/4, whose spectrum is only made of eigenvalues equal to 0 or 1. The total number of required qubits is therefore n + 2r + 1 = 6.
To proceed with the QMS algorithm we need two basic ingredients: to be able to perform a QPE and to determine the set C of unitary operators to be used for generating trial states. As previously discussed, the set C has to ensure ergodicity (in this context, by ergodicity we do not mean the possibility to actually explore the whole Hilbert space, but only the whole set of energy eigenstates) and if U ∈ C, also U † has to belong to the same set (reversibility).
The naive choice of single spin flips 1⊗1⊗σ x , 1⊗σ x ⊗1 and σ x ⊗1⊗1 does not ensure ergodicity, as evident from the observation that all these operators commute with H. An ergodic choice can be obtained instead by replacing the spin flip operator σ x in the previous expressions by the Hadamard gate 1 √ 2 (σ z + σ x ) [8]. Since these transformations are involutions, reversibility is also automatically satisfied. With this choice, we obtain an average acceptance probability in the Metropolis step that goes from 99% for β = 0.1 to 91% for β = 1.0.
For what concerns the QPE, this actually boils down to a real-time evolution (see Sec. III A) and, for the frustrated triangle, time evolution is particularly simple: since all terms in H commute pairwise and their square is proportional to the identity, we can easily construct the time evolution operator as a product of terms like We will come back later to the systematics associated to the QPE when such an analytical expression is not available. Now we have all the ingredients required for the application of the QMS algorithm and, in order to test its effectiveness, we start by studying the dependence on the temperature of the internal energy H , comparing numerical data with exact analytical results. Such a comparison is performed in Fig. 1 and we can see that there is excellent agreement between the two determinations.
We now move to the computation of thermal averages of operators not commuting with the Hamiltonian. We have tested several Hermitian operators, but here we only report the results obtained for The outcomes obtained in all other cases are equivalent. After some preliminary test (see Sec. V A for a thorough analysis) we decided to use 200 rethermalization steps after each measure of A, and in Fig. 1 we show the results obtained with a statistics of about 10 7 measures. Even in this case, the agreement between the numerical estimates and the exact results is remarkable; this is a direct proof of the fact that the algorithm properly samples the configurations and effectively solves the sign problem of this model. After having checked that the QMS algorithm correctly samples the states of the system, we can investigate some aspects related to the efficiency of the method. In particular, a possible weak point of the QMS algorithm is the iterative procedure required in the step 4 (see Sec. III B) to project back a rejected update: if a very large number of iterations is required, the algorithm obviously loses its effectiveness. However we found that, at least for the frustrated triangle studied in this work, the number of reverting steps required for each iteration is relatively small, and its distribution can be well approximated by an exponential, as shown in the histogram of Fig. 2 for some values of β. The distribution of the number of reverting steps gets broader for small β values, signaling that the algorithm is less efficient in this regime; on the other hand, this is mitigated by the fact that at low β the acceptance is higher, so that reverting steps occur less frequently. A possible way of reducing the impact of the reverting steps on the simulation time would be to look for a set C of moves with larger acceptance probability, so to reduce the occurrence of the iterative step.
One last comment must be done on the number of quantum gates required to perform a single Metropolis step. In order to perform the QPE and apply the operator U , i.e. before measuring the acceptance register, we need ∼ 100 gates. If we accept the new configuration we need additional ∼200 gates to go back to step 1, otherwise, in case of rejection, we need ∼100 gates in order to perform a single reverse attempt. In the counting of the quantum gates all the classical measures and all the quantum gates controlled by a classical measure have been excluded. We see that, even for a very simple system like the frustrated triangle, the number of quantum gates used is large, therefore, in order to run the algorithm on a real quantum computer, a very high fidelity on the quantum gates is required.

V. SYSTEMATICS OF THE QMS ALGORITHM
The QMS algorithm generally presents a certain number of systematic errors. In the previous section, we already discussed the rethermalization updates to be performed after each measure of an observable not commuting with the Hamiltonian. For the frustrated triangle this was the only source of systematic error in the QMS algorithm. Here we discuss the systematics that have to be taken into account, beyond rethermalization, when using QMS to simulate a generic quantum system.
Digitalization of the energy: in the QMS update, energy has to be stored in two quantum registers, each composed of r qubits, and if the spectrum of the Hamiltonian cannot be exactly mapped to an r−bit binary number, one introduces a truncation error in the energies and thus a systematic in the sampling. Even for toy models, this is clearly unavoidable whenever energy eigenvalues are incommensurable with each other, like for a three-state system with energy levels 0, 1/ √ 2, and 1. Digitalization of the state: this problem is present for systems with continuous degrees of freedom (like gauge theories associated to continuous Lie groups). In such cases we need to choose a proper digitalization of the state, i.e. a way of approximating the state by using only n qubits. For more details, see Refs. [39][40][41][42].
QPE implementation: as briefly recalled in Sec. III A, the implementation of the QPE algorithm requires a realtime unitary evolution on the system register using e −iHt . In practice, such an evolution cannot be performed exactly, but has to be approximated by a Trotterization, thus introducing a systematic error related to the finite time-step (see, e.g., Refs. [43][44][45][46][47][48] for recent discussions about this systematic effect in similar contexts).
The effect of all these systematics can be reduced at the expense of increasing the computational complexity of the algorithm: by increasing the number of used qubits, the digitalization of both energies and states can be made arbitrarily accurate and similarly the Trotter time-step can be reduced at will, making the QPE more and more accurate. While this is true in theory, in practice one is interested in understanding how these systematics scale with the computational complexity, in order to understand whether the algorithm is useful or not for the problem at hand.
The last two systematics of the list (i.e. digitalization of the state and QPE implementation) are not specific of the QMS, but are common to all QC algorithms. For this reason, they have already been discussed in the literature and we will not further investigate them, thus concentrating on the rethermalization and the digitalization of energy.

A. Rethermalization
We have performed several simulations with different numbers of rethermalization steps, to quantify the number of needed rethermalization updates. As an example, the results obtained for the average values H and A [defined in Eq. (16)] are shown in Fig. 3 for β = 1. The observed behavior is analogous to that seen in the thermalization of classical MC simulations: average values depend on the number of rethermalization steps n th adopted, but they converge to an asymptotic value when n th increases.
This asymptotic value is the correct thermal average, which can however be reached at different values of n th for different observables, as clearly seen in Fig. 3; this is just a consequence of the fact that different observables have different thermalization times. In order not to introduce systematics in the final estimates, n th has to be chosen large enough that all the observables have reached their plateau values with an accuracy that is larger than the target precision of the simulation.

B. Digitalization of the energy
In order to assess how much the QMS algorithm relies on an accurate estimation of the energies, we considered a simple system made up of two qubits representing an Hamiltonian H with the following four energy levels: With this setup, energy levels can never be exactly digitalized, and we can also investigate the effect of having two nearby levels (1/ √ 2 and 3/4) that cannot be resolved if the number of used qubits is too small.
As a preliminary test, we isolated the QPE step: in Fig. 4 we show the results of a single application of the QPE procedure Φ to the eigenstate of H associated to the eigenvalue λ 1 = 1 √ 2 , using a variable number r of qubits for the energy register. Since λ 1 cannot be represented by a fraction of powers of 2, a measure in the energy register does not always collapse to the same state, but has a certain probability to collapse to a state with neighboring energy, so that each measure is affected by a sampling error, which however vanishes exponentially in r.
More interesting are the effects of errors in the estimates occurring during the full QMS algorithm, as can be observed from the average of the energies measured, shown in Fig. 5 for three different digitalizations of the energy (r = 4, 6, and 8). Using 4 qubits seems not to give a sufficiently accurate estimate along the explored range in β, while 6 and 8 qubits show reasonable accuracy, even if the convergence behavior is not always clear, as apparent, for example, at β = 1. Nevertheless, comparing the energy probability distribution of the samples with the exact one (Fig. 6), we can observe a clearer convergence to the exact weights at increasing number of qubits r. Notice that a proper discrimination between neighboring eigenvalues in the spectra (as 1 if the number of qubits is sufficient.

VI. CONCLUSIONS
In this paper we have described, implemented, and tested a quantum algorithm able to measure thermal averages for quantum mechanical systems. Our results are for a very simple system of frustrated spins, that, if rewritten in terms of Euclidean path integrals exhibits the infamous sign problem, that makes the use of clas- sical Monte Carlo simulations impossible. In this sense, our toy model shares the same problems of more complex and physically meaningful problems, such as the study of the phase structure of Lattice QCD at finite temperature and density. It therefore conceptually paves the way to quantum simulations of those systems, as soon as appropriately powerful quantum computers become available in the future.
The algorithm chosen for our study is the Quantum Metropolis Algorithm introduced in Ref. [16], which has been briefly reviewed in Sec. III. Even if inspired by the standard Metropolis algorithm, it cannot be considered a Markov chain in the classical sense, since the measurement of observables interferes with the Markov process itself, thus destroying equilibration and leading to the need for rethermalization after each measurement is taken (unless the observable commutes with the Hamiltonian). Of course this is not a peculiar feature of the algorithm, but rather something which one expects in general, because of the quantum nature of the system and of the algorithm.
We have verified that our quantum algorithm reproduces, within statistical errors, exact results for the thermal averages of the energy and of an additional operator that does not commute with the Hamiltonian. We have quantified the size of the statistical and sampling errors and their behavior with respect to the number of simulation steps. We have also discussed and measured the systematic errors due to several different effects, such as the accuracy of the phase estimation process, the sampling accuracy as a function of the qubit register-size, as well as some conceptually unavoidable performance losses associated to the reverting steps of the algorithm.
All these systematic errors, together with the computing depth of the complete algorithm, will ultimately play a key role to define the actual viability of this quantum algorithm as one moves from simple toy models to the study of real-life complex physical systems. This perspective suggests several avenues for further studies, ranging from the implementation of the present model onto actual QC machines, to the implementation of the algorithm for more complex (possibly gauge invariant) models, to the analysis of the optimal balance point between computational accuracy and efficiency versus register-size and algorithm complexity: work is in progress along these directions. Here we present a derivation of Eq (4), which was used in Sec. II to better understand the structure of the pathintegral representation of the frustrated triangle. We start from the Taylor expansion of the Euclidean evolution operator: It is then easy to check that the Hamiltonian (1) satisfies the property H 2 = (3J 2 1 + 2JH), therefore we have, for any integer n, a relation of the form H n = α n 1 + β n H .
To fix the values of α n and β n we can diagonalize both members of this equation and remember that the eigenvaules of H are just −J and 3J. In this way, we immediately get