Quantum Computation of Finite-Temperature Static and Dynamical Properties of Spin Systems Using Quantum Imaginary Time Evolution

Developing scalable quantum algorithms to study finite-temperature physics of quantum many-body systems has attracted considerable interest due to recent advancements in quantum hardware. However, such algorithms in their present form require resources that exceed the capabilities of current quantum computers except for a limited range of system sizes and observables. Here, we report calculations of finite-temperature properties including energies, static and dynamical correlation functions, and excitation spectra of spin Hamiltonians with up to four sites on five-qubit IBM Quantum devices. These calculations are performed using the quantum imaginary time evolution (QITE) algorithm and made possible by several algorithmic improvements, including a method to exploit symmetries that reduces the quantum resources required by QITE, circuit optimization procedures to reduce circuit depth, and error mitigation techniques to improve the quality of raw hardware data. Our work demonstrates that the ansatz-independent QITE algorithm is capable of computing diverse finite-temperature observables on near-term quantum devices.


I. INTRODUCTION
Quantum computers have long been considered as a potential tool to simulate quantum many-body systems [1][2][3].While near-term quantum devices have made rapid progress in simulating ground-state properties [4][5][6][7][8][9][10] and dynamics [11][12][13][14][15], the study of finite-temperature physics on quantum computers is less understood and established [16].Early works on digital quantum simulation of finitetemperature physical systems involved thermalizing the quantum simulator by coupling to a bath comprised of ancilla qubits [17][18][19] or sampling thermal states using quantum versions of the Metropolis algorithm [20,21].These schemes require prohibitively large numbers of qubits and deep circuits and are hence out of reach for near-term quantum devices.
More practical variational algorithms have been proposed in recent years, such as protocols to construct thermofield double states [22,23] and machine learning-based methods to construct Gibbs states [24][25][26][27][28].However, the accuracies of these variational schemes depend on the quality of the ansatz.While other non-variational alternatives exist, they are subject to various assumptions.For example, the minimal effective Gibbs ansatz [29] algorithm generates a minimal ensemble of pure states but presumes the eigenstate thermalization hypothesis.
Recently, the quantum imaginary time evolution (QITE) algorithm was introduced [30].Compared to variational-based algorithms of imaginary time evolution on quantum computers [31][32][33], QITE is ansatz-independent.The QITE algorithm approximates imaginary time evolution with unitary operators over a domain of qubits and is able to reach the ground states of systems with a few sites.QITE can also be used to calculate finite-temperature quantities, for instance by combining with sampling techniques such as the minimal entangled typical thermal states (METTS) algorithm [34,35], together denoted as the quantum METTS (QMETTS) algorithm.However, the original work on QITE [30] focused on the general formalism, while reduction and optimization of quantum resources were not thoroughly investigated.Subsequent development of QITE [36][37][38][39] proposed several variations of the original algorithm, but the practical evaluation of finite-temperature properties on existing quantum devices remains largely unaddressed.
Here, we report QITE-based calculations of finitetemperature static and dynamical properties of onedimensional spin systems with up to four sites on fivequbit IBM Quantum devices.The computed observables include finite-temperature energies, static and dynamical correlation functions, and excitation spectra.These calculations are made possible by several algorithmic improvements.First, we exploit symmetries in the spin Hamiltonians to reduce Pauli strings in the QITE unitaries, thus reducing the overall required quantum resources.Second, circuit optimization procedures including gate decomposition and circuit recompilation are used to further reduce circuit depth.Third, error mitigation techniques, namely post-selection, readout error mitigation and phase-and-scale correction, are used to improve the quality of raw hardware data.Our work demonstrates that with efficient use of quantum resources and effective error mitigation strategies, the ansatz-independent QITE algorithm is capable of com-puting diverse finite-temperature observables on nearterm quantum devices.
This paper is organized as follows.In Sec.II we review the QITE algorithm and propose a quantum circuit to evaluate finite-temperature dynamical correlation functions.In Sec.III we introduce the algorithmic improvements including Pauli string reduction, circuit optimization and error mitigation that enabled us to obtain accurate results from hardware.Section IV presents the results of our two-site and four-site calculations.Finally, we conclude and suggest directions for future studies in Sec.V. m=1 ĥ[2m] corresponds to the even-odd Trotterization used in one-dimensional tensor network calculations of quantum many-body systems [40].We consider first-order Trotterization [41] of the full imaginary time propagator e −β Ĥ : where n β is the number of imaginary time steps and ∆τ = β/n β .
The QITE algorithm approximates each imaginary time propagator e −∆τ Ĥ[l] by a unitary operator where x[l] µ are real coefficients and σ µ are Pauli strings.
Here we use the notation σ 0 = I, σ x = X, σ y = Y, σ z = Z to denote the identity and the Pauli matrices, so that each Pauli string can be written in the form σ µ = N −1 j=0 σ µj where σ µj acts on qubit j and µ j ∈ {0, x, y, z}.The Pauli strings σ µ are chosen from the set where P ĥ[m] is the set all Pauli strings over a domain of D qubits larger than or equal to the support of ĥ[m].To apply the QITE unitaries, without an efficient decomposition scheme each unitary needs to be further Trotterized as The coefficient vector x[l] is found by minimizing the square of the difference between the unitarily evolved state e −i∆τ Ĝ[l] |Ψ and the imaginary-time-evolved state 2 .This minimization results in a linear system where In our implementation, we expand the exponential e −∆τ Ĥ [l] in b[l] and c[l] to second order in ∆τ : To construct the linear systems, given the terms in Eqs. ( 6), ( 8) and ( 9) we measure operators of the form The QITE algorithm is carried out by iterating the procedure of constructing the circuit from the QITE unitaries obtained at the previous imaginary time steps, measuring the operators in Eq. (10), constructing the linear system in Eq. ( 5), solving for x[l], and propagating the state with the new unitary e −i∆τ Ĝ [l] .

B. Finite-temperature dynamical correlation functions
Finite-temperature static observables have been previously computed on quantum hardware using the QMETTS algorithm by averaging over the observable evaluated from each METTS sample state [30].In this work, we show that finite-temperature dynamical observables, in particular finite-temperature dynamical correlation functions, can be computed using a similar averaging procedure as for finite-temperature static observables.
On quantum computers, dynamical correlation functions can be calculated using the circuit reported in Refs.[42,43]; this circuit has been recently used to compute neutron scattering cross-section [14] and magnon spectra [15] on quantum hardware.To obtain finitetemperature dynamical correlation functions, we insert the QITE circuit into the dynamical correlation function circuit, resulting in the circuit shown in Fig. 1.The ancilla qubit is initialized in |0 and the system qubits are initialized in |Ψ .Define |Ψ(τ ) = e −τ Ĥ |Ψ /||e −τ Ĥ |Ψ || as the state initialized in |Ψ and evolved to imaginary time τ , and |Φ(τ ) as the QITE-evolved state that approximates |Ψ(τ ) .Let subscript a (s) denote quantities on the ancilla qubit (system qubits).To evaluate finitetemperature observables at an inverse temperature β, we evolve the initial state by QITE to β/2 so that the joint ancilla-system density operator prior to measurement is where and Ũ ( Ṽ ) is the controlled-U (controlled-V ) gate.
Measuring X (Y ) on the ancilla yields the real (imaginary) part of the dynamical correlation function on a single QITE-evolved basis state: If the initial states are the METTS sample states, an unweighted average over the initial states yields the finite-temperature dynamical correlation function U (t)V β .In this work, we consider trace evaluation in the exact expression of an observable Ô at finite temperature: The numerator trace and the denominator trace are either evaluated by full sampling over the entire Hilbert space, denoted as full trace evaluation, or by random sampling over a subspace of the Hilbert space, denoted as stochastic trace evaluation.If Ô = U (t)V , Eq. ( 16) yields the finite-temperature dynamical correlation function U (t)V β .Ô can also be a static observable, in which case Eq. ( 16) yields the static observable at finite temperature.

A. Pauli string reduction by Z2 symmetries
If we include all 4 D Pauli strings over each domain consisting of D qubits, each QITE unitary e −i∆τ Ĝ[l] applied as in Eq. (4) yields O(N 4 D ) multi-qubit rotation gates by the standard rotation gate decomposition [44], which results in a circuit too deep on near-term quantum devices even for D = 2.Because of this prohibitive resource overhead, we describe a systematic method to reduce the number of Pauli strings in the QITE unitaries when the Hamiltonian and initial state have Z 2 symmetries.
Z 2 symmetries on qubit Hamiltonians have direct parallels with the stabilizer formalism in quantum errorcorrecting codes [45].Suppose the Hamiltonian has d Z 2 symmetries, i.e.Ĥ commutes with elements of a group isomorphic to Z d 2 generated independently by d Pauli strings, and the initial state is in the +1 eigenspace of all d generators.If we regard the symmetry group Z d 2 as the stabilizer S, the symmetry sector of the initial state corresponds to the stabilizer subspace V S .
In stabilizer codes, the normalizer of the stabilizer N (S) includes all Pauli strings that commute with elements of the stabilizer S, and all valid operations on the code space are in the quotient group N (S)/S.Intuitively, to preserve Z 2 symmetries, among all Pauli strings from P Ĥ[l] the QITE unitaries should only include those from the quotient group N (S)/S.We now show that the original QITE algorithm subsumes the requirement that the Pauli strings should be chosen from N (S)/S because the action of the unitary e −i∆τ Ĝ [l] with the Pauli strings from the unreduced set P Ĥ[l] is the same as the action with Pauli strings from the reduced set P Ĥ[l] ∩ N (S)/S.This result is stated as the following proposition, proved and discussed in Appendix A.
Proposition.Suppose QITE is applied to approximate the imaginary time propagator e −∆τ Ĥ[l] on the state |Ψ .If there exists a stabilizer S such that every element of S commutes with Ĥ[l] and Further reduction in the number of Pauli strings can be achieved by recalling from Ref. [30] that when the Hamiltonian and the initial state are real in the computational basis, the state after imaginary time evolution must be real.Thus, only Pauli strings with an odd number ofw Y need to be included in the QITE unitaries.Since Z 2 symmetries and the conditions of a real Hamiltonian and initial state are independent, when both conditions are satisfied, the number of Pauli strings can be reduced using both conditions, in which case the reduced set of Pauli strings is modified to P Ĥ[l] ∩ N (S)/S ∩ {σ µ : j δ µj ,y ≡ 1 (mod 2)}.In both panels the energy trajectories using reduced numbers of Pauli strings match the energy trajectories without reduction, which also match the energy trajectories from exact imaginary time evolution.
In practice the Proposition is used inductively on the Trotter terms Ĥ[l], which implies the stabilizer need to be chosen such that every element of the stabilizer commutes with all Ĥ[l], or equivalently with Ĥ.For spin Hamiltonians, the stabilizer generators are usually global Z 2 symmetries such as Z ⊗N and X ⊗N .For general Hamiltonians, the Z 2 symmetries can be found by Gaussian elimination on the parity check matrix formed from the Hamiltonian terms [46].
To confirm our Pauli string reduction scheme, we compare the QITE energy trajectory as a function of imaginary time simulated without noise or measurement sampling on a single initial state with and without reduction of the Pauli strings in the QITE unitaries by Z 2 symmetries.The Hamiltonians we study include the transversefield Ising model (TFIM) Hamiltonian and the Heisenberg XXZ Hamiltonian with open boundary conditions assumed for both.In Fig. 2 we plot energy versus imaginary time calculated with QITE on a single initial state.The unreduced set of Pauli strings only includes Pauli strings with odd numbers of Y because the Hamiltonian and initial state are real in the computational basis.We choose a sufficiently small imaginary time step size ∆τ to ensure that the Trotter errors from expansion in Eq. ( 4) are negligible.Figure 2(a) plots the energy trajectory for the initial state |0001 in the four-site TFIM with J = h = 1.The Hamiltonian and the initial state have a Z 2 symmetry Z 0 Z 1 Z 2 Z 3 .By combining reduced Pauli strings from all three D = 2 domains, we obtain six Pauli strings in the QITE unitaries The 120 Pauli strings in the QITE unitaries without reduction is reduced to the 6 Pauli strings In both panels of Fig. 2, the energy trajectories using reduced numbers of Pauli strings match the energy trajectories without reduction, which also match the energy trajectories from exact imaginary time evolution.

B. Circuit optimization
Even with reduction of Pauli strings in the QITE unitaries by Z 2 symmetries, applying the QITE unitaries as in Eq. ( 4) may still result in a circuit too deep to be implemented on current quantum hardware.In this section we describe circuit optimization techniques that further reduce circuit depth.
In two-site calculations, both the QITE circuit and the real time evolution circuit can be optimized to constant depth with a standard one-and two-qubit gate set, regardless of the number of imaginary and real time steps.For example, in two-site TFIM there is only one Pauli string X 0 Y 1 in the QITE unitaries after reduction by the Z 2 symmetry Z 0 Z 1 .Suppose the unitary applied to the state at the k-th imaginary time step is e −i∆τ x k X0Y1 .Then the unitaries at all imaginary time steps can be multiplied into a single two-qubit rotation gate e −iθX0Y1/2 where θ = 2∆τ k x k .For real time evolution, the two-qubit operator e −i Ĥt is decomposed by the KAK decomposition [47][48][49][50] into six single-qubit gates and two CNOT gates.
In four-site calculations, neither the QITE circuit nor the real time evolution circuit is of constant depth.If we Trotterize the QITE unitaries as in Eq. ( 4) and similarly for the real time propagator, the circuit is too deep to be accurately implemented on existing quantum devices.Therefore, we recompile the circuit by fitting each QITE unitary e −i∆τ Ĝ [l] or the real time propagator e −i Ĥt to a parametrized circuit [51][52][53] using an open-source tensor network-based quantum simulation library [54].In Fig. 3, we show the recompiled four-site QITE circuit, where the U 3 gate is a generic single-qubit gate defined as The four U 3 gates at the left constitute the base gate round.Each additional gate round consists of a layer of CNOT gates and a layer of single-qubit gates.The additional gate rounds alternate between even-odd and odd-even pairs of qubits.Let the target unitary be U targ and the recompiled unitary be U rec (θ), where θ is a composite variable denoting all the angles.Given a reduced density operator ρ on the finite domain acted on by the target unitary, write ρ targ = U targ ρ U † targ and ρ rec (θ) = U rec (θ)ρ U † rec (θ).The optimal recompiled unitary is found by performing a gradient descent to maximize the fidelity [55] Since the QITE unitaries are real, we use the oneparameter single-qubit gate R y (θ) = U 3 (θ, 0, 0) in the recompiled circuit for QITE, while for real time evolution we keep the U 3 gate as the parametrized single-qubit gate.

C. Error mitigation
To mitigate the effect of hardware noise on the measurement results, we post-process our hardware data by error mitigation methods including post-selection, readout error mitigation and phase-and-scale correction.Post-selection and readout error mitigation are applied to the measurement outcomes at each imaginary time step in the QITE subroutine; phase-and-scale correction is applied to the final computed finite-temperature dynamical correlation function as a single-step post-processing.
Post-selection is performed on Z 2 symmetries discussed in Sec.III A. When the Hamiltonian and the initial state have Z 2 symmetries, the final state after imaginary or real time evolution should have the same stabilizer parities as the initial state.However, during execution of the circuit, gate errors and qubit decoherence can induce nonzero overlap of the qubit state with the subspace of the wrong parity.Post-selection can mitigate these undesirable effects by discarding measurement outcomes with the wrong parity [56,57].
We specifically consider the symmetry from a single stabilizer generator.If the operator to be measured is an ancilla operator, we can simply measure the stabilizer generator on all the system qubits and read off the parity without interfering with measurement of the ancilla.If the operator to be measured acts on system qubits, we need to simultaneously measure the operator and the stabilizer generator, which is possible because all operators in Eq. ( 10) commute with the stabilizer generator by our choice of Pauli strings in the QITE unitaries in Sec.III A. Specifically, each operator and the stabilizer generator can be simultaneously measured by using Clifford gates to transform the Pauli string components of the operator and the stabilizer generator until they are qubit-wise commuting, so that their expectation values can be read off on different qubits [58][59][60][61].
Figure 4 shows the circuit to simultaneously measure the stabilizer generator Z 0 Z 1 Z 2 Z 3 and a Pauli string that commutes with it in a four-site QITE calculation.The sequence of CNOT gates after the QITE circuit in Fig. 4 transforms Z 0 Z 1 Z 2 Z 3 to Z 3 .Since Z 3 acts on a single qubit, it necessarily qubit-wise commutes with the transformed Pauli string.In practice, we stop applying CNOT gates when the transformed Pauli string becomes qubitwise commuting with the transformed stabilizer generator, and therefore the number of CNOT gates added to the end of the circuit ranges from zero to three.
To account for noise in the final measurement, the built-in readout error mitigation routine in Qiskit [62] is applied to each measurement outcome.Because of the small size of the systems we study, for N qubits we carry out full calibration on all 2 N initial states.Application of the inverse of the calibration matrix to the raw measurement counts is performed by the default least-square fitting method.
We assess the effectiveness of applying post-selection and readout error mitigation at every imaginary step of QITE by simulating the finite-temperature energies of two-site and four-site TFIMs using full trace evaluation with measurement sampling and the noise model from ibmq rome.In both panels of Fig. 5, QITE is applied with Pauli strings reduced and circuits optimized.In particular, four-site QITE unitaries are of domain size D = 2 and recompiled with three rounds of gates.From Fig. 5 we can see that both readout error mitigation and post-selection shift the raw data toward the exact data, confirming the effectiveness of both schemes in reducing the effect of noise.Furthermore, a combination of readout error mitigation and post-selection is observed to be most effective in mitigating the errors, which is not apparent on two sites in Fig. 5(a) presumably because of the small size of the system but clearly evident on four sites in Fig. 5(b).
In calculations of finite-temperature dynamical correlation functions, the ancilla qubit is in the state |+ before entangling with the system qubits.When there is a long sequence of gates in the real time propagator e −i Ĥt , decoherence of the ancilla qubit such as amplitude damping to the qubit ground state |0 and depolarization will significantly affect the X and Y measurement results on the ancilla.To mitigate the effect of ancilla decoherence, we apply phase-and-scale correction [14,15] as a single-step post-processing to the result at the end of the calculation.The only finite-temperature dynamical correlation function considered in this work is Z 0 (t)Z 0 β , which is equal to 1 analytically at t = 0. Hence, we apply phase-and-scale correction by dividing the raw hardware Z 0 (t)Z 0 β at each t by the raw hardware Z 0 (t = 0)Z 0 β to enforce the condition Z 0 (t = 0)Z 0 β = 1.

IV. RESULTS
Experiments of computing finite-temperature observables were conducted on IBM Quantum devices ibmq bogota and ibmq rome [63], both of which consist of five qubits arranged on a chain with nearest-neighbor interactions and similar error rates.IBM's open-source library Qiskit [64] was used to implement our algorithms on the devices.In each calculation, the N system qubits 0, ..., N − 1 are arranged adjacent to each other and the ancilla is closest to system qubit 0.
The systems we study are sufficiently small that we apply QITE to approximate the full imaginary time propagator e −∆τ Ĥ at each imaginary time step, which is equivalent to setting L = 1 and Ĥ[1] = Ĥ in Eq. ( 1).The QITE linear systems in Eq. ( 5) are solved by a conjugate gradient method.Because hardware noise and measurement sampling lead to ill-conditioned A matrices in the QITE linear systems, we add a regularizer of 0.2 to the diagonal elements of each A matrix in the four-site calculations.
Each calibration circuit used for readout error mitigation is repeated 1000 times; each Pauli string measurement circuit used to construct the QITE linear systems is repeated 8000 times.Error bars from full trace evaluation result only from measurement sampling and are the size of the markers in most figures; error bars from stochastic trace evaluation originate from both measurement sampling and initial state sampling.A detailed description of error bars in full and stochastic trace evaluation is given in Appendix B.

A. Two-site calculations
We study the two-site TFIM defined in Eq. ( 17) by setting h = 1 and varying J.The finite-temperature energy E β and static correlation function X 0 X 1 β are calculated on the Hamiltonians with J = ±1, ±3, while the finite-temperature dynamical correlation function Z 0 (t)Z 0 β and excitation spectra are calculated on the Hamiltonian with J = 3.In all calculations, finitetemperature observables are calculated by full trace evaluation and the circuits are optimized according to the circuit optimization procedures in Sec.III B.
Figure 6 shows the finite-temperature energy E β and static correlation function X 0 X 1 β of the twosite TFIM with J = ±1, ±3 from β = 0 to β = 2.In both Fig. 6(a) and Fig. 6(b) the finite-temperature observables obtained on hardware are in good agreement with exact values.Further, if we regard each finite-temperature variable as a function of J, analyti- cally it can be shown that E β (J) = E β (−J) and X 0 X 1 β (J) = − X 0 X 1 β (−J).This relation is satisfied in the hardware data.In each of the observables we calculated, the mean absolute percentage error between hardware and exact results averaged over all β ranges from 1% to 4%.
Next, finite-temperature dynamical properties were calculated on the two-site TFIM with J = 3, h = 1.The dynamical correlation function Z 0 (t)Z 0 β is evaluated from β = 0 to β = 2 and at real time from t = 0 to t = 8π.Figures 7(a) and 7(b) show the real and imaginary parts of Z 0 (t)Z 0 β at β = 0.2 and β = 1.8 up to t = 4π.From Figs. 7(a) and 7(b) we see that even without phase-and-scale correction, the real and imaginary parts of Z 0 (t)Z 0 β agree well with the exact results at both small and large β, presumably due to the constant and shallow depth of the real time evolution circuit.
The spectral density S(ω) is obtained by a discrete Fourier transform of the dynamical correlation function Z 0 (t)Z 0 β .Specifically, at each β where n t is the total number of points in the time series, t m = m∆t, and ω k = 2πk/n t ∆t.With this definition of Fourier transform, the peaks at positive (negative) frequencies correspond to emissions (absorptions) of excitations of the system.In Fig. 7(c) we plot the excitation spectra of the twosite TFIM at β = 0.2.The exact excitation spectrum is obtained by a Fourier transform of the exact Z 0 (t)Z 0 β at the same points in real time as the Z 0 (t)Z 0 β obtained on hardware.From the plot, we can see that the hardware excitation spectrum agrees well with the exact excitation spectrum.The frequencies ω = 0, ±5.94, ±7.18, at which the peaks in the excitation spectra are located, correspond to the excitation frequencies of ω = 0, ±6.00, ±7.21 from exact diagonalization of the Hamiltonian.The deviation of the frequencies in the excitation spectra obtained from hardware Z 0 (t)Z 0 β compared to the frequencies obtained from exact diagonalization is due to the finite real time domain in our hardware calculations.
To analyze the evolution of the excitation spectra across different temperatures, we plot the amplitudes at the two emission frequencies versus β in Fig. 7(d).Analytically, the amplitude of the transition from an initial state , where E f is the energy of the final state and Z is the partition function.In the two-site TFIM, the only allowed transitions are between the two states in each of the two-dimensional eigenspaces of Z 0 Z 1 with eigenvalues ±1.The frequency ±7.18 corresponds to a transition in the +1 eigenspace, where the ground state lies, and the frequency ±5.94 corresponds to a transition in the −1 eigenspace, where the first excited state lies.As the temperature decreases from infinite temperature (β increases from 0), the populations in the two lowest states first increase until the ground state population dominates over that of the first excited state at around β = 0.4, a trend reproduced by the amplitudes obtained from hardware data in Fig. 7(d).Thus, Fig. 7 shows that quantum hardware accurately captures the finite-temperature dynamics of the two-site TFIM across a wide range of temperatures.

B. Four-site calculations
We next proceed to four-site spin systems.We study the four-site TFIM defined in Eq. ( 17) with J = 3, h = 1.Full trace evaluation is employed unless otherwise specified.
First, let us consider the gate count in the four-site circuits.For the four-site TFIM with D = 2, after reduction by the Z 2 symmetry Z 0 Z 1 Z 2 Z 3 there are six weighttwo Pauli strings, which are the ones given in Eq. ( 19).If we Trotterize the QITE unitaries as in Eq. ( 4), each unitary requires 12 CNOT gates by the standard rotation gate decomposition [44], which becomes unfeasible on near-term quantum hardware after the first few imaginary time steps.When D = 4, even after reduction by one Z 2 symmetry there are still 28 Pauli strings in each QITE unitary.After Trotterization and rotation gate decomposition, each QITE unitary requires more than 50 CNOT gates for a single imaginary time step.Hence, when D = 2, we compute finite-temperature observables both by Trotterizing and by recompiling the QITE unitaries with three gate rounds; when D = 4, we only recompile the QITE unitaries with three gate rounds.To obtain dynamical correlation functions, we additionally recompile the real time propagator e −i Ĥt with five gate rounds.The number of gate rounds is chosen so that the fidelity in Eq. ( 22) is at least 0.999 on average in each calculation.
Figure 8 shows the finite-temperature energy E β and static correlation functions X 0 X 1 β , X 0 X 2 β , X 0 X 3 β of the four-site TFIM.From the figure, we can see that the finite-temperature observables calculated with Trotterized D = 2 QITE unitaries deviate from those calculated with D = 2 recompiled QITE unitaries or D = 4 recompiled QITE unitaries after β = 0.1.This deviation is due to the deep circuit resulting from 12 layers of CNOT gates per imaginary time step, compared to 3 layers of CNOT gates per imaginary time step in the recompiled circuit.Moreover, even the recompiled QITE unitaries are not able to track the exact finite-temperature observables after the first few β.In particular, the slope is reversed compared to the exact result after β = 0.5.QITE up to β = 0.5 corresponds to 5 imaginary time steps and hence 15 layers of CNOT gates, which is almost at the limit of circuit depth on these quantum devices.
We examine more closely the calculations with recompiled QITE unitaries and focus on the data at β ≤ 0.5.For both D = 2 and D = 4, the recompiled QITE unitaries apply the same number of layers of gates at each imaginary time step, so the difference in calculated observables is formally caused by the difference in domain sizes.Since the domain size should grow with correlation length and hence with imaginary time [30], we should expect the hardware results to be closer to the exact results with D = 4 than with D = 2 especially at large β.However, this hypothesis only holds for X 0 X 2 β and X 0 X 3 β but not for E β and X 0 X 1 β .Failure of this hypothesis is likely due to that fact that the larger 28dimensional linear system arising from D = 4 incorporates more errors from hardware noise compared to the smaller 6-dimensional linear system arising from D = 2.
To explore the scalability of our method, we compare stochastic trace evaluation with full trace evaluation in calculating the finite-temperature energy of the four-site TFIM.Stochastic trace evaluation is performed by uniformly selecting initial states in the full trace evaluation result with recompiled D = 2 QITE unitaries.In Fig. 9, we plot the stochastic trace evaluation results with 10 and 20 samples along with the full trace evaluation and exact results; the inset shows the running average of E β versus number of samples n samples .As can be seen from the figure, random sampling with 10 samples already reproduced the results from full sampling on all 16 initial states, indicating that using scalable sampling schemes is a promising approach to studying larger systems.Finally, in Fig. 10 we show the dynamical properties of the four-site TFIM with J = 3, h = 1 at β = 0.2.The calculation is implemented by recompiling D = 2 QITE unitaries with three gate rounds and real time propagation with five gate rounds.Figure 10(a) shows the real and imaginary parts of Z 0 (t)Z 0 β after phaseand-scale correction.With this correction, both the real and the imaginary parts show good agreement with the exact result.Figure 10(b) shows the excitation spectra obtained by Fourier transforming the exact and phaseand-scale-corrected hardware Z 0 (t)Z 0 β at the same points in real time.The excitation spectrum from hardware data accurately reproduces not only the frequencies ω = 0, ±4.90, ±6.37, 7.84 but also the peak amplitudes.
The favorable agreement of the hardware Z 0 (t)Z 0 β with the exact result is in contrast with the deviation of finite-temperature static observables from the exact results in Fig. 8.In fact, the raw hardware Z 0 (t)Z 0 β at t = 0 is 0.821 + 0.397i, which is far from the exact value 1, indicating that phase-and-scale correction has a significant effect in correcting raw hardware data.Even though phase does not enter the static observables we computed, lack of a scale correction scheme for the static observables may explain their large deviation from the exact values compared to dynamical observables.More- over, even though the recompiled circuit in Fig. 10 includes up to 11 gate rounds with the QITE and real time evolution gates combined, the ancilla is initialized after the QITE circuit and hence only experiences 5 gate rounds prior to measurement.The relatively shallow circuit applied to the ancilla may be another reason for the good performance of the quantum device for calculating the finite-temperature dynamical observable Z 0 (t)Z 0 β .

V. CONCLUSION AND OUTLOOK
Our work demonstrates that finite-temperature physics of quantum many-body systems is accessible with near-term quantum hardware.With methods to reduce required quantum sources and mitigate errors in raw hardware data, QITE enables the practical calculation of finite-temperature energies, static and dynamical correlation functions, and spectral densities of excitations.
On two sites, static and dynamical observables for a wide range of temperatures are accurately captured by quantum hardware.An important factor underlying this accuracy is the constant depth of the circuit in both QITE and real time evolution.Constant depth in QITE allowed us to extend QITE-based finitetemperature calculations from a single site [30]; constant depth in real time evolution allowed us to reproduce exact finite-temperature dynamical correlation functions on quantum hardware without phase-and-scale correction as compared to previous studies [14,15].
On four sites, finite-temperature static observables calculated on quantum hardware with circuit recompilation are in reasonable agreement with exact results at β ≤ 0.5.We were also able to accurately reproduce the finitetemperature dynamical correlation function using phaseand-scale correction at a high temperature β ∼ 0.2.However, accurate determination of observables at lower temperatures still appears challenging using the current recompilation scheme where the QITE unitaries are recompiled separately at each imaginary time step.Therefore, treating larger systems and lower temperatures will require additional reduction of circuit depth such as recompilation with merged imaginary time steps [38] or lower error rates on quantum devices either from efficient error mitigation for imaginary time or from improvements in hardware.
Scalable methods need to be employed for systems of larger size.For this reason, we examined how stochastic trace evaluation performs in calculating finitetemperature observables compared to full trace evaluation.We found that on four sites stochastic trace evaluation reproduced full trace evaluation results accurately in the temperature regime we studied.Compared to the previously proposed QMETTS algorithm, stochastic trace evaluation has zero autocorrelation time.A detailed comparison of QITE-based computation of finitetemperature observables with different sampling schemes is certainly a topic worth exploring.Furthermore, with the availability of more qubits [65,66], trading increased computational time due to sampling for an increased number of qubits via constructing density matrix purification states [22,23] may be another feasible direction for studying finite-temperature physics on near-term quantum hardware.
where the quantities with single primes are indexed by µ such that σ µ ∈ N (S) and those with double primes are indexed by µ such that σ µ / ∈ N (S).By setting x[l] to 0, the linear system is reduced to A To show that the set of σ µ can be reduced from N (S) to N (S)/S, suppose σ µ and σ µ belong to the same coset in N (S)/S, then σ µ = ±σ µ s for some s ∈ S. In the QITE unitary e in the sum is a power of −i∆τ times a product of the form ν (x[l] ν σ ν ).If a product term contains x[l] µ σ µ , the action of this term on |Ψ is proportional to In the product over ν , each σ ν ∈ N (S), so Since this applies to every pair of Pauli strings in the same coset, Ĝ[l] can be written as Ĝ where µ is chosen such that σ µ ∈ P Ĥ[l] ∩N (S)/S, x[l] µ = µ η µ x[l] µ , η µ = ±1 and µ is chosen such that σ µ ∈ σ µ S.
Since all Pauli strings on the exponent of e −i∆τ Ĝ[l] commute with elements of S, e −i∆τ Ĝ [l] commutes with elements of S and hence e −i∆τ Ĝ[l] |Ψ ∈ V S .
Our Pauli string reduction scheme is related to the qubit encoding scheme that removes redundant qubits by exploiting Z 2 symmetries reported in Ref. [46].In the qubit encoding scheme, a Hamiltonian over some number of qubits is transformed to another Hamiltonian over a smaller number of qubits by a series of Clifford gates.Our Pauli string reduction scheme coincides with the qubit encoding scheme when the domain size D equals the total number of qubits N , in the sense that the reduced set of Pauli strings in our scheme exactly corresponds to all Pauli strings in the encoded Hamiltonian with redundant qubits removed in the qubit encoding scheme.
However, because the weight of a Pauli string can change during the Clifford transformation, the two schemes differ when D < N .On the one hand, some Pauli strings can decrease in weight after encoding.If we include all Pauli strings with domain size D in the encoded Hamiltonian, these Pauli strings might include those with domain size D > D in the original Hamiltonian, thus increasing the total number of Pauli strings.On the other hand, some Pauli strings can increase in weight after encoding and result in an increased cost of the QITE algorithm.As an example, consider performing QITE on a Hamiltonian with periodic boundary condition and the Z 2 symmetry Z 0 Z 1 Z 2 Z 3 .One of the D = 2 Pauli strings is X 0 Y 3 .In the qubit encoding scheme, the symmetry operator Z 0 Z 1 Z 2 Z 3 is transformed to Z 3 so that qubit 3 can be eliminated, but the weight-two Pauli string X 0 Y 3 is transformed to the higher-weight Pauli string X 0 X 1 Y 2 , thus requiring a larger QITE domain and increasing the overall cost of the algorithm.Therefore, in the present work, we use Z 2 symmetries to reduce the number of Pauli strings in the QITE unitaries rather than eliminate redundant qubits.

APPENDIX B: ERROR BARS IN TRACE EVALUATION
We describe calculation of error bars of finitetemperature observables in full and stochastic trace evaluation.Here we use E(Q) and Var(Q) to denote the mean and variance of a quantity Q.The error in Q is the square root of its variance.
A finite-temperature observable Ô β ≡ O has the expression O = i P i O i / i P i , where P i = || |Ψ i (β/2) || 2 is the (unnormalized) probability and O i = Ψ i (β/2)| Ô|Ψ i (β/2) is the expectation value of the observable after imaginary time evolution on the ith basis state |Ψ i .On quantum computers, each probability P i is built up from the energy expectation value at each imaginary time step: imaginary time evolution (QITE) We begin by reviewing the QITE algorithm in the context of a general Trotterization scheme of the imaginary time propagator.Consider imaginary time evolution on N qubits under a Hamiltonian Ĥ = M m=1 ĥ[m], where each ĥ[m] acts on a local set of qubits.Since the local terms ĥ[m] are not commutative, we need to Trotterize the imaginary time propagator e −β Ĥ by grouping local terms ĥ[m] into Trotter terms Ĥ[l] such that each Ĥ[l] is a sum of local terms ĥ[m] and Ĥ = L l=1 Ĥ[l].For example, for a two-local Hamiltonian where each local term ĥ[m] acts on qubits m − 1 and m, setting L = 2, Ĥ[1] = M/2 m=1 ĥ[2m − 1] and Ĥ[2] = M/2

FIG. 1
FIG.1Quantum circuit to calculate the finite-temperature dynamical correlation function U (t)V β .The ancilla qubit is initialized in |0 and the system qubits are initialized in |Ψ .Measuring X (Y ) on the ancilla yields the real (imaginary) part of U (t)V on the QITE-evolved initial state.Performing a thermal average over all initial states yields U (t)V β .

FIG. 2
FIG. 2 Energy E versus imaginary time β simulated without noise or measurement sampling on a single initial state with and without reduction of the Pauli strings in the QITE unitaries by Z2 symmetries.(a) Four-site TFIM with J = h = 1 and initial state |0001 .The imaginary time step size is set to ∆τ = 0.01.The number of Pauli strings from three D = 2 domains is reduced from 16 to 6 by one Z2 symmetry Z0Z1Z2Z3.(b) Four-site Heisenberg model with J = ∆ = 1 and initial state (|0101 + |1010 ) / √ 2. The imaginary time step size is set to ∆τ = 0.03.The number of Pauli strings on the single D = 4 domain is reduced from 120 to 6 by two Z2 symmetries Z0Z1Z2Z3 and X0X1X2X3.In both panels the energy trajectories using reduced numbers of Pauli strings match the energy trajectories without reduction, which also match the energy trajectories from exact imaginary time evolution.

compared to 16
Pauli strings without reduction by Z 2 symmetries.Figure 2(b) plots energy versus imaginary time of the initial state (|0101 + |1010 )/ √ 2 on the foursite Heisenberg model with J = ∆ = 1.The Hamiltonian and the initial state have two

FIG. 3
FIG.3Four-site recompiled circuit.The four U3 gates at the left constitute the base gate round.Each additional gate round includes a layer of CNOT gates and a layer of single-qubit gates as shown in the dashed box.The additional gate rounds alternate between even-odd and odd-even pairs of qubits, so that the circuit shown consists of three gate rounds.

FIG. 4
FIG. 4 Measurement of a Pauli string in a four-site QITE calculation with post-selection on the stabilizer generator Z0Z1Z2Z3.The appended CNOT gates achieve simultaneous measurement of the Pauli string with the stabilizer generator by transforming Z0Z1Z2Z3 to Z3 acting on a single qubit, from which the stabilizer parity is read off.The other qubits are measured in X-, Y -or Z-basis depending on the Pauli string measured.Measurement outcomes with the wrong parity are discarded.

FIG. 6
FIG. 6 (a) Finite-temperature energy E β and (b) static correlation function X0X1 β of the two-site TFIM with J = ±1, ±3 and h = 1 versus inverse temperature β.The imaginary time step size is set to ∆τ = 0.1.In each of the calculated observables, the mean absolute percentage error between hardware and exact results averaged over all β ranges from 1% to 4%.

FIG. 7
FIG. 7 Finite-temperature dynamical properties of the two-site TFIM with J = 3, h = 1.The imaginary time step size is set to ∆τ = 0.1.(a) Real and (b) imaginary parts of the finite-temperature dynamical correlation function Z0(t)Z0 β at β = 0.2 and β = 1.8 versus real time t.(c) Finite-temperature excitation spectra |S(ω)| 2 versus frequency ω.Positive (negative) frequencies correspond to emissions (absorptions).(d) Amplitudes of the two emission peaks at ω = 7.18 and ω = 5.94.The hardware data capture finite-temperature dynamics of two-site TFIM across a wide range of temperatures.

FIG. 8
FIG. 8 (a) Finite-temperature energy E β and static correlation functions (b) X0X1 β (c) X0X2 β (d) X0X3 β of the four-site TFIM with J = 3, h = 1 versus inverse temperature β with different QITE unitaries.The imaginary time step size is set to ∆τ = 0.05.The D = 2 QITE unitaries are either Trotterized as in Eq. (4) or recompiled, while all D = 4 QITE unitaries are recompiled.The results with recompiled QITE unitaries are closer to exact results than the results with Trotterized QITE unitaries due to circuit depth.Between the calculations with recompiled unitaries, D = 4 is not necessarily closer to exact results than D = 2 for all observables possibly due to the increased influence of hardware noise in the larger linear systems.

FIG. 10
FIG. 10 Finite-temperature dynamical properties of the four-site TFIM with J = 3, h = 1 at β = 0.2.QITE is performed with a time step of ∆τ = 0.05 and recompiled D = 2 unitaries.(a) Real and imaginary parts of the finite-temperature dynamical correlation function Z0(t)Z0 β versus real time t.Raw hardware data are post-processed by phase-and-scale correction.(b) Finite-temperature excitation spectra obtained by Fourier transform of exact and phase-and-scale-corrected hardware Z0(t)Z0 β at the same points in real time.The hardware Z0(t)Z0 β and excitation spectrum after phase-and-scale correction are in good agreement with the exact results.
stochastic trace evaluation, we need to consider initial state sampling on top of measurement sampling.Define the numerator N = n −1 samples n samples i=1 P i O i and the denominator D = n −1 samples n samples i=1 P i so that O = N/D.By first-order expansion of O in N and D and assuming = Φ i (k∆τ /2)| Ĥ|Φ i (k∆τ /2) and n β/2 = β/2∆τ ; each O i is the expectation value of the observable on the QITE-evolved state |Φ i (β/2) .Note that here both the exact imaginary-time-evolved state |Φ i (β/2) and the QITE-evolved state |Ψ i (β/2) are consistent with the definitions in Sec.II B. In full trace evaluation O is regarded a function of the P i and O i , which are random variables because of measurement sampling on quantum computers.Var(O) can be evaluated by expanding O to first order in all P i and O i and assuming all P i and O i are independent, which then gives [E(P i ) 2 Var(O i ) + (E(O i ) − E(O)) i=1 2 Var(P i )]