Algorithms for quantum simulation at finite energies

We introduce two kinds of quantum algorithms to explore microcanonical and canonical properties of many-body systems. The first one is a hybrid quantum algorithm that, given an efficiently preparable state, computes expectation values in a finite energy interval around its mean energy. This algorithm is based on a filtering operator, similar to quantum phase estimation, which projects out energies outside the desired energy interval. However, instead of performing this operation on a physical state, it recovers the physical values by performing interferometric measurements without the need to prepare the filtered state. We show that the computational time scales polynomially with the number of qubits, the inverse of the prescribed variance, and the inverse error. In practice, the algorithm does not require the evolution for long times, but instead a significant number of measurements in order to obtain sensible results. Our second algorithm is a quantum-assisted Monte Carlo sampling method to compute other quantities which approach the expectation values for the microcanonical and canonical ensembles. Using classical Monte Carlo techniques and the quantum computer as a resource, this method circumvents the sign problem that is plaguing classical Quantum Monte Carlo simulations, as long as one can prepare states with suitable energies. All algorithms can be used with small quantum computers and analog quantum simulators, as long as they can perform the interferometric measurements. We also show that this last task can be greatly simplified at the expense of performing more measurements.


I. Introduction
The advent of quantum simulators [1,2] opens many exciting opportunities to probe and understand fundamental problems in physics, ranging from condensed matter to high energy physics and quantum chemistry [3][4][5].Feynman's original proposal in 1982 was to build a universal digital quantum computer that can imitate any physical systems.Although tremendous progress has been made, building a universal quantum computer that will fulfill Feynman's vision is still a long-term task.However, both near term (noisy) quantum computers and analog quantum simulators can already help us to address some of those problems.The latter, where the interaction is engineered directly according to the physical Hamiltonian under investigation, are particularly advanced in different platforms, like cold atoms in optical lattices [6,7], trapped ions [8], Rydberg atoms [9], quantum dots [10], superconductors [11], photons [12], etc.In particular, very controlled experiments can be carried out with around 50 qubits [13][14][15][16][17][18][19][20][21][22][23][24][25][26] and it is expected that this number will be significantly increased in the coming years.
There are many questions that crave answers from quantum simulators, especially physical properties of ground, non-equilibrium, and finite temperature states.Most of the theoretical work on quantum simulations has focused on the dynamics of many-body quantum systems, as well as on their properties at zero temperature.Since the first algorithm [2] that showed how the dynamics could be efficiently simulated, large improvements have been achieved leading to a very economic algorithm [27].In practice, in analog quantum simulators the dynamics are naturally implemented by letting the system evolve according to the engineered Hamiltonian [3].For ground state problems, the situation is quite different since determining its properties is very demanding and, in general, it requires exponential time in N , the number of qubits to be simulated [28].A quantum simulator can still be of big help since the corresponding classical simulator requires exponential resources both in time and memory, whereas the quantum one achieves a moderate speed-up albeit with polynomial memory.The first algorithms [29][30][31] used quantum phase estimation to project onto an eigenstate of the Hamiltonian, and have been successively improved [32][33][34][35].In particular, in [34] a cosine-filtering operator is used to prepare a state close to the ground state with a very small variance.It is similar in its conception to quantum phase estimation, but has a better scaling for that purpose.This idea has also been used in the context of tensor networks [36] to estimate the amount of entanglement required to achieve small energy variances along the whole spectrum of a Hamiltonian.Although all these quantum algorithms were originally designed for scalable quantum computers, proposals to use them with analog quantum simulators have recently been put forward [37].This can be very convenient for small simulators, as the exponential scaling of the resources still limits their applicability to large systems.Other heuristic algorithms, like adiabatic [38][39][40] and variational [41,42] strategies, can be very useful and overcome the exponential scaling in certain cases [4].
Quantum algorithms for excited states or finite temperature are more scarce.In [43] it is shown how to realize the imaginary time evolution operator to produce a Gibbs state, whereas other algorithms propose sampling techniques [44][45][46].Phase estimation can also be directly used to prepare states at different energies, and thus address quantum statistical questions in the microcanonical ensemble.All those algorithms may work well in practice.However, as for classical ones, they require an exponential time in N , although only polynomial memory resources.Additionally, it remains challenging to implement most of them with the existing small quantum computers or analog quantum simulators.
In this paper we introduce and test two different types of quantum algorithms to determine physical properties in an energy interval or at finite temperature (see Fig 1 for a graphical summary).The idea underlying our proposal relies on the cosine-filter of [34,36] to target states with small energy variance at selected energies, in which the observables could be measured.Preparing those states may nevertheless be challenging in practice.We overcome this obstacle by showing that observations obtained after running quantum simulators for different stroboscopic evolution times are sufficient to determine the values of interesting quantities, without the need to prepare the filtered state at all.This is in the spirit of the time series approach introduced in [47] and later used for the related question of estimating the (binned averages of) Hamiltonian eigenvalues [48] (see also extensions based on Hamiltonian querying and Chebyshev expansions [49,50]).Our algorithms explicitly target finite energy properties, such as microcanonical expectation values, for which we can demonstrate efficiency.The quantities required in our method can be obtained with interferometric measurements, which involve the conditional evolution depending on the state of a single qubit, and are specially suited for quantum simulators.More concretely: 1. Our first result is an efficient hybrid quantumclassical algorithm that targets the physical properties of states in an energy interval around that of any state that can be efficiently prepared (see Fig. 1 for an illustration).We prove that this algorithm can be carried out in time that is polynomial in N , the inverse error, and the inverse width of the filtering operations.The later is related to the energy variance of the targeted state.Up to our knowledge, there is no classical algorithm achieving this polynomial scaling.
2. Our second algorithm combines the quantum simulation with classical Monte Carlo methods, and provides practical methods to obtain both microcanonical and canonical expectation values of observables (see Fig. 1 for an illustration).These quantum assisted Monte Carlo methods use quantum simulators to compute the sampling probabilities that are required in Quantum Monte Carlo methods.Remarkably, they circumvent the sign problem [51], the main obstacle of applying such methods to many physics problems, so long as one can prepare (product or other kind of) states of suitable energies.
Let us briefly mention that our proposal is differ-ent from other emerging families of quantum methods that make use of hybrid schemes, such as variational quantum eigensolvers [41,42] and, in particular, hybrid ansatzes [52] or Quantum Subspace Diagonalization [53] (see also [54][55][56][57]).Different to our algorithms, such methods typically target ground state problems and, while some of them [53] also use a superposition of a state evolved to different times, the coefficients of the superposition are variational parameters that need to be optimized, unlike in our proposal, where the coefficients are fixed by the filter.An even more significant difference is that these methods are not proven to be efficient.
Finally, the interferometric methods required for the algorithms presented here have been used in the context of the Loschmidt echo [58] in NMR [59].They have also been proposed for ions [60], atoms [61] and, more recently, to perform phase estimation with such systems [37].We will give several alternative procedures to simplify that task, that can be applied in different situations.First, we will show that one can replace phase estimation by the ability of both preparing cat-like states [62] and having access to two additional internal states.This capability already exists in different platforms [26,[63][64][65][66][67][68][69][70][71].Then we will show that one can perform a similar procedure but without the requirement of the two additional states.Then, we will give a procedure that does not require the preparation of cat-like states, but proceeds with a sequence of measurements.Finally, we will give an even simpler method that works for Hamiltonians possessing certain symmetries, like XY or Hubbard models.The last two methods have an advantage with respect to the previous methods since no cat-like state needs to be evolved, and thus the method is more resilient against decoherence.However, they require a larger number of measurements.
The structure of this paper is as follows.In Section II we introduce the models and the basic idea of the cosinefilter.In Section III, we present the first algorithm and show how one can use a quantum simulator to efficiently compute certain expectation values around fixed energies.In Section IV we give more practical methods for the same purpose and we test them numerically.In Section V, we present a family of quantum assisted Monte Carlo algorithms for microcanonical and canonical observables that combine classical Monte Carlo and quantum simulators.We explore numerically their performance and demonstrate that they are robust against certain noises.In the appendices we present the methods to replace interferemotric measurements, describe the specific model we used in our numerical study, give details of the proof regarding the polynomial scaling of our algorithm, and investigate how much one has to decrease the variance in order to converge to the microcanonical and canonical results for a non-exactly solvable model.

Efficient hybrid quantum-classical algorithms for microcanonical observables (Sec. III-IV)
Given an easy to prepare state , an observable < l a t e x i t s h a 1 _ b a s e 6 4 = " j v g / m T 9 A and < l a t e x i t s h a 1 _ b a s e 6 4 = " e G 7 K T g L j M j 8 A i e S P A g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e G 7 K T g L j M j 8 A i e S P A g = = < / l a t e x i t > Classical CPU 1. Determine required times

NISQ devices
For each : 1. Prepare 2. Evolve 3. Measure j q i g z N q a y D c F b f n m V t C / q 3 l X d u 7 + s N a p F H C U 4 h S q c g w f X 0 I A 7 a E I L G K T w D K / w 5 m T O i / P u f C x a 1 5 x i 5 g T + w P n 8 A W S T k d E = < / l a t e x i t >

| i
< l a t e x i t s h a 1 _ b a s e 6 4 = " I G 5 v u y Y U 5 G l + f l u X 8 r s 5 3 A 0 a y A s j q i g z N q a y D c F b f n m V t C / q 3 l X d u 7 + s N a p F H C U 4 h S q c g w f X 0 I A 7 a E I L G K T w D K / w 5 m T O i / P u f C x a 1 5 x i 5 g T + w P n 8 A W S T k d E = < / l a t e x i t >

| i
< l a t e x i t s h a 1 _ b a s e 6 4 = " e C g E K y q P T i 0 Z o U K g 5 a l 1 t W W G 9 W w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L C 2 C p 5 K I V I 8 F L x 4 r 2 g 9 o Q 9 l s N + 3 S z S b s T o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 < l a t e x i t s h a 1 _ b a s e 6 4 = " e C g E K y q P T i 0 Z o U K g 5 a l 1 t W W G 9 W w = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L C 2 C p 5 K I V I 8 F L x 4 r 2 g 9 o Q 9 l s N + 3 S z S b s T o Q S + h O 8 e F D E q 7 / I m / / G b Z u D t j 4 Y e L w 3 w 8 y 8

NISQ devices
< l a t e x i t s h a 1 _ b a s e 6 4 = " s a h l / f 6 y q L G u

Heuristic quantum-assisted Monte Carlo algorithms for microcanonical and canonical observables (Sec. V)
For inverse temperature , an observable A and < l a t e x i t s h a 1 _ b a s e 6 4 = " e G 7 K T g L j M j 8 A i e S P A g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e G 7 K T g L j M j 8 1. Graphical summary of the concept and main results of the paper.Left: Our algorithms compute properties of a state with narrow energy variance (schematical spectral distribution shown in orange), which would result from applying a filter onto an easy to prepare state, e.g. a product state (green curve).For local Hamiltonians, for which the density of states in the thermodynamic limit approaches a Gaussian distribution of width proportional to √ N (blue curve), the accessible energy densities can lie on the tails of the spectrum.The box under the graph shows the actual energy spectrum for a N = 10 Ising chain.Right: The main contribution of our paper is the proposal of two sets of algorithms, schematically summarized in this figure.The first one (above) is a provably efficient hybrid quantum-classical algorithm for computing expectation values of observables in filtered states as illustrated on the left.It is based on repeated preparation, evolution and measurements on an easy to prepare state, run by a quantum device, and postprocessing by classical computation.The second method is a set of quantum assisted sampling algorithms that can be used to compute microcanonical and canonical properties (summarized below for the canonical case).They perform a classical Monte Carlo importance sampling where the quantum device is used to determine the sampling probabilities efficiently.They require shorter coherence time than the first algorithm, at the price of more measurements.

A. Setup
We consider N spins on a lattice and a Hamiltonian and denote by E min and E max the minimum and maximum eigenvalues of H.We will assume that The Hamiltonian H could be local in any spatial dimension, i.e. the term h n only acts on the n-th lattice site and its neighbors.However, we emphasize that the algorithms proposed here can also be applied to more general setups, where h n has long-range interactions or even with more complicated Hamiltonians, although some of our estimations will rely on the original form (1) or in its locality.The main requirement is that the evolution generated by H can be efficiently implemented with the quantum simulator.More specifically, that given an initial state, ψ, and an observable A, one can efficiently determine with a sufficiently small error.Furthermore, if one can measure those quantities, one also has access to . The value of (3) can be readily measured using a quantum computer (see, e.g., Ref. [72]).For analog quantum simulators, this may not be possible.In Appendix A we give a series of alternatives to obtain such quantities.In general, by repeating the experiment L times, the error will be additive and scale as L −1/2 .

B. Initial states
In order for the simulation algorithm to be efficient, we will choose as ψ a state that can be prepared by the simulator.The simplest are product states where p n are normalized states, e.g., for qubits, In general, we denote the mean energy and variance of H in the state ψ by In case H is local and the state ψ has finite correlation length, both E ψ and σ 2 ψ will scale as N .The energies of states that can be efficiently prepared will determine the range of energies that our algorithm can efficiently explore.If we restrict ourselves to product states, the mean energy E p does not cover the whole spectrum of H; there are energies E p,min and E p,max such that we can always choose a state p with E p in the interval [E p,min , E p,max ] but never outside (see Fig. 1 for an illustration).It has been shown that the range of the interval can be extensive in N for local Hamiltonian (See Appendix C for more details).Finding product states within this interval amounts to solving a mean-field problem, thus can be done efficiently on classical computers.In order to access energies outside this interval, one can consider other states, ψ, that are still easy to prepare but can cover a wider range of energies.In particular, we could consider products of spin blocks, matrix product states [73][74][75], or states obtained through adiabatic evolution or variational methods.Alternatively, the state ψ could be prepared by starting from a product state and running the quantum simulator with a different Hamiltonian for some time.

C. Cosine Filter
Following [34], we define the cosine-filtering operator where we use ... 2 to indicate the nearest even integer.
Here δ can take arbitrary values, including decreasing (eg δ ∼ 1/N ) or constant (δ ∼ 1) with N , that will be considered in the following sections.In order to interpret the action of this operator, it is useful to approximate [34] as long as the spectrum of the operator that appears in the argument of the cosine lies in the interval [−π/2, π/2] (in fact, this is also true in a bigger interval, see Appendix E).Thus, it basically projects out the eigenstates of H that have an energy E with |E − E| δ, and thus acts as a filter around E [76].By definition, 0 < P δ (E) ≤ 1 1.
As in [34], we approximate up to an error (in operator norm) bounded by 2e −x 2 /2 for X ∞ ≤ 1, and where We can use this expansion to express the cosine-filter (7) in terms of the evolution operator e −iHt for certain times t.For |E| ≤ N/2, we take X = (H − E)/N so that where The idea will be, as in [36], to apply (9) to certain states in order to filter them around some energy E, and then to obtain expectation values of observables with the resulting states.However, instead of preparing the state, we first express the desired values in terms of (2), and then use the quantum simulator to measure those values separately.The number of measurements will be 2R times the number of repetitions required to obtain a prescribed accuracy.Each of the runs of the simulator will be for a time t ≤ 2x/δ.In the end, we perform the multiplications and sum classically.

D. Ising Model
In order to benchmark our algorithms, we will use a rather trivial model for which we can obtain numerical results for values of N ∼ 100 qubits, which should be attainable in present or planned quantum simulators.Let us take an even number, N , of fermionic modes and a Hamiltonian where a n are annihilation operators of the vacuum |vac and we have chosen periodic boundary conditions for the fermions: a N +1 = a 1 .Through the Jordan-Wigner transformation, this corresponds to the Ising Hamiltonian where σ x,z are Pauli operators, with appropriate boundary conditions.Defining the operators in momentum representation where k = −N/2 + 1, . . ., N/2, we will perform some computations with (Fock) states of the form where k 1 < k 2 . . .< k .They form an orthonormal basis and, even though they are not eigenstates of H, they are easy to deal with for large values of N .In appendix B we provide the analytical formulas that allow us to obtain exact numerical results for this model.The minimal eigenvalue of H, E min , and the minimal energy attained by a state of the form ( 16), E p,min can be easily computed with those formulas.Even though we are interested in energies above E min and finite temperatures, let us notice that at zero temperature the Hamiltonian (13) features a phase transition at g = h.For h g the ground state is the vacuum, whereas for h g it is a superposition of the vacuum with states where excitations occur in pairs of momenta ±k.

III. Efficient quantum algorithm for observables at finite energy
Given a state ψ and an observable A, we define Both quantities are related to the microcanonical expectation value of A. In particular, if E|ψ = 0, where |E is the eigenstate of H corresponding to the energy E, (17b) converges to that value in the limit δ → 0.
The numerator and denominator of ( 17) can be expressed in terms of a A,ψ (t n , t m ).Thus, the quantum algorithm uses the quantum simulator to determine those quantities up to the required precision, computes classically c m and exp(iEt m ), and then performs (classically) the required sums and multiplications.We will show that both (17a,17b) can be efficiently computed using a quantum simulator that has access to (2).But for that, we have first to explain what we mean by "efficiently" and also formulate the problem more precisely.
We say that a state ψ can be efficiently prepared if, for any prescribed error, > 0, we can obtain a state ϕ with ϕ − ψ 2 < in a time Furthermore, we say that the quantum simulator can efficiently measure A if it can perform measurements to obtain a A,ψ (t , t) (with t ≤ t) with an error smaller than in a time Note that this basically requires an efficient procedure to evolve according to the Hamiltonian and the possibility of performing interferometric measurements.
Result: If a quantum simulator can efficiently prepare ψ and measure A, with A ∞ ≤ 1, then for any , δ > 0, one can always find , so that one can obtain (17) up to an additive error in a time including the cost of finding the value E.
Note that for δ = poly(1/N ), the result can still be obtained in polynomial time.Furthermore, since σ ψ ≤ N/2, E will differ from E ψ by at most a constant times σ ψ log 1/2 (N ).If σ ψ ∝ O( √ N ) as it occurs in states with finite correlation length and local Hamiltonians, this difference will only scale as √ N log N .The reader may wonder why we need to introduce an interval, instead of simply fixing E to some value, for instance E = E ψ .The reason is that the spectrum of H is discrete and it may well be that ψ|P δ (E)|ψ is exponentially small in N .As we will show in Appendix C, this issue can be avoided if we are allowed to vary E in a small interval.We are not aware of any classical algorithm that can achieve this scaling.
The result can be proven by expressing the numerator and denominator of (17a,17b) as a function of (2) by means of (9), and then showing that both, as well as their quotient, can be computed with the required accuracy in a time (21).Here, we show it explicitly only for (17a), but it can be done in the same way for (17b).We define the notation and denote by ∆p, ∆q the bounds on their error; that is, if the measured values are p, q, they fulfill |p − p| < ∆p (and analogously for q).Let us first argue that if we require ∆p to scale polynomially with N −1 , δ and , we can reach it with the quantum simulator in a time (21).The reason is that, using ( 9) and ( 12) we will only need to determine the 2xN/δ values a A,ψ (t m ), and each of them will require to run the quantum simulator for a time ≤ 2x/δ.Furthermore, we will have to repeat the procedure a number of times (21) in order to reduce the error.The time to perform each of those tasks also scales in the same way, given that, by assumption, the quantum simulator can efficiently measure A. Thus, we conclude (21).An analogous argument applies to ∆q.For a given E, the total error |p/q − p/q|, will be upper bounded by ∆p + p∆q/q q − ∆q ≤ ∆x + 2∆q/q q − ∆q , as long as ∆q < q, and where we have used that p ≤ 2.
One can readily check that if ∆p = q/3 and ∆q = q 2 /6, then the error will be bounded by .Thus, the problem is reduced to proving that q = poly(N, 1/δ), (24) since in this case, ∆p and ∆q will scale polynomially with N −1 , δ and which, as we argued, can be accomplished with (21).As for the cost of finding the value E, we show in Appendix C that there always exists an interval of energies of size ∆E ≥ δ 2 /6N within the interval (20) where . Since σ ψ ≤ N/2, we have that in that interval q fulfills (24).Thus, the procedure consists in dividing the interval (20) in 24N rσ ψ /δ 2 equal slices and picking an energy E in each of them.At least one of them is then guaranteed to fulfill (24) and thus we will be able to determine (17) with an error smaller than in a time (21).

IV. Practical computations
In practice, one can use the quantum simulator much more efficiently than what has been presented in the previous section, and also employ it to access other physical properties.In this section we propose and analyze several algorithms to compute different quantities related to the microcanonical and canonical quantum statistical ensembles.We will also illustrate them with some examples for the model of Section II D.

A. Local density of states
The simplest quantity is a broadened version of the local density of states, which (up to a factor) converges to that quantity in the limit δ → 0. Using (9), we can express with (12).As before, our algorithm uses the quantum simulator to determine a ψ (t m ).The method can be made more efficient by noticing that norm of X 0 = (H − E ψ )/(rσ ψ ) will be bounded by one in the subspace where the state ψ has most of its weight if we choose r ∼ 1.Thus, we can use the expansion (9) with X = X 0 and M = r2 σ 2 ψ /δ 2 to obtain (25), but now with R = xrσ ψ /δ and t m = 2m/(rσ ψ ).For σ ψ ≤ √ N the number of required measurements will significantly decrease with respect to (12).In that case, we denote rσ ψ = r √ N so that In Fig. 2 we have plotted D δ,k (E k ) in logarithmic scale for N = 100 spins, δ = 0.1, and the Hamiltonian (13) with g = 1, h = 2, for 50 randomly generated states |k (16).We have taken x = 3 and compared the results of the original (12) (circles) and the optimized alternatives (27) (crosses for r = 0.4, plus symbols for r = 1).For this value of δ, the improved method yielding (27) with r = 0.4 corresponds to 120 measurements of a ψ (t) and a maximum value of t = 60.For δ = 1 it requires 12 measurements with a maximum value of t = 6, which is very reasonable for present experiments.We observe that for r = 1 one already obtains an error of the order of 10 −3 , whereas for r = 0.4 it is about 10 −2 , which is what one could expect with imperfect devices.

B. Microcanonical observables
The very same simplified program can be applied to (17), as we can also choose different values of r to make the procedure much more efficient.In Fig. 3 we have plotted (17a) for the model ( 13) with g = 1, h = 2, and we have chosen the energy as an observable, i.e.A = H/N .Note that the microcanonical expectation value of A = H at E is only trivial (i.e., equals to E exactly) in the limit δ → 0. Here, our purpose is to investigate the convergence of H δ,ψ (E) with respect to δ.We have again chosen 50 random states |k , and subtracted their mean energy E k /N to optimize the visualization, since We have taken δ = 0.1 for the triangles, while δ = 1 for plus symbols and circles.Only for the latter we have chosen the more efficient version (27), with r = 1.We observe a clear difference (at the percent level) for different values of δ, however the value r = 1 is sufficient to obtain reliable results.An interesting question in this context is to what extent one can recover the microcanonical expectation value by decreasing δ.This makes sense if the system is sufficiently large and fulfills the eigenstate thermalization hypothesis (ETH) [77,78], which ensures the convergence of the procedure in the thermodynamic limit.The question of how narrow the energy support of a pure state needs to be in order to recover microcanonical expectation values has been analyzed in Ref. [36,79] for one-dimensional systems.In [79] it is concluded that, for generic systems, δ needs to decrease with N , while [36] gives evidence that for local observables it may suffice that δ ∼ 1/ log(N ).In Appendix D we analyze this question for a non-integrable model and up to N = 28 spins, and give further evidence for the need to decrease δ with N .This can be qualitatively understood as follows.Let us denote by |E the eigenvectors of H with energy E. In case the ETH applies, the diagonal matrix elements of physical observables (which can also include, for instance, correlation functions.) in the energy basis rapidly converge to the microcanonical expectation value, while off-diagonal elements vanish exponentially fast [77].However, when one has a superposition of an exponential number of eigenstates around some energy, the sum of the off-diagonal terms does not need to converge even though the diagonal terms do.This is also the reason why δ must decrease with N , since P δ (E)|ψ contains superpositions of |E and thus the expectation value depends on off-diagonal elements.

V. Quantum assisted Monte Carlo algorithms for microcanonical and canonical observables
The considerations at the end of the previous section suggest an alternative strategy that can be used for microcanonical and canonical observables.At a high level, the algorithms we propose in this section can be seen as a quantum version of the classical quantum Monte Carlo algorithms: they use classical Monte Carlo to sample different initial states, while the quantum device assists with the computation of sampling probabilities and the measurements of observables.
Different to the algorithm presented in Sec.III, the quantum assisted sampling discussed here is not proven to be efficient.However, when compared to the first algorithm from an experimental point of view, these methods offer the potential advantage of requiring shorter coherence times, at the cost of increased number of measurements in the quantum device.Moreover, the sampling probabilities are always positive and the sign problem in classical quantum Monte Carlo [51] is circumvented.We also present numerical experiments for a N = 100 Ising Hamiltonian, to demonstrate the viability of these methods on near-term devices.The results show that physical quantities can be obtained accurately even in the presence of certain level of noise.

A. Microcanonical observables
The quantity may converge to the microcanonical expectation value even for constant δ since, by definition, P δ (E) is diagonal in the eigenbasis of H and thus no off-diagonal element appears in the expectation value.This is indeed observed in Ref. [80], where tensor network techniques are used in order to compute quantities closely related to (28) with energy resolution corresponding to larger values of δ than those required by (17).In fact, that a constant δ suffices also follows from the fact that, if the thermodynamic limit N → ∞ exists for the observable A, then in that limit the expectation value of A in almost all eigenstates of H should coincide if the corresponding energies fulfill The intuitive reason is that in that limit, only intensive quantities matter.Therefore, what is relevant is δ/N , so that as long as it vanishes, one should obtain the thermodynamic value.
In principle, one could obtain (28) in the very same way as in (17a).This can be seen by noticing that where Φ is a maximally entangled state of each spin with an auxiliary one That is, just one has to add an auxiliary qubit for each existing one and prepare an entangled state of each pair.Thus, from (29,17a) it follows that one could compute the numerator and denominator independently with the help of the quantum simulator, and then the quotient.However, we face here the problem that the denominator will typically decrease exponentially with N .For local Hamiltonians, the reason is that σ Φ ∝ √ N and therefore, for any extensive value of the energy, E = eN , Φ|[P δ (E) ⊗ 1 1]|Φ ∼ exp (−cN ) for some c = O(1).This makes this procedure impracticable already for N > ∼ 20.In the following, we propose an algorithm to circumvent, at least in part, this issue.Let us re-express (28) with the help of an (over)complete basis of states fulfilling where dµ ψ is a measure in the basis set.For instance, we can take an orthonormal basis of product states |p 1 , p 2 , . . ., p N or all product states (5), in which case Inserting ( 31) in ( 28) and using definitions (17a) and ( 28), we can rewrite This quantity can be computed using Monte Carlo algorithms so long as one is able to compute D δ,ψ (E) and A δ,ψ (E).But these two can be computed using the quantum simulator.More concretely, one can implement importance sampling according to a distribution proportional to D δ,ψ (E).For instance, using a basis of product states, as mentioned above, this can be achieved by a Metropolis-Hastings algorithm in which once a move is proposed, the quantum simulator is used to estimate the corresponding value of D δ,ψ (E), and thus determine the acceptance rate.The Monte Carlo algorithm circumvents not only the burden of summing over all states, but also the need to measure D δ,ψ (E) with an exponentially small accuracy, because D δ,ψ (E) needs to be evaluated only for states ψ for which it is not negligible.We also emphasize that if the observable A is chosen so that A|ψ = λ|ψ , and λ can be classically computed, then one only has to determine D δ,ψ (E) with the quantum simulator.
We have tested this algorithm with the Hamiltonian (13) and the simplest Monte Carlo algorithm that changes one spin at a time for the sampling.Notice that our goal is to demonstrate the viability of the approach, rather than the competitiveness of the algorithm, since that would require optimizing the sampling methods and other parameters.We have chosen the Hamiltonian so that we can compare the results with an exact calculation for N ∼ 100 spins.We have implemented a Metropolis algorithm that takes a random state |k (16), according to a probability proportional to D δ,k (E).We have then computed the magnetization A = M with as a function of E. Since the model is exactly solvable, we have also computed A δ (E) directly using the method of Appendix B. This direct numerical calculation requires the computation of very large and small numbers, so that one can easily run into precision problems.In fact, for some plots we can only provide the exact result for some values of E, since otherwise our exact method did not give consistent values.In Fig. 4(a) we plot M δ (E) for N = 20 spins, g = 1, h = 2 and δ = 1, 4 (red and blue lines) obtained with the numerical computation.The symbols are obtained with the Monte Carlo method with δ = 1 (circles) and δ = 4 (squares).We have checked that the results for δ = 1 and δ = 0.1 are almost indistinguishable, and this is why we plot only the first ones.We have sampled 10 5 times for each point.In Fig. 4(b) we plot the same for g = 1 and h = 2, also showing very good results.Notice that the lower curve terminates at around E k ∼ −25; the reason is that the lowest energy that can be reached with the states ( 16) is E p,min = −20, −14.39 for Fig. 4(a,b) respectively, and thus, at lower energies, D δ,ψ (E) becomes very small for any state ψ.Thus, not only the Monte Carlo method but also the exact one encounters convergence problems for such low energies and this is why they are not plotted.We have used an exact summation to obtain the dotted line.This is only possible because we only have N = 20 spins in this figure.
In Fig. 5 we have considered larger systems, with N = 100.In Fig. 5(a) we depict the magnetization for g = 1 and h = 2.For a more convenient graphical representation of the data, we plot in the inset the modified quantity M δ (E) = M δ (E) − 1/2 − 0.004E/N as a function of E. The solid lines indicate the exact results for δ = 1, 4 (red and blue), while the symbols show the Monte Carlo results (circles and crosses for δ = 1, squares and plus symbols for δ = 4).As before, we have sampled 10 5 times for each point.The cross and plus symbols in the inset indicate the results when a cutoff of 10 −2 is set for D δ,k (E).That is, in the Metropolis method, as soon as we compute it and obtain a smaller value, we set it to zero.In practice, this intends to resemble an experiment where this quantity has been obtained to that precision.The solid lines are exact results, which we have only computed up to some values of E, since otherwise we encountered precision problems.As one can gather from the plots, the Monte Carlo results resemble well the corresponding values.For E < −60 the exact method encounters precision problems, while this happens for the Monte Carlo method only for E < −100 = E p,min .When decreasing δ even further, the results are practically indistinguishable from those of δ = 1.However, the program gets unstable for energies E < ∼ E p,min since, as expected, the values of D δ,k (E) become extremely small due to the Gaussian dependence.

B. Canonical observable:
Now we show how, using similar ideas, one can also compute canonical observables, i.e., where β is the inverse temperature.We use the fact that for sufficiently small δ, we can approximate where E rhs)0,1 = E 0,1 ∓yδ with E 0 (E 1 ) the lowest (highest) eigenvalue of H, and y 1.This motivates the definition which converges to the canonical value for δ → 0. Using (33), we have The quantum algorithm proceeds in the same way as before, by using the quantum simulator to recover a ψ (t m ) and a A,ψ (t m ), and then a classical computation to do the rest.In particular, Monte Carlo samples both the states ψ and the energies E with a probability proportional to e −βE D δ,ψ (E).For observables that are diagonal in the chosen basis, it is also possible to sample only over states, with the integrated probability dEe −βE D δ,ψ (E), which can be reconstructed using the quantum simulator in the same way as for the microcanonical observables, since the energy dependence can be integrated analytically.
We have performed Monte Carlo computations and shown the results in Fig. 6.We have plotted the magnetization (34) as a function of the inverse temperature β for the Hamiltonian (13) with N = 100 spins and g = 0.3, h = 0.8 (lower curve) and g = 0.4 and h = 0.4 (upper curve).The solid line is the exact result (B13) computed with the formulas given in Appendix B. Note that here there is no problem with the precision, as most of the products in the numerator and denominator cancel and one ends up with a simple sum (B14), and that the result is independent of δ [since the definition (B13) is, too].The symbols are obtained with the Monte Carlo simulation for δ = 1 and x = 3.We have discretized the integral in energy appearing in (38) by taking E from −N to N in steps of 0.5, and sampling only those values, although the Metropolis algorithm only took values around certain energies, E±5, for each value of β.For the Monte Carlo methods we took E 0,1 = ±3N/2, but those values were never reached in the Metropolis sampling.From the figure we conclude that the result converges very well already for δ = 1, and also that the performance of the Monte Carlo method is very good.One can observe that for larger values of β, there is a little bias towards smaller values of the magnetization.In the inset of Fig. 5, we also display the results carried out with the Monte Carlo computation with a cutoff 10 −2 , still rendering competitive results.We note that the fact that the upper curve corresponds to the critical point g = h does not have relevant consequences at the temperatures considered here.For lower temperatures, the program does not give reliable results as it has to scan energies for which D δ,ψ (E) becomes very small.

VI. Summary and Outlook
We have proposed and analyzed two types of quantum algorithms to characterize quantum many-body states in finite energies intervals and temperatures.They are based on the cosine-filter, which when applied to a state, reduces its variance to a predetermined value around a given energy.However, instead of preparing the filtered state, our algorithms compute different quantities that allow one to reconstruct expectation values of observables or other magnitudes.The algorithms we proposed are quite flexible -they can either be implemented on digital quantum computers (e.g., based on superconducting qubits, trapped ions, or Rydberg atoms) and on analog quantum simulators (e.g., based on cold atoms, trapped ions, photons), as long as they can perform certain kinds of interferometric measurements.We have also presented in Appendix A more practical approaches which, in some cases, possess clear advantages with respect to such measurements.
We have shown how our first algorithm can be used to to efficiently compute expectation values over filtered states, as long as we can prepare the initial state.The algorithm requires a time that scales polynomially with the system size, N , the inverse variance, and the inverse error.We have also shown how it can be simplified in practice, leading to a practicable method for present or planned quantum simulators.The simulator has to be run for times ≤ 6/δ (where we took x ∼ 3 and r ∼ 1) and therefore for δ = O(1) this should be feasible for existing devices.The price one has to pay is that one has to perform many more measurements.This number can be estimated by taking into account that we need to compute sums of the form (26), where we add 2R terms so that to have a total error of the order of , this will require of the order of R −2 ≈ 3 √ N /δ 2 measurements.For N = 100, δ = 1 and = 10 −2 this yields of the order of 3 • 10 5 measurements.Actually, by taking into account that c m becomes small, using smaller x ∼ 1, exploiting symmetries and other optimizations, we expect that it is possible to reduce significantly this number, to the order of 10 4 .
The second algorithm we introduce proposes a way to combine the first algorithm with Monte Carlo simulations in order to obtain the expectation value of microcanonical and canonical observables.The resulting methods also rely on the ability to prepare certain states (e.g.product states) efficiently, to run the quantum simulator for times of the order of 6/δ and to perform interferometric measurements.Our numerical simulations of these algorithms for a simple model indicate that with several tens of thousands of samplings one can obtain reliable results, at least in some energy regimes.Those algorithms can also be significantly improved by using standard Monte Carlo strategies to speed up convergence.These quan-tum assisted sampling algorithms open new possibilities to study thermal properties of quantum many-body systems with near-term quantum devices.
There are other modifications that may help to improve the algorithms.First, one could take a different expansion of the filter.For instance, one can use instead of the cosine-filter, Chebyshev expansions for a quantum computer [34,36,49,50], a low-pass filter, or choose c m differently for an analog simulator to better adapt to the specific errors in which it may incur.Second, in the sampling, there is information that can be collected to investigate other physical questions.For instance, the spins or energy samples that are used in the Monte Carlo methods can be used to estimate other quantities directly (like higher moments), without the need to make other computations.Another possibility is to use the adiabatic or variational algorithms to prepare states with small variance to start with, which would make the algorithms more efficient, or allow us to access energies which are out of the scope of product states.
We have formulated most of our results for a manybody spin Hamiltonian with local interactions.However, they can be equally applied to fermionic systems, systems with longer-range interactions, or disordered setups, as long as the quantum device can emulate the corresponding Hamiltonians.Also, the ideas developed here can be easily adapted to measure other quantities of interest, like Green functions or structure factors.
A. Ways to retrieve a A,ψ (t) In this appendix, we describe different methods to retrieve a ψ (t) and a A,ψ (t) from analog quantum simulators.Those typically possess severe limitations in the available operations as compared with quantum computers, and thus they may not be able to carry out the interferometric measurements that are needed to obtain those quantities.The methods adapt to different situations and have different requirements.
We start out by briefly reviewing the standard method based on conditional dynamics where, during the evolution, all the qubits have to be coupled to an extra one, called the control qubit.Then, we analyze a closely related method where the interaction with the control qubit only needs to occur at the beginning and at the end of the evolution, although at the expense of using other internal states.The third method does not require a control qubit but the possibility of creating certain cat-like states, as well as an extra internal state; the latter can be avoided for certain kind of Hamiltonians for which one can prepare one eigenstate efficiently, like, for instance, Hamiltonians of Heisenberg type.The fourth method builds on the previous one and does not need cat states.This is thus much simpler to implement in practice and may be more robust against decoherence.However, it requires more measurements.The last method applies to Hamiltonians with some special symmetries, like XY or Hubbard models, and can be very practical.

Conditional dynamics
The conceptually simplest way to obtain a A,ψ (t) is to carry out conditional dynamics [72] depending on the state of one of the qubits, the control qubit, which we call c.This corresponds to the operation If we denote by H c the Hadamard transformation on the control qubit, then one can first produce the state By then measuring the control qubit in the computational basis, and the observable A on the rest, one can retrieve (2).This method requires the ability to couple all the qubits to the control one during the evolution in order to apply (A1), which may be difficult in practice.However, it can be very naturally applied in trapped ion simulators [37,60], since all ions are coupled to the same phonon bus, and in Rydberg atoms in optical lattices, as one atom in a Rydberg state can influence the dynamics of the rest [81].In fact, very efficient techniques have been proposed to perform this kind of dynamics and measurements using that implementation [82][83][84].
which is the vacuum of spin up, and an arbitrary Fock state of spin down modes, and define Then, an orthogonal basis of common eigenvectors |ψ can be formed as all possible products of one of these factors for each pair of sites.

B. Ising model
In this appendix we give some formulas that we have used in our numerical illustrations regarding the Ising model.Let us consider N fermionic modes and a Hamiltonian (B1) where a n are annihilation operators of the vacuum |vac and we have taken periodic boundary conditions for the fermions, a N +1 = a 1 .Through the Jordan-Wigner transformation, this corresponds to the Ising-like Hamiltonian where σ x,z are Pauli operators, with appropriate boundary conditions.
As usual, we first perform a Fourier transform with (we assume even N ) The new Hamiltonian is where , and Note that tr( Hk ) = 0.
We now proceed as follows: for each value of k = 1, . . ., N/2 − 1 we write H k = Hk as a 4 × 4 matrix in the basis and also define H 0 = H0 + HN/2 and write it as a 4 × 4 matrix in the basis We obtain where and σ z = −1 0 0 1 .The direct sum structure corresponds to the subspaces {|1 , |2 } and {|3 , |4 }.Thus, the problem is reduced to N/2 non-interacting 4-level systems.
Although we are interested here in finite energies and temperatures, we now briefly discuss the zero temperature behavior (i.e., the ground state).In the basis introduced here, it can be written as with eigenvalue E 0 = − k z k .The coefficients α k and β k depend on g and h.In particular, for h g we have α k ∼ 1 and β k ∼ 0, whereas for h g they change with k, and thus one gets superpositions of many configurations.At g = h, z N/2−1 → 0 as N → ∞, so that the gap closes and there is a quantum phase transition.

Microcanonical average
We compute now some of the quantities that are used in the numerical illustrations.For the sake of simplicity, let us assume that where A k only depends on b ±k for k = 0 and on b 0,N/2 for k = 0. Mapping back to qubits, such operator A is that can be decomposed into a sum of local operators, of which the magnetization we used in Eq. ( 34) is an example.
Let us start with where c m are defined with respect to δ and R is given in (12) or (27).We have which can be readily computed using (B9).We can also compute (28) directly, where tr(e −i2mH k /N ). Defining we can write Using (B9), (B9b) we have while the values of s m,k depend on the observable.

Canonical average
In a similar way we may compute the canonical average

C. Efficient computation
Here we prove the statement that was used to show that one can compute (17a) efficiently.Namely, we show that for any state ψ, and for any δ ≤ N/ √ 2, there exists an interval of width at least δ 2 /6N contained in an energy window |E − E ψ | ≤ rσ ψ (20) such that, for any E in that interval, and We emphasize that for our purposes we do not need a tight bound, so that we will be very rough when bounding different quantities with the goal of obtaining simple expressions.Moreover, although n(E) involves the cosine filter (7), it is simpler to bound the result of the Gaussian approximation (8).Specifically, we prove that there is an energy interval of radius at least δ2 /3N , contained in the window (20) and with its center inside the spectral limits of H, such that for any energy inside this interval, we will have that n(E) ≥ ñ(E).Because H ∞ = N/2, and the center of the interval is within the limits of the spectrum of H, we can ensure that H − E ≤ N + δ 2 /6N ≤ 13N/12, and thus (C4) is satisfied for E, so that bounding the Gaussian approximation is enough for the desired result.
Let us denote by First, we will show that for (C2), T (r) is upper and lower bounded by some quantity, which will indicate that there exists some E in the interval such that ñ(E) is sufficiently large.Then, by upper bounding the derivative of ñ(E), we will conclude that there is a neighbourhood of that E where ñ(E) fulfills the required condition.We write T (r) = T 0 − T 1 (r), where and T 1 (r) = T (r) − T 0 .Given that ñ(E) ≤ 1, we have where the last step uses that for (C2), r ≤ 2/π.We can perform the integration in (C6) explicitly using (C3) with (C2), to get Thus, putting things together we obtain the lower bound ) where we have chosen (C2) and used that r 2 /2 > log[2(1 + σ 2 ψ / δ2 ) 3/2 ].In order to bound T (r) from above, let us denote by E 0 the value where ñ(E) attains its maximum within the interval (20).Note that E 0 must be within [E min , E max ] (since, if it was outside, choosing the closest extreme eigenvalue of H would yield a larger value for ñ).We then have It is now possible to show that in a neighbourhood of E 0 of radius δ 2 /3N , ñ takes sufficiently large values, ñ(E) ≥ ñ(E 0 )/2.For that, we just have to bound the derivative of ñ within the interval (20).Let us call ñ max the maximal value of the derivative of ñ(E) with respect to E in that interval, which occurs at some E = E 1 .Although E 1 could lie outside the limits of the spectrum of H, if that is the case, and using that x e −x 2 /2 has a maximum at x = 1, it is easy to see that the maximum cannot occur more than δ away from the edges, so that H − E 1 ∞ ≤ N + δ ≤ 3N/2.Then we can bound, Let us finish this appendix with a remark on the difficulty of finding the state ψ given a prescribed energy.As we mentioned at the end of Sec.II B, this can be analyzed, for instance, through mean field theory for product states, or with matrix product states in one-dimensional problems.Regarding rigorous bounds, for product states it was shown by Lieb in [86] that for any local Hamiltonian with spectrum contained in [−N, N ], one can efficiently find a product state with an energy −DN , where D > 1/9.While this is a theoretical bound, we expect that for most relevant Hamiltonians D will be much larger (in fact, one can find much tighter bounds for specific models, like those used in this paper).Furthermore, as emphasized in the main text, one does not necessarily have to use product states, which gives access to even larger values of D.

D. Convergence to the microcanonical and canonical values for a non-integrable model
In this appendix we numerically investigate the convergence to the microcanonical and canonical values of the quantities defined in the main text using exact diagonalization.In particular, we show how a polynomially decreasing δ ∼ poly(1/N ) seems to be enough for the quantity A δ,ψ (E) (17b) to converge to the true microcanonical expectation value.
We consider the Ising model in a tilted field, described by the Hamiltonian which is in general non-integrable, except in the limits g = 0 (classical) and h = 0 (transverse field Ising model).In the following, we choose a strongly nonintegrable point h = 0.5, g = −1.05[87], and we have taken J = 1.Notice that this corresponds to a different normalization for H/N than the one used in the main text.However, it is enough to ensure that (7) filters out energies much farther than the width δ for all the states analyzed here (see appendix E), so that it allows us to study how the microcanonical values are approached as the width decreases.We consider real translationally invariant product states, which can be parametrized as |Ψ(θ) = |p(θ) ⊗N , where |p(θ) = cos θ|0 + sin θ|1 .In the thermodynamic limit, these states have energy density E/N = cos 2 (2θ) + h cos(2θ) + g sin(2θ), ranging over most of the energy band (see Fig. 8).
We choose three values of θ corresponding to states in the lower part of the energy band, namely θ 1 = π/4, for which E 1 /(JN ) = g = −1.05,θ 2 = π/3, for which E 2 /(JN ) = −0.909,and θ 3 = π/6, with E 3 /(JN ) = −0.409.For each of these states, and for several local observables, we compute exactly the expression A δ,ψ (E = E ψ ) from Eq. (17b) for different values of δ.In order to check that it is enough to decrease δ polynomially with the system size, we run the calculations for system sizes 10 ≤ N ≤ 28, and δ ∝ N s , for s = 0, −1, −2.The results, for observables A = σ [N/2],z ⊗σ [N/2+1],z and A = σ [N/2],x , are shown in Fig. 9.As reference, we estimate the microcanonical expectation values in the thermodynamic limit using uniform MPS [74] (more concretely, we approximate the canonical ensemble at the same energy density in the thermodynamic limit as a matrix product operator, in which the observables can be easily computed, and use the fact that in this limit, both ensembles are equivalent).
Our results indicate that, although a constant value of δ is not enough for A δ,ψ to approximate the microcanonical value, when δ decreases as 1/N 2 , the expectation values indeed converge.For δ ∝ 1/N we observe that the values are reasonably close, and they would be compatible with a slower convergence.

E. Approximating the Cosine filter as Gaussian
We want to bound the absolute value of the difference f M (x) = e −M x 2 /2 − cos M x. (E1) Since both terms are even, we can consider only x ≥ 0. If |x| < π/2, e −x 2 /2 ≥ cos x, and both terms are positive, so also 0 ≤ cos M x ≤ e −M x 2 /2 , and The bound actually holds for slightly larger x, as long as e −x 2 /2 ≥ | cos x|, which is true up to x 1 ≈ 0.566π.Actually, the Gaussian form approximates the cosine also beyond this value, since for x 1 ≤ x ≤ π − µ, it holds |f M (x)| = cos M x − e −M x 2 /2 ≤ cos M x ≤ cos M µ, (E3) where we have assumed that M is even, as in the text.Thus the difference decreases exponentially with M up to x = π − µ.
At very small x it is more useful to use the Taylor expansion.There we can notice that at very small |x| the difference vanishes as f M (x) ≈ M 12 x 4 + O(x 6 ), and, in fact, f M (x) ≤ M 12 x 4 .
b 9 a K O M p w B j W 4 A A 9 u o A n 3 0 I I 2 M J j A M 7 z C m 5 M 6 L 8 6 7 8 7 F s L T n F z C n 8 g f P 5 A 3 z M j 5 I = < / l a t e x i t > {t i } < l a t e x i t s h a 1 _ b a s e 6 4 = " b 3 n w I 6 d N o 2 V / P I 5 C a D P Z U Y F 7 0 s I 8 9 I W y a i f n e y S R P N 8 2 / h 8 h + S 7 m B B c O F 8 6 J n m J g A w / s 3 6 5 6 t b

FIG. 2 .
FIG.2.Local density of states as a function of the mean energy for the Hamiltonian(13) for N = 100 spins, g = 1, h = 2, δ = 0.1, and 50 randomly chosen states |k .The circles have been computed with(12) whereas the cross and plus symbols correspond to(27) with r = 0.4 and r = 1, respectively.In all cases x = 3.

FIG. 3 .
FIG.3.Expectation value of the energy for the same model and parameters as in Fig.2.The data show the results of (12) for δ = 0.1 (triangle) and 1 (circles), and (27) for δ = 1 with r = 1 (plus symbols).

FIG. 4 .
FIG.4.Magnetization computed according to(28) for Hamiltonian(13) and N = 20 spins computed exactly (solid lines) and by the Monte Carlo method described in the text with 10 5 samples per point (symbols).The Monte Carlo simulations were done for the full expression (12) with x = 3, for the different values of δ shown in the legend.Upper figure: g = 1, h = 2; Lower figure: g = 2, h = 1.

FIG. 5 .
FIG.5.Magnetization as a function of energy, as in Fig.4but with N = 100.For the figure in the inset, we have subtracted 1/2 + 0.004E/N in order to make the curves more visible.There, the crosses and plus symbols indicate the values obtained by setting a cutoff of 10 −2 in the sampling procedure to resemble an experiment.The upper figure corresponds to g = 1, h = 2, whereas the lower to g = 2, h = 1.

8 FIG. 6 .
FIG.6.Magnetization as a function of the inverse temperature for the Hamiltonian(13) for N = 100 spins, g = 0.3, and h = 0.8 (lower curve) and g = h = 0.4 (upper curve).The solid line represents the exact value M (β), whereas the symbols are obtained with the Monte Carlo method, with 10 5 samples per point and δ = 1 and x = 3.We have discretized the values of the energy from −N to N in intervals of 0.5.The symbols '+' are obtained by setting a cutoff equals to 10 −2 in D δ,k (E).

H 2 FIG. 8 .
FIG.8.Energy density in the non-integrable Ising chain (D1) in the thermodynamic limit.The solid line indicates the energy density of translationally invariant real product states, while the dashed horizontal lines indicate the energy densities of the ground and maximally excited states (estimated numerically with MPS).The colored symbols indicate the states chosen in our numerical simulations.

7 FIG. 9 .
FIG.9.Convergence of Eq. (17b) to the microcanonical values for several local observables, in the non-integrable Ising model (D1).Each column corresponds to a real translationally invariant product state (in order of increasing energy density), determined by θ = π/4 (left), π/3 (center) and π/6 (right).The upper row illustrates the results for the 2-site observable σz ⊗σz, and the lower row for σx, measured, in both cases, in the center of the chain.The dashed lines indicate the microcanonical values in the thermodynamic limit corresponding to the same energy density.Our results show that δ ∝ 1/N 2 converges fast to the microcanonical expectation value, and for δ ∝ 1/N the values are reasonably close.