Robust Extraction of Thermal Observables from State Sampling and Real-Time Dynamics on Quantum Computers

Simulating properties of quantum materials is one of the most promising applications of quantum computation, both near- and long-term. While real-time dynamics can be straightforwardly implemented, the finite temperature ensemble involves non-unitary operators that render an implementation on a near-term quantum computer extremely challenging. Recently, [Lu, Ba\~nuls and Cirac, PRX Quantum 2, 020321 (2021)] suggested a"time-series quantum Monte Carlo method"which circumvents this problem by extracting finite temperature properties from real-time simulations via Wick's rotation and Monte Carlo sampling of easily preparable states. In this paper, we address the challenges associated with the practical applications of this method, using the two-dimensional transverse field Ising model as a testbed. We demonstrate that estimating Boltzmann weights via Wick's rotation is very sensitive to time-domain truncation and statistical shot noise. To alleviate this problem, we introduce a technique that imposes constraints on the density of states, most notably its non-negativity, and show that this way, we can reliably extract Boltzmann weights from noisy time series. In addition, we show how to reduce the statistical errors of Monte Carlo sampling via a reweighted version of the Wolff cluster algorithm. Our work enables the implementation of the time-series algorithm on present-day quantum computers to study finite temperature properties of many-body quantum systems.

Figure 1: Equilibrium properties from quantum dynamics and the noise-induced infrared catastrophe.Bottom panel: Energy spectrum of a 12-qubit Ising chain in a transverse field hx = 0.5.In principle, the spectrum contains sufficient information to compute properties like thermal phase transitions.Top panel inset: A quantum computer or simulator can be used to efficiently compute the Loschmidt amplitude G(t) = ⟨ψ|e −iHt |ψ⟩ which contains information about the local density of states via its Fourier transform.In practice, determining G(t) suffers from shot noise (orange).Top panel: Collecting data from suitably chosen initial states yields the density of states D(w).While it is easy to obtain a reasonable estimate for said density (solid lines), one needs to weight D(w) with a Boltzmann factor (gray dashed) to obtain thermal observables.Any small error induced by shot noise is amplified exponentially (dash-dotted orange vs. blue), leading to significant artefacts when simulating lowenergy physics.The present paper introduces a method to address this problem.
To overcome such challenges, algorithms for quantum computers have been proposed, such as the quan-1 Thermal equilibrium observables from the time-series algorithm All information about a system in thermodynamic equilibrium can be obtained from the expectation value ⟨ Ô⟩ = Tr ρ Ô , (1) where Ô is the operator representing the observable of interest and ρ is the density operator describing the statistical ensemble with the specified thermodynamic constraints.In particular, we are interested in the canonical ensemble, where the temperature of the system 1/β is fixed.This is described by the density operator In the time-series algorithm [28][29][30], thermal expectation values are evaluated using a complete orthonormal basis of states |ψ⟩ where the Boltzmann weights W ψ and local observables O ψ are defined as Because e −β Ĥ is positive definite for any Hamiltonian Ĥ, the Boltzmann weights are non-negative numbers that can be normalized and interpreted as a probability distribution The expression in Eq. ( 3) is then estimated by sampling M states from this distribution via some Monte Carlo method and averaging over the local observable of the sampled states The error of this estimator is O(1/ √ M ).In particular, the error only scales with the inverse square root of the number of independent samples and not the size of the Hilbert space [33].Nevertheless, when samples are correlated, the correlation length can generally depend on the size of the physical system.
In "classical" quantum Monte Carlo methods [34], the Boltzmann weights are calculated by splitting the operator e −β Ĥ into n smaller pieces e −β Ĥ/n and inserting resolutions of the identity in between.The resulting terms are then calculated and summed over also by means of Monte Carlo methods.However, unlike the original Boltzmann weights, the individual terms can generally be negative because they involve transition matrix elements between different basis states ⟨ψ ′ |e −β Ĥ/n |ψ⟩, leading to the infamous sign problem.To avoid this issue, the Boltzmann weights in the time-series algorithm are instead calculated by evaluating their real-time counterpart, known as Loschmidt echos, on a quantum computer G ψ (t) := ⟨ψ|e −iHt |ψ⟩ (9) and then numerically performing a Wick's rotation to the imaginary-time axis t → −iβ.In practice, this is done with the help of the local density of states which is nothing but the Fourier transform of Loschmidt echos Boltzmann weights are then calculated as The last two relations can be readily verified by expressing each of the quantities in terms of the eigenstates of the Hamiltonian |n⟩ and their energies The last equation also explicitly shows the positivity of the Boltzmann weights W ψ (β).Note that as opposed to the rotation from the imaginary-time axis to the real-time axis, known as the analytic continuation problem [35][36][37][38][39][40], obtaining the density of states from Loschmidt echos is a well-defined problem since Fourier transform is just a unitary operation.Nevertheless, obtaining Boltzmann weights from the density of states remains a difficult problem (see Fig. 1), as will become evident later.The numerator of the local observables O ψ can be similarly calculated with the help of a quantum computer.In the special case where the basis |ψ⟩ diagonalizes the observable Ô, its local observables can be directly and efficiently evaluated on a classical computer.In this work, we focus on the latter case.proposed by Refs.[28] and benchmarked in Refs.[29,30] for calculating finite temperature properties with the help of quantum computers.In this work, we provide improvements on the highlighted steps of Monte Carlo sampling and Wick's rotation.
Directly sampling the states from the Boltzmann distribution is generally not possible.
Instead, Metropolis-Hastings algorithm is a general Markov chain method which enables sampling from arbitrary distributions.The price we pay in Markov chain methods compared to direct sampling is that these samples are correlated, and we generally need more samples to achieve the same level of accuracy.The Metropolis-Hastings algorithm can start from any initial state, but it helps to choose a state that already has a large Boltzmann weight.The algorithm then repeatedly proposes a new state |ψ ′ ⟩ from the current one |ψ⟩ using some probability distribution P ψ→ψ ′ .The proposed sample is accepted or rejected using an acceptance ratio If accepted, the current state |ψ⟩ is replaced by the new one |ψ ′ ⟩.Otherwise, the current state |ψ⟩ is used.
In principle, one can use any proposal strategy as long as it guarantees ergodicity, i.e., it allows going from any state to any other state in a finite number of steps.
In practice, one should strive to use a proposal distribution that closely resembles the target distribution, here the Boltzmann distribution, in order to minimize the correlation between the samples.In Fig. 2, we show a flow chart of the full method.The main loop of the algorithm starts by proposing a new state.Then it submits circuits for calculating a time series of its Loschmidt echos to the quantum computer and waits for the resulting shots.From these shots, Loschmidt echoes are estimated and used to calculate the Boltzmann weight of the new state.We then use that to accept or reject the proposed state.Finally, at the end of an iteration, we store the state for later post-processing and calculations of When implementing the time evolution quantum circuits, the depth is reduced by executing the gates corresponding to the bonds in each of the three different orientations in parallel.For example, in the 10-site lattice, the following sets of bonds can parallelized: {(0, 1), (2, 3), (5,6), (7, 8)}, {(1, 2), (3,4), (6,7), (8,9)} and {(0, 9), (2, 7), (4, 5)}.
observables.In the rest of this paper, we will discuss each of these steps in detail, taking the transversefield Ising model as a specific example.

Test Case: Transverse-Field Ising Model
The transverse-field Ising model (TFIM) consists of a set of quantum magnetic dipoles of spin 1/2 on a lattice with nearest neighbor interactions.The interaction is determined by the spin alignment along some axis (usually taken to be the z-axis) in addition to an external magnetic field along a perpendicular axis (usually the x-axis).Its Hamiltonian on a general lattice reads where ⟨i, j⟩ indicates summation over neighboring lattice sites, J is the coupling constant, which we set in the following to J = 1, specifying ferromagnetic interactions.h x is the strength of the external filed, and Ẑi and Xi are Pauli matrices representing, respectively, the z-and x-components of the spin-1 2 operator at lattice site i.This model is one of the simplest quantum spin models exhibiting non-trivial physics that cannot be described classically.On the other hand, its dynamics can be relatively straightforward to implement on gate-based quantum computers, which makes it an ideal candidate for testing the algorithm.The observable of interest is the magnetization at different temperatures β and different values of the external field parameter h x .
The TFIM can be solved exactly in 1D, where it shows a phase transition at zero temperature but no order at finite temperatures [41].In higher dimensions, the phase diagram is divided by a secondorder critical line into a ferromagnetic ordered region with ⟨ Ẑ⟩ ̸ = 0 and paramagnetic disordered one with ⟨ Ẑ⟩ = 0 [42].The endpoints of the critical line are characterized by the critical temperature β c of the classical Ising model (i.e., h x = 0) and a critical value of the field parameter at zero temperature.In this work, we use a two-dimensional honeycomb lattice with 10 and 16 sites (see Fig. 3).The classical critical temperature of this lattice is 3) 2 in the thermodynamic limit [43,44].As a sampling basis set |ψ⟩, we choose the product states of the Z i operators.These are the easiest states to prepare on a quantum computer and can be directly implemented using X gates.Besides, they are the eigenstates of the magnetization operator; therefore, we can easily calculate their local magnetization on a classical computer and only need the quantum computer for evaluating Loschmidt echos.

Quantum Circuits
Calculating the real and imaginary parts of Loschmidt echos can be achieved using the so-called Hadamard test.This involves two quantum circuits, one for the real part and one for the imaginary part, where each circuit contains one ancilla qubit and a time evolution circuit controlled by the ancilla as shown in Fig 4 .By measuring the ancilla qubit, the desired real and imaginary parts can be calculated as the expectation value of the observable Ẑancilla , i.e., the difference between the probability of measuring the value 0 and the probability of measuring the value 1.Moreover, by additionally measuring Ẑancilla ⊗ Ô, i.e. correlating a measurement of Ô on the system qubits with a measurement of Ẑ on the ancilla, we can obtain ⟨ψ| Ôe −i Ĥt |ψ⟩ additionally to the Loschmidt echos from the same shots [28].This enables the evaluation of Eq. (6) for arbitrary observables.As mentioned before, here we focus on the case where this is not necessary, since |ψ⟩ are eigenstates of Ô.
The time evolution circuits are obtained via secondorder Trotterization.
where the Hamiltonian is split into two noncommuting parts Ĥ = Â + B. The first part Â contains the Xi terms, while the second B contains the mutually commuting Ẑi Ẑj terms.The time evolution of each part can be achieved easily using one-and twoqubit gates.In particular, we add one R ZZ (−2Jt) for each pair of interacting sites and an R X (2ht) gate for each site.To reduce the depth of the time-evolution circuits, it is important to order the Ẑi Ẑj terms such that bonds with disjoint sites can be executed in parallel.For the honeycomb lattice, this is achieved < l a t e x i t s h a 1 _ b a s e 6 4 = " f k X 8 c f u 0 t t 6 m 1 5 0 / 4 f F H g e S p a W R S l X V 3 d N V V F p 5 j u P f n e j S 5 S t X r + 1 d 7 9 6 4 e e v 2 n d 7 x 7 1 t r T p x Z t 9 s q q y 8 8 o Z 4 5 Y f w Z / g 3 c 9 J M r B t X s n 1 8 7 r 3 H 9 n F a S G E w D P 9 6 / q 3 b d + 7 e 2 7 g f P H j 4 6 by grouping the bonds along each of the directions 0, π/3, 2π/3 (e.g., see the caption of Fig. 3).The Hadamard test requires controlling the time evolution circuit, which is achieved by controlling the R ZZ and R X gates.To control an R ZZ gate, we first decompose it into an R Z gate sandwiched between two CNOT gates and then control only the R Z gate.CNOT gates do not need to be controlled because they always come in pairs.The resulting controlled R X and R Z gates are implemented each using a single R ZZ gate, which is native to ion-trapped quantum computers, as shown in Fig 5b .Consequently, a single time slice of Â requires a number of two-qubit gates that equals the total number of lattice sites, while implementing a single slice of B requires three times the number of bonds in the lattice.An example of the controlled time-evolution circuit is given in Fig. 5a.
To reduce the number of two-qubit gates in the controlled circuit, we rewrite the second-order Trotterization as following By implementing e −i Ât/2n as part of the state preparation circuit, Â terms appear only n times in the controlled time evolution circuit instead of the n + 1 times needed using Eq. ( 16) .In total, implementing a single controlled second-order Trotter step requires 43 two-qubit gates for the 10-site lattice and 73 for the 16-site one.

Monte Carlo Sampling
The most straightforward strategy for sampling product states of the TFIM is the single spin flip update.This strategy picks a site randomly from the lattice < l a t e x i t s h a 1 _ b a s e 6 4 = " I W z c y   and reverses its spin.The proposal probability is symmetric and equals P ψ ′ →ψ = P ψ→ψ ′ = 1/L.Therefore, the acceptance ratio is determined solely by the Boltzmann weights: However, the old and new states differ only by one spin, which means that this local strategy is slow in exploring the Hilbert space, leading to a high correlation at low and moderate temperatures.
A more effective approach that makes global moves is Wolff's cluster update algorithm [32].It involves randomly selecting a cluster of aligned spins and flipping all of the spins in the cluster at once.The cluster is constructed by choosing a random lattice site (the cluster's root) and adding with probability 1 − e βJ each of the neighboring spins that point in the same direction.The parallel neighbors of the newly added sites are again added with that probability, and the cluster keeps expanding stochastically until there are no new parallel neighbors.An example on a square lattice is given in Fig. 6.
This algorithm is designed to sample from the Boltzmann distribution of the classical Ising model, and it was shown to be efficient even near the critical temperature [32].Here, we adapt this algorithm by using it as a proposal strategy within the Metropolis-Hastings algorithm.To this end, we need the ratio of the proposal probabilities, which is equal to the ratio of the classical Boltzmann weights of the new and old states (see Eq. ( 6) of Ref. [32]).Given that the energies of the product states in the classical Ising model are the same as those of TFIM, this ratio can be written as  which involves only energies of the product states and can be calculated efficiently on the classical computer.We want to emphasize here that using the cluster update as a proposal strategy does not affect the correctness of the Monte Carlo simulations.As a direct application of the Metropolis-Hastings algorithm, these simulations satisfy the detailed balance condition for the quantum TFIM.In the Metropolis-Hastings algorithm, the detailed balance condition is satisfied as long as the proposed samples are accepted with the proper acceptance ratio (Eq.( 14)).One is then free to choose any ergodic proposal strategy as long as the correct ratio of proposal probabilities is used in the acceptance ratio.When using the cluster update, the acceptance ratio reads In Fig. 7, we compare the single flip update with the cluster update for the 16-site honeycomb lattice with h x = 0.5 at the classical critical temperature β c .The plot demonstrates how the cluster update thermalizes almost instantly, while the single flip update needs hundreds of iterations and suffers from much longer correlation times.The efficiency of the cluster update depends on the strength of the field h x and inverse temperature β, because the classical Boltzmann distribution, which is used to propose samples in the cluster update, deviates more from the target quantum Boltzmann distribution at lower temperatures and higher h x .This can be seen by expanding classical and quantum Boltzmann weights as a function of inverse temperature When h x ̸ = 0, the classical and quantum moments of product states do not match ⟨ψ| Ĥ|ψ⟩ n ̸ = ⟨ψ| Ĥn |ψ⟩ for n ≥ 2, and the difference is proportional to h n x .Consequently, sampling with the cluster update becomes less efficient for lower temperatures and higher h x .

Wick's Rotation
Wick's rotation is the most crucial component of the time-series Monte Carlo algorithm.The procedure for performing it is formally straightforward: Compute the local density of states D ψ (ω) via Fourier transform of Loschmidt echos [cf.Eq. (11)], then integrate it with e −βω factor to get Boltzmann weights [cf.Eq. (12)].In practice, however, obtaining Boltzmann weights reliably is difficult, especially for lower temperatures.The problem comes from the exponential factor e −βω , which makes Boltzmann weights sensitive to the errors made in estimating the local density of states.Therefore, we need to find a way to obtain the Fourier transform with a high accuracy using only a limited number of noisy data points.
In this section, we first discuss the issue of only having a limited number of data points using the method of Gaussian filtering introduced by Ref. [30].We propose an automatic way of determining its parameters and study its error behavior.We then introduce the non-negative least squares (NNLS) method, which gives superior results, especially for a small number of data points.Next, we address the issue of statistical shot noise which is present in any quantum computer (even a fault-tolerant one) and show how to stabilize NNLS in the presence of this noise.Finally, we benchmark our recipe using the 10-site TFIM.

Effect of Time-Series Truncation
Given that the spectrum is bounded, then according to the Nyquist-Shannon theorem [45], we only need to calculate Loschmidt echos at discrete times spaced at an interval The spectral bound max {|E i |} depends on the system under study.For example, the 10-and 16-sites TFIM with h x = 1 has the following values 13.4478 and 22.6223, respectively.Obtaining the exact value of this bound requires solving for the lowest and highest eigenenergies of the Hamiltonian Ĥ, which is a difficult problem in general.However, in many cases, one can already obtain an upper bound analytically.For example, for the TFIM, we find the following upper bound where N bonds is the number of lattice bonds and N sites is the number of lattice sites.
Although it is sufficient to sample Loschmidt echo at the finite rate of Eq. ( 22), we still need, in principle, to evaluate them at an infinite number of times points.Otherwise, truncating the time series at a maximum time T max is equivalent to convoluting its Fourier transform with a sinc function An illustrative example is shown in Fig. 8.This convolution is problematic for three reasons.On the one hand, the new truncated density D Tmax ψ is not band limited, so we need to choose δt smaller in order to evaluate its values at higher |ω|.On the other hand, even with a full knowledge of D Tmax ψ , the integral of the Boltzmann weights will not converge because 1/ω decays much slower than the exponential increase of e −βω for large negative frequencies.Finally, even for integrals that do converge analytically, e.g., when β = 0 or for the microcanonical algorithm [28], the oscillations of the sin function make detecting the convergence numerically more challenging.

Gaussian Filter Method
To circumvent these issues, Ref. [30] proposed to multiply the time series by an (unnormalized) Gaussian function of width This is equivalent to convolving its Fourier transform with a Gaussian of the reciprocal width δ  The bottom panel shows in blue its Fourier transform when sampled at a rate 1/∆t = 64/π and truncated at Tmax = 4π.The exact density of states, a linear combination of delta functions located at the eigenenergies, is plotted in black.For ease of visualization, the delta functions are replaced by Gaussians whose width equals the grid spacing ∆ω = 1/8.Note the negative component of the density of states obtained from the truncated time series.
Inserting this ansatz into Eq.(12), we find that the effect of this broadening on the Boltzmann weights is simply a multiplication by a scale factor W ψ,δ (β) := dω e −βω D ψ,δ (ω) = e β 2 δ 2 /2 W ψ (β).(27) Since this scale factor is independent of the specific state |ψ⟩, the normalized Boltzmann weights p ψ stay invariant.While the broadened density in its exact form produces an equivalent result, it is more amenable to approximations.For large enough δ ≫ 1/T max , the Gaussian filter smoothens out the oscillations in the region of interest in D ψ,δ (ω), which allows cutting off the remaining artificial oscillations in the tail.Ref. [30] proposes setting to zero all density values below a certain positive cutting threshold D cut ψ .However, it does not discuss how to choose such a threshold or what error the truncation introduces in the resulting Boltzmann weights.
Ideally, we want to choose the minimum cut such that it filters out the artificial oscillations but nothing more.Determining this value exactly is hard because it requires disentangling the original peaks from the artificial ones.Instead, we propose, as a heuristic, to use the maximum absolute value of the negative part of the density to determine the magnitude of those Figure 9: Schematic plot of premature truncation and interference effects.Using a finite sampling frequency π/∆t, the Fourier transform is a periodic summation of the original density (green) with period 2π/∆T .As δ increases, the broadened Fourier transform (red) inside the window [−π/∆t, π/∆t[ increasingly deviates from the true broadened density (blue) due to interference with shifted copies outside this window.The finite sampling frequency π/∆T , thus, puts a practical limit on the width of the Gaussian filter.
oscillations.We then choose the cut proportional to this value We found that using C cut = 2 gives stable results in all examined cases.The next question to address is choosing the parameter δ given a fixed maximum time T max .If we know the truncated time series analytically (i.e., knowing its values to infinite precision at all times prior to T max ), then higher values of δ are always better, and the truncation error of the Gaussian filter will decrease monotonically.In practice, however, there are two limitations.The first limitation is that, unlike the original density D ψ , the broadened one D ψ,δ is not band limited but decays as a Gaussian.Therefore, one must choose a large enough ω window outside which the weight is negligible.This window increases as δ increases, and thus to maintain the same level of accuracy, we need to decrease the time step ∆t.Failing to do so will lead to premature truncation of the density and an interference with shifted copies of itself (see Fig. 9).
The second limitation is the precision of floatingpoint arithmetic.For small enough 1/δ, the tail of the Gaussian filter becomes numerically zero, such that the filtered time series G ψ,δ is effectively truncated even further instead of being reweighted.Therefore, a limit to α := δT max is determined by the floating-point precision.To show this effect, we plot in Fig. 10a the relative error in estimating Boltzmann weights as a function of δ for the time series of Fig. 8.The error decays exponentially with increasing δ up to a point α = 8, after which the error plateaus or increases (depending on temperature).In Fig. 10b, we show the exact same calculation with single-precision floating numbers.In this case, the optimal value of α reduces to α = 5.The increase in the error after that point  is attributed to the aforementioned reduction in the effective maximum time T max .The dependence of the relative error on T max while fixing α is shown in Fig. 11a.From this plot, we see that the error decreases exponentially with increasing 1/T max until it plateaus when hitting numerical accuracy.Finally, in Fig. 11b, we show the error as a function of inverse temperature β, which exhibits a similar exponential increase with β.The dependence of the error on the different parameters can be summarized in the limit of large α as follows: where C is some positive constant.This expression, derived from some simplified assumptions in Appendix A, is only valid when β/T max < C and matches the observed behavior in the previous plots.An interesting observation is the duality between temperature and maximum time, where lowering the temperature requires increasing the maximum time proportionally to maintain the same level of accuracy.Another important observation is that using a fixed α = δT , the  error does not vanish even in the limits T max → ∞ or β → 0. However, by choosing the appropriate value of α, values close to the numerical accuracy can be reached, as evident from Figs. 11a and 11b.
It should be emphasized that the results discussed so far used a sampling frequency much higher than the Nyquist frequency in order to minimize the aforementioned interference effects on the broadened density.When the sampling frequency is relatively low, we cannot set δ to its optimal value T max /α.Instead, it should be chosen such that δ ≪ π/∆t − max {|E i |}, which increases the error significantly as illustrated in Fig. 12.In practice, the need for high sampling frequency puts yet an additional burden on the quantum resources required for an accurate reconstruction of Boltzmann weights.

Non-Negative Least Squares
In the previous section, we have shown how to systematically reduce the error of the Gaussian filter method by increasing both the maximum time and the sampling frequency.When either is low, however, it remains problematic to extract Boltzmann weights, especially at relatively low temperatures.Relative Error Figure 12: Same calculations as Fig. 11 but using time series sampled at the lower rate 1/∆t = 16/π.To limit the effect of this low sampling rate on the broadened density (cf.Fig. 9), a cap is put on the filter parameter: δmax = (π/∆t − E0)/2, where E0 = −13.477(8) is the ground state energy.Note that curves of Tmax = π and Tmax = 2π are on top of each other.
In the following, we introduce the main contribution of our work: a method for obtaining much better results using the same limited data set.The critical idea is using the fact that physical densities of states must be non-negative.So instead of taking the time series and directly computing its Fourier transform, we look among the non-negative densities for the one whose inverse Fourier transform best fits the time series.This can be found using the non-negative least squares (NNLS) method, written formally as: where χ 2 measures how well the density fits the time series: where n t is the number of time points.In practice, we construct a frequency grid in the interval [−π/∆t, +π/∆t[ with enough points to resolve the details of the density, and we use it to build the inverse Fourier transform matrix M i,j = e −iωj ti .Then we solve the linear system G = MD under the constraint   8 using non-negative least squares method.For ease of comparison, both the exact density and NNLS are convoluted with a Gaussian whose width equals the grid spacing ∆ω = 1/8.D ≥ 0, where G and D are vector representations of the time series and the density, respectively.The solution is obtained numerically using the active set algorithm by Lawson and Hanson [31].
In Fig. 13, we show the NNLS solution of the time series in Fig. 8, and the result matches the exact density on the specified ω grid.This is a significant improvement over the direct Fourier transform and translates into several orders of magnitude reduction in the error when calculating Boltzmann weights.In Fig. 14, we show the relative error of NNLS as a function of maximum time and temperature.These results were obtained using the lower sampling rate 1/∆t = 16/π and are superior to the ones obtained from the Gaussian filter method using the same data (cf.Fig. 12).
When obtaining Boltzmann weights from NNLS densities, we recommend to additionally use a very mild Gaussian filter of width proportional to the spacing of the frequency grid.Without it, NNLS has to choose between two neighboring frequencies when placing a delta function (or a very narrow peak).This makes low-temperature results sensitive to how close one of the grid points is to the exact location of the delta function.Using the suggested filter, the broadened density is easier to represent on the discretized grid.It should be emphasized that, unlike in the original Gaussian filter method, the role of the broadening is only to regularize the effect of frequency discretization and can be made arbitrarily small using a finer ω grid.

Influence of shot noise on NNLS: Quantile Filter and Discrepancy Principle
The discussion so far has been about limited but numerically exact data.Using quantum computers, however, the time series are estimated by averaging over a finite number of shots from the quantum circuits.As a result, the time series contains statistical noise that scales as 1/ √ N shots .In particular, the shot noise on the real and imaginary parts of Loschmidt Relative Error

echo G have zero mean and variances
N shots , The effect of this noise on the NNLS density is introducing small, randomly placed peaks with a total weight proportional to the noise level as illustrated in Fig 15 .When these spurious peaks occur at a frequency lower than the first true peak, they introduce a disproportionally larger error in the Boltzmann weights, especially at low temperatures.The first step in amending the effect of noise is redefining the fit function χ 2 of NNLS in Eq. (31) such that each term is weighted by the inverse of its variance: This drives NNLS to fit accurate data points better than noisier ones.For example, Loschmidt echo at  8 using non-negative least squares method.The shot noise is approximated as Gaussian noise of zero mean and variances given by Eq. (32).Here N shots = 1000 is used.t = 0 is known exactly and corresponds to the normalization constraint G(0) = dωD ψ (ω) = 1.This constraint can then be approximately enforced during NNLS by using a relatively large weight for G(0) terms in the fit function.
To limit the effect of the noise on the Boltzmann weights, we introduce, as a post-processing step, a quantile filter where the density outside a specific quantile range [q, 1 − q] is truncated.The quantile of a normalized density function is defined as ) is its cumulative distribution function (CDF).If value of q is chosen appropriately, leading and trailing spurious peaks will be eliminated.The main advantage of the quantile filter, as opposed to the simple cutting procedure used by the Gaussian filter, is that the proper value of q is relatively stable regarding the broadening of the density.
The optimal truncation value q can be determined automatically using the discrepancy principle [46,47].According to this principle, the truncation level is chosen such that the truncated density Dq ψ (ω) fits the time series up to the expected amount of noise on that series.When N shots is large, the shots noise is approximately Gaussian, and thus the reweighted fit function of Eq. (33) follows the chi-squared distribution with 2n t degrees of freedom, where n t is the number of time points.The mean value of this distribution equals its degrees of freedom, so the discrepancy principle suggests choosing the quantile value q such that χ 2 [ Dq ψ (ω)] = 2n t .A larger value can be used to avoid accidental over-fitting of the noise or if other sources of systematic error are expected (e.g., Trotter error).
In Fig. 16, we plot the relative error in Boltzmann weights as a function of the number of shots for the time series of Fig. 8.We simulated the effect of shot noise using Gaussian random variables of appropriate variances given by Eq. (32).The density was then reconstructed using NNLS and truncated with a quantile filter, where the truncation level was determined automatically using the discrepancy principle.The plot shows that using the quantile filter in combination with the discrepancy principle effectively regularizes the impact of shot noise and allows systematically improvable results.The error on the Boltzmann weights scales linearly with the error on the time series, and thus it scales as 1/ √ N shots .In particular, no additional hyperparameter is introduced as the algorithm only takes the number of shots as an additional input parameter.

Benchmarking NNLS
To test the overall effectiveness of our Wick's rotation strategy (NNLS + quantile filtering + discrepancy principle), we calculate the exact Boltzmann weights W ⋆ ψ for all the states of the 10-sites TFIM and compare them with the weights W ψ estimated via Wick's rotation of the time series of these states.As a global measure of the error made during Wick's rotation, we need a statistical distance quantifying how different the estimated Boltzmann distribution One option, which we employ here, is the relative entropy (also known as Kullback-Leibler divergence) of the estimated Boltzmann distribution with respect to the exact one It is worth noting that we also checked the L 1 -distance (also known as Kolmogorov distance) and it shows similar trends to the entropy, so we restrict the discussion below to the later.In Fig. 17, we plot the error for different values of N shots at h x = 1 and different β.We find a slight improvement of these results when enforcing constraints on the first few moments of the density, in this case, the mean energy E ψ := ⟨ψ|H|ψ⟩ = dωD ψ (ω)ω, and the energy variance ⟨ψ|(H − E ψ ) 2 |ψ⟩ = dωD ψ (ω)(ω − E ψ ) 2 .These moments can be calculated efficiently and accurately on  the classical computer for product states.The truncation performed by the quantile filter will generally violate these constraints.However, they can be reintroduced by shifting and rescaling the truncated density to fix its moments to their exact values.The plot shows the relative error in Boltzmann weights using adjusted moments which makes a slight but noticeable improvement for low numbers of shots and high temperatures.Fig. 18 shows the number of shots required to reach a specific level of accuracy in the Boltzmann distribution for different temperatures.The number of shots scales exponentially with the inverse temperature.This exponential scaling comes from the fact that errors in the estimated density gets amplified by the Boltzmann factor e −βω which itself scales exponentially with the inverse temperature.
In Fig. 19, we also show the dependence of the relative entropy error on the strength of magnetic field parameter h x .As h x increases, Wick's rotation becomes more difficult.This can be explained by the fact that the density of states gets broader with increasing h x (the energy variance equals h 2 x times the number of sites).With wider support, more of the shot noise is located inside the main body of the den- sity and is not removed by the quantile filter, which cuts only the tails.This observation suggests that sampling entangled states (instead of simple product states) with lower energy variances can lower the effect of noise on Boltzmann weights.Note that in the limit when the sampled states are eigenstates (e.g.h x = 0), moments-adjusted densities can become exact even in the presence of shot noise.This is a result of enforcing the moments of the density by hand and contingent on using exact values of the eigenenergies.
If the moments themselves, however, contain errors, then these errors will be reflected in the Boltzmann weights and scale exponentially in inverse temperature.For example, let ψ 1 , ψ 2 be two eigenstates with energies E 1 , E 2 , respectively.Then the ratio of their Boltzmann probabilities is e −β(E2−E1) .If an error ϵ is made in estimating E 2 − E 1 , then the relative error in the ratio is |1 − e −βϵ |, which is exponential in inverse temperature.Finally, we investigate the effect of Trotterization on Wick's rotation.We expect that once the Trotter error is below the statistical shot noise, Trotterization has no noticeable effect on Boltzmann weights.We used a Trotter step δt := ∆t/n Trotter , where n Trotter represents the number of Trotter steps per sampling period ∆t.The results are plotted in Fig. 20 for different numbers of shots and confirm our expectations.For a high number of shots, the time series is accurate enough such that increasing the number of Trotter steps systematically decreases the error in Boltzmann weights.On the other hand, using a relatively low number of shots using one Trotter step per sampling period is already enough.

Simulation Results
To demonstrate the robustness under the combined effects of Trotter error, shot noise and Monte Carlo  sampling errors, we implemented the whole algorithm using quantum circuits simulated via the IBM's Aer-Simulator with no hardware errors.We performed calculations for the 16-sites TFIM with h x /J = 1 at different temperatures.Loschmidt echoes were evaluated on n t = 8 times points equally spaced between zero (excluded) and JT max = 1 (included).We used a variable Trotter step with a maximum value of 0.25.This implies that the first two time points use one Trotter step, the third and fourth points use two steps and so on.We sampled the output of each circuit using 10000 shots.With this number of shots, Trotter error is still appreciable compared to the shot noise.
To regulate its effect on the density D(ω), we applied the discrepancy principle with a higher value χ 2 [ Dq ψ (ω)] = 5n t .We used four Monte Carlo chains for each temperature value, each with 512 samples and a different random seed.Two chains started from the all-spin-up state, and the other two started from the all-spindown state.We discarded the first 32 samples of each chain as a burn-in period.To save on the number of executed quantum circuits, we implemented a caching mechanism that allows reusing the shots when the same quantum state is proposed again.The savings gained from this trick depends on the length of the Markov chain and the temperature.At lower temperatures, most samples are low energy states, while at high temperatures, the Boltzmann distribution is more uniform with a lower probability of repeated sampling.The ratio of unique states for β c /3 chains was around 94%, for β c was around 63%, and for 8β c /3 was around 23%.
In Fig. 21, we show how the average squared magnetization evolves with iterations of different Monte Carlo chains at β c .The exact result falls within two standard deviations from the final estimated average.On the other hand, the classical result (i.e., for h x = 0) is well outside the error bars, so we could not have obtained these results by simply using the classical Boltzmann weights e −β⟨ψ| Ĥ|ψ⟩ .This comparison is relevant because the quantile filter truncates the tails of the density of states, making it look closer to the density of states of the classical Ising model (which is a delta function at the mean energy).Therefore, in the presence of noise, the estimated quantum Boltzmann weights tend to get biased towards the classical Boltzmann weights, and the results here show that, at this temperature, the truncation does not significantly bias the estimates towards the classical model.
In Fig. 22, we compare the estimated values of squared magnetization using AerSimulator (plotted in green) with the exact ones (plotted in black) at inverse temperatures β c /3, 2β c /3, ..., 8β c /3.We see that the bias is limited to 7.4% even below the critical point at inverse temperature 4β c /3.As expected, the bias gets more prominent as the temperature decreases.To disentangle the source of this bias, we also performed simulations using Loschmidt echos obtained numerically via QuSpin Python package.These Loschmidt echos were calculated at the same time points as in the Aer simulations, but they are numerically exact, i.e., without Trotter error or shot noise.Accordingly, the Boltzmann weights were estimated using NNLS without requiring the quantile filter.The results using these exact Loschmidt echos (plotted in blue) show no observable bias beyond the statistical error bars.This indicates that the bias in AerSimulator results can be mainly attributed to Trotter error and the shot noise.Reducing the bias involves taking more shots and smaller Trotter steps.We also notice that statistical error bars get larger for lower temperatures.This is explained by the lower efficiency of the cluster update at lower temperatures, which has been discussed earlier at the end of Sec. 4. Reducing these error bars is a matter of taking more Monte Carlo samples.

Summary and Outlook
Time-series quantum Monte Carlo is a promising hybrid algorithm for calculating the thermal properties of quantum materials.The algorithm relies on calculating quantum Boltzmann weights via Wick's rotation of real-time dynamics to imaginary time.In this work, we studied in detail how to calculate Boltzmann weights from noisy and truncated time series.We revisited the Gaussian filter method and discussed an automatic way of determining its cut parameter.Under simplified assumptions, we derived an asymptotic of the error made on Boltzmann weights and found that lower temperatures require proportionally longer times to maintain the same level of accuracy.We then proposed to use the non-negative least squares method, which enforces the non-negativity of the density of states.We found it to reduce the error in the Boltzmann weights by several orders of magnitudes compared to the Gaussian filter method using the same data.Numerical results suggest that the duality between maximum time and temperature still holds using this method.
To regularize the effect of shot noise, we suggested using a quantile filter for truncating noisy tails of the density.The amount of truncation is determined automatically via the discrepancy principle.Benchmarking results show that the number of shots required to maintain a specific level of accuracy scales exponentially with inverse temperature.It also suggests that entangled states with lower energy variances can be less susceptible to shot noise.To test the robustness of the recipe in realistic settings, we fully implemented the time-series quantum Monte Carlo algorithm for the two-dimensional TFIM, including quantum circuits and an improved sampling strategy based on the cluster update algorithm.The simulations indicate that, despite the aforementioned sensi-tivity to shot noise, the algorithm is able to obtain relatively accurate results up to the critical temperature using affordable resources.Hardware errors aside, the parameters used in these simulations are within the capabilities of the current hardware.This demonstrates that the method can be of practical use on current NISQ devices.
While the ferromagnetic TFIM model investigated here is relatively easy to solve on classical computers, the improvements made in Wick's rotation are general.In particular, we expect the method to be useful for frustrated quantum spin models, e.g., the antiferromagnetic Heisenberg model on a kagome lattice, which is hard to simulate classically.Also, the proposed reweighted cluster update can be similarly applied in that situation by proposing samples from an efficient classical sampling algorithm, e.g. the KBD algorithm [48,49], and reweighting with the ratios of quantum to classical Boltzmann weights.Whether the reweighted algorithm remains efficient at relevant low temperatures and large system sizes remains to be investigated.

A Error Scaling in Gaussian Filter Method
There are various sources of error in the estimation of Boltzmann weights in the Gaussian filter method.First, there is an interference error resulting from the finite sampling frequency π/∆t because the broadened density has infinite support (see Fig. 9).There is also a quadrature error associated with performing the integral numerically with a finite number of ω points.These two sources of error are controlled by using finer time and frequency grids, and we assume that ∆t and ∆ω are chosen small enough to make the resulting errors negligible.More importantly, there is an error introduced by cutting values below D cut ψ and errors on the other values remaining from having a finite maximum time T max .These errors are dominated by the error resulting from setting the values at the lowest frequencies to zero.Let ω cutoff be the smallest frequency below which all values of the density are truncated.That error can then be calculated as To simplify the analysis, we focus on the case where the exact density is a single delta function at frequency ω 0 and discuss implications for the general case later.In this simple case, the broadened density is a Gaussian and the exact value of its Boltzmann weight is W ψ,δ = exp(β 2 δ 2 /2 − βω 0 ).The relative error can then expressed in terms of the error function as following: The distance between the peak and cutoff frequency ω cut − ω 0 will depend on the maximum time T max and broadening parameter δ and can be estimated as follows.The density of the truncated time series reads where erfi is the imaginary error function and we de- The first term gives the exact density of the full time series, while the second one approximates the error made due to the time-series truncation.The oscillations start to dominate when the exponential factors e a 2 and e b 2 have similar magnitudes.This allows us to infer how the frequency cut ω cut scales where C is some positive constant of order 1.Substituting back in Eq. (36), the relative error reads (40) In the same limit as above, δ T max → ∞, and assuming β < CT max , we can expand the error function erf(z) at negative infinity to order O(z −3 e −z 2 ) and approximate the relative error as (41) Setting α := δ T max , Eq. ( 29) is obtained.Now we discuss the general case of a linear combination of delta functions.Let µ be the mean of the density and σ is its standard deviation.Assuming δ is much larger than the spacing between the lowest and highest frequencies, for the sake of error scaling, the broadened density could be approximated by a Gaussian of width σ+δ and mean µ.With this assumption, the previous error analysis holds by setting ω 0 to the mean µ and replacing δ by the extended width σ + δ.
ψ⟩ ← | ψ 0 ⟩ Post-Processing < l a t e x i t s h a 1 _ b a s e 6 4 = " c e 5 L N 1 8 0 b J w 1 a h A 8 r A K C 9 G J c D v 8 = " > A A A C U H i c b Z B N T x s x E I Z n U 1 p o + h X a I x e L q F J P 0 W 6 F a A 9 F Q n D p B U G l B p C y 6 W r W 8 R I r 3 r V l z 1 a N r P 0 z / J p e y 6 3 / p B c E z k c l A o x k + d U 7 M x 7 P k x s l H c X x 3 6 j 1 Z O 3 p s / W N 5 + 0

Figure 2 :
Figure2: Flowchart of the time-series Monte Carlo algorithm proposed by Refs.[28] and benchmarked in Refs.[29,30] for calculating finite temperature properties with the help of quantum computers.In this work, we provide improvements on the highlighted steps of Monte Carlo sampling and Wick's rotation.

Figure 4 :
Figure 4: Hadamard test circuits for calculating the real (top) and imaginary (bottom) parts of Loschmidt echo.The box e −i Ĥt represents the time evolution circuit (see Fig.5a), H is the Hadamard gate and S † is the adjoint of the phase gate.
x 2 P b s P L o 7 7 w d f + l / P P 3 e H 3 1 s 8 9 8 o a 8 I z 0 S k A E Z k h / k j I w J J 3 f O o U O d t + 7 A j d y Z m 2 6 k r t P W v C Z b 4 a o / 3 u v R L g = = < / l a t e x i t > RZZ( ✓/2) RZ(✓/2) < l a t e x i t s h a 1 _ b a s e 6 4 = " q d D E 8 X a b 6 5 u S N C i 7 l w z 9 O 6 Z C 9 w 1 S s / F D e 9 l P m Z J m 0 b F z J 9 t E 5 9 8 M + T n I p D I b h P 8 9 / 8 v T Z 1 v P t n e D F y 1 e v d z t 7

Figure 5 :
Figure 5: (a) Example of the controlled time-evolution circuit for the transverse field Ising model.Here is shown a single second-order Trotter step for a 4-sites ring.(b) Decomposition of the controlled RZ and RX gates using RZZ gate and other single-qubit gates.

Figure 6 :
Figure 6: Example of a cluster update step taken by Wolff's algorithm.Arrows represent sites of spin up or down, and black ones denote those already added to the cluster.Dashed boxes designate the most recently added sites, while gray boxes designate new neighbors with the same spin as the cluster.

Figure 7 :
Figure 7: Evolution of the average squared magnetization using different proposal strategies.Each point represents the average of the chain up to the specified iteration.The first 50 iterations (burn-in period) are excluded.Both chains start with the all-spin-down state.

Figure 8 :
Figure 8: Illustration of the effect of truncating the time series on its Fourier transform.The top panel shows the time series of the all-spin-up state of 10-sites TFIM with hx = 1.The bottom panel shows in blue its Fourier transform when sampled at a rate 1/∆t = 64/π and truncated at Tmax = 4π.The exact density of states, a linear combination of delta functions located at the eigenenergies, is plotted in black.For ease of visualization, the delta functions are replaced by Gaussians whose width equals the grid spacing ∆ω = 1/8.Note the negative component of the density of states obtained from the truncated time series.

Figure 10 :
Figure 10: Dependence of the relative error on the Boltzmann weights of the Gaussian filter method on the product δTmax for the time series of Fig. 8.The time series was sampled at a rate 1/∆t = 64/π up to maximum time Tmax = 4π.Calculations were performed using double-precision (top panel) and single-precision (bottom panel) floating-point representations.

Figure 11 :
Figure 11: Dependence of the relative error of the Gaussian filter method for estimating the Boltzmann weights on the maximum time Tmax (top panel) and inverse temperature β (bottom panel) for the time series of Fig. 8.The calculations were performed using double precision.The time series were sampled at a rate 1/∆t = 64/π, and the filter parameter δ was chosen such that α = δTmax = 8.

Figure 13 :
Figure 13: Fourier transform of the time series of Fig.8using non-negative least squares method.For ease of comparison, both the exact density and NNLS are convoluted with a Gaussian whose width equals the grid spacing ∆ω = 1/8.

Figure 14 :
Figure 14: Dependence of the relative error of the nonnegative least squares method on the maximum time Tmax (top panel) and inverse temperature β (bottom panel) for the time series of Fig. 8.The time series were sampled at a rate 1/∆t = 16/π as in Fig. 12.

Figure 15 :
Figure15: Fourier transform of a noisy version of the time series of Fig.8using non-negative least squares method.The shot noise is approximated as Gaussian noise of zero mean and variances given by Eq.(32).Here N shots = 1000 is used.

Figure 16 :
Figure 16: Dependence of the relative error of the nonnegative least squares method on the number of shots for the time series of Fig. 8.The time series were sampled at a rate 1/∆t = 16/π with maximum time Tmax = 4π.

Figure 17 :Figure 18 :
Figure 17: Dependence of the relative entropy error on the number of shots for 10-sites TFIM with hx = 1.The time series were sampled at a rate 1/∆t = 16/π with maximum time Tmax = 4π.

Figure 20 : 1 .
Figure 20: Dependence of the relative entropy error on the number of Trotter steps n Trotter for 10-sites TFIM with hx = 1.The time series were sampled at a rate 1/∆t = 16/π with maximum time Tmax = 4π.Note that n Trotter represents the number of Trotter steps per sampling period ∆t, so that the total number of Trotter steps for time t is n Trotter × (t/∆t).

Figure 21 :
Figure 21: Evolution of average squared magnetization for four different Monte Carlo chains of the 16-site TFIM with hx = 1 and β = βc using noisy Loschmidt echos via Aer-Simulator.Each colored point represents the average of the chain up to the specified iteration.The first 32 iterations (burn-in period) are excluded.The classical refers to the value of the equivalent Ising model (i.e., with hx = 0).

Figure 22 :
Figure 22: Average squared magnetization of the 16-site TFIM with hx = 1 at different inverse temperatures β.Estimates refer to the simulation results of the time-series Monte Carlo algorithm using either noisy Loschmidt echos (via AerSimulator) or exact ones (via QuSpin).Error bars represent 95% confidence (two standard deviations).The classical curve represents the values for the equivalent Ising model (i.e., with hx = 0).