Classical algorithms for many-body quantum systems at finite energies

We investigate quantum inspired algorithms to compute physical observables of quantum many-body systems at finite energies. They are based on the quantum algorithms proposed in [Lu et al. PRX Quantum 2, 020321 (2021)], which use the quantum simulation of the dynamics of such systems, as well as classical filtering and sampling techniques. Here, we replace the quantum simulation by standard classical methods based on matrix product states and operators. As a result, we can address significantly larger systems than those reachable by exact diagonalization or by other algorithms. We demonstrate the performance with spin chains up to 80 sites.


I. Introduction
Computing the properties of quantum systems in equilibrium is of central interest in many-body physics. For a system at finite temperature, there exists a wide spectrum of techniques that are used in practice. For large systems, where exact solutions are unreachable, they typically approximate the corresponding Gibbs state using variational, sampling, or series expansion methods [1][2][3][4][5][6][7]. For systems at a finite energy, e.g. in the microcanonical ensemble, methods are more scarce [8][9][10].
A possible approach consists in simulating the dynamics and extracting the equilibrium properties from there. For instance, one can use spectral filters in order to retrieve expectation values of an observable O, by averaging them at different times. In this way, one obtains results connected to the diagonal ensemble corresponding to the initial state [11], namelȳ O = n |c n | 2 E n |O|E n (1) where |E n are the energy eigenstates and c n the coefficients of the initial state in that basis. For local Hamiltonians and initial states with finite correlation length (like product states), the values of c n are significant if |E n − E| < δ = O( √ N ), where E is the mean energy of the initial state. Under the eigenstate thermalization hypothesis (ETH) [12][13][14],Ō converges to the equilibrium value in the thermodynamic limit, and it is not necessary to average the results at different time but just to wait for a sufficiently long time. Quantum computers and analog quantum simulators are very well suited for that task, since they can deal with the dynamics of many-body quantum systems in a very natural way [15,16]. Classical methods to simulate the dynamics typically suffer from the linear growth of entanglement [17], which gives an exponential cost with time. Even if the (weak) ETH applies, the thermalization time can be very long [18,19], and this severely restricts the applicability of such classical algorithms.
In this work we propose and analyze a classical algorithm to compute expectation values of the form (1). In particular, for c n that are Gaussian functions of |E n − E| with a variance δ, this can be achieved by simulating the dynamics for a time O(1/δ). This allows us, for instance, to reach δ = O( √ N ) by just using standard time-evolution techniques for tensor networks for very short times, when the entanglement is still very small and thus the techniques work well. One can also reach values of δ = O(1) with modest computational resources. The algorithm is inspired by a quantum algorithm presented in [10] that allows one to compute expectation values of the form (1). This algorithm combines classical sampling (Monte Carlo) techniques with time series and Loschmidt echo-like measurements [20,21] that can be obtained by quantum simulation of the dynamics. Our main modification is to replace the latter by a classical simulation using tensor network states [22][23][24][25][26][27]. This allows us to compute (1) for times O(1/δ) instead of the thermalization time, thus circumventing the problem of entanglement growth. This is done at the expense of having to sample, which just involves the repetition of the whole procedure until convergence. We apply the algorithm to one dimensional systems, sample over a basis of product states and use matrix product states (MPS) and operators (MPO) to simulate the evolution [28]. We illustrate the performance of the method for a non-integrable Ising chain, for which we obtain convergence to the microcanonical values for systems up to 80 sites, far larger than what is possible with exact diagonalization.
Apart from that, in [10] another quantum algorithm was proposed to compute physical observables in a state where an energy filter of width δ is applied. For local Hamiltonians, the computational time also scales as O(1/δ). Here we also analyze a classical algorithm inspired in that method. We notice that this method typically requires much narrower δ (thus longer times) to approach thermodynamic quantities. However, the filtering achieved with a limited evolution time can be optimized if the initial state is chosen with already a reduced energy width. Here we demonstrate this possibility by applying the classical version of the first algorithm on matrix product states found by minimizing the energy variance.
The rest of the paper is organized as follows. In section II we briefly review the concept of energy filters, introduce the filter ensemble and discuss its applications arXiv:2204.09439v1 [quant-ph] 20 Apr 2022 to determine microcanonical and diagonal properties. We also review briefly the quantum algorithms [10] that motivate this work. In section III we discuss the details of a TNS simulation of the quantum algorithms, the different possibilities and associated parameter choices. Section IV presents our numerical results for Ising chains, for each of the algorithms implemented. The paper is closed with the discussion in section V.

II. Filters and quantum algorithms
We start by recalling the definition and properties of the energy filters that are at the basis of the algorithms in [10] and this work.

A. Energy filters
The main tool used by the finite energy algorithms discussed here is a filtering operator that suppresses energy eigenstates outside a target energy interval. In particular, given the Hamiltonian H, we define a Gaussian filter centered in energy E and of width δ as the following op-eratorP Notice that, up to normalization,P δ (E) is a diagonal ensemble in the energy basis. We refer to it as the filter ensemble . ( The corresponding expectation values are precisely of the form (1), with coefficients |c n | 2 distributed according to a Gaussian of width δ. For local Hamiltonians and large systems, product states have that kind of spectral decomposition, with δ ∼ O( √ N ) [29,30] but, since they contain coherences in the energy basis, the corresponding expectation values are very different. An ensemble as (3) could nevertheless be obtained from a product state, but only after evolving and averaging over a long time.
As the width δ is reduced, the filter approaches the microcanonical ensemble. Thus we can make use of the filter to access the microcanonical properties of the quantum system in the following two different ways.
a. Filtering a state. Given a state |ψ , its local density of states (LDOS) is defined as D ψ (E) = ψ|δ(E −Ĥ)|ψ . A broadened version can be computed with the filter as We can also use the filtered state to explore the microcanonical ensemble expectation value O micro (E) of an observableÔ. In generic cases in which ETH is satisfied, and for a value E at which the LDOS does not vanish, will converge to O micro (E) in the limit δ → 0. b. Filtering the whole spectrum. The filter ensemble itself converges to the microcanonical ensemble as the width is reduced. Hence, we can also use it without specifying a state, but directly taking its trace. In particular, the density of states (DOS) of the Hamiltonian H, defined as D(E) = tr δ(E −Ĥ) , can be approximated through the broadening of the δ functions as Moreover, the expectation values in the filter ensemble O δ (E) = tr Ô ρ (E,δ) = tr ÔP δ (E) / tr P δ (E) (7) will converge to the microcanonical values as δ → 0. Whereas in generic cases one can in principle approach the microcanonical expectation values with either (5) or (7), the convergence of the second with δ is much faster, since the filter ensemble is diagonal in the energy basis, while (5) contains contributions from off-diagonal matrix elements [31][32][33], which can converge much slower than the diagonal part [14,34].

Implementing the filter
For the purpose of numerical, but also quantum, simulations of the filter, it is convenient to substitute an approximation for the Gaussian filter. In this paper, following the quantum algorithms in [10], we focus on the cosine filter [35,36] defined aŝ where α is a parameter (with dimensions of energy) that controls the validity of the approximation, and M = α 2 /δ 2 2 with · · · 2 giving the closest smaller even integer. The approximation is valid if Ĥ − E ≤ απ/2, when the spectrum ofĤ lies in one period of the cosine function [37]. (8) can be further approximated by a truncated series of evolution operatorŝ where R = xα/δ , t m = 2m/α, x is a constant that bounds the truncation error in operator norm as With the cosine filter the problem is turned into evolving states or computing the traces of time evolution operators, which leads itself to a natural implementation in quantum simulators.

B. The quantum algorithms
In Lu et al.'s paper [10], two hybrid classical-quantum algorithms corresponding to the two different ways of applying the filter were introduced. We sketch them here for completeness.
The first one computes (5) for a state |ψ that can be easily prepared. Suppose the quantum device can efficiently obtain the following quantities then (5) can be determined by classical postprocessing as without explicitly preparing the filtered state. The required time scale is proved to be a polynomial of system size N , the inverse of the width of filter 1/δ and the inverse of the error, provided the state |ψ can be prepared efficiently, for a value of E in a small interval around the mean energy of |ψ . In the second algorithm (quantum-assisted Monte Carlo), importance sampling is applied to compute (7). Let us rewrite that expression as where { | φ } is an (over-)complete basis, with dµ φ the appropriate measure to ensure the closure relation dµ φ |φ φ| = 1 (a simple choice is for instance the computational basis). D δ,φ is the LDOS defined in (4) and Both D δ,φ (E) and O δ,φ (E) can be obtained by measuring the quantities defined in (11), as long as we can run the first algorithm with the quantum device for the states in the basis { | φ }. Then a Metropolis-Hastings step can be applied classically with regard to the probability distribution D δ,φ (E)/ dµ φ D δ,φ (E), and the value of O δ (E) can be estimated. Because D δ,φ (E) is positive, this method does not encounter a sign problem. Given δ, this second algorithm provides access to observables in the filter ensemble at the cost of simulating time evolutions for times O(1/δ), at the expense of repeating the procedure until the sampling converges. This is especially remarkable because one could obtain a similar result from the time evolution of an initial state with the same distribution of coefficients, but this would require evolving for as long as the thermalization time. As a particular application, if one chooses δ = o( √ N ), the expectation values of intensive quantities will (under the ETH) already be equivalent to those in the Gibbs ensemble at the same mean energy (see III B 1), and therefore this algorithm is an inexpensive way of accessing thermal properties.
If we are interested in microcanonical expectation values, we can use either algorithm, but we need to reduce the width of the filter. Since the trace quantities converge faster in δ to the microcanonical values, a shorter evolution time is required with this second algorithm. In exchange, the procedure needs to be repeated over many states, to perform the classical Monte Carlo sampling.

Extreme values of energy
Both algorithms above rely on the evolution of easily preparable states (a requirement for the single initial state in the first algorithm or the whole sampling basis in the second). A most practical choice is that of product states. The mean energies of such states are contained in an extensive but generally restricted interval within the spectrum [38], such that values close to the edges of the spectrum may be out of reach. As indicated in [10], the accessible range of energies can be extended by considering larger sets of states. In particular, MPS can be used to circumvent the limitation.
For the first algorithm, the initial state |ψ can be found as an MPS with a small bond dimension such that its energy expectation value is close enough to E, as MPS serves as a good representation of the ground state and low-lying excited states. This can be done, for instance, by first finding the MPS with that bond dimension and minimal energy (E min ), and then changing its parameters until the desired energy E > E min is reached. Another possibility is to find the MPS minimizing (H − E) 2 [39,40].
For the second algorithm, a different basis for Monte Carlo sampling can be chosen, where we start from any state |φ 0 whose mean energy is close to E, obtained in the same way, and apply a random Pauli matrix σ x , σ y or σ z on a random site in each proposed move. This strategy gives a complete basis set, as (15) for any state |φ 0 . The change in mean energy is O(1) in each move, and thus this choice of basis ensures enough states for sampling.

III. Classical simulation
The methods that we study in this paper replace the quantum simulation of the dynamics in the quantum algorithms of [10] by classical simulations using tensor networks. The longest evolution time required in (9), which is the deciding factor for the efficiency, is determined by the width of the filter δ as t max = t R ≈ 2x/δ. A quantum simulator should be able to deal efficiently with times t = O(poly(N )), which gives access to δ = Ω(poly(1/N )) [41]. This should be sufficient for both (5) (in the case ETH is satisfied [34]) and (7) to converge to the microcanonical values. On a classical computer, TN techniques provide the possibility to simulate the time evolution of a local Hamiltonian [22,23,28], but the bond dimension required to do so can increase exponentially with time. Thus, starting from a product state, we can simulate times t max ∝ log N with a bond dimension polynomial in system size, which would allow us to efficiently perform classical simulations of the algorithm for δ = Ω(1/ log N ). For the actual implementation of the classical simulation of the dynamics, there exist several options, some of which we discuss in this section.

A. Tensor network implementation
There are multiple different approaches to simulate time evolution with TN techniques (see [28] for a recent review). Some of the most commonly used methods are based on a Suzuki-Trotter approximation of the time evolution operator. One possibility is then to repeatedly apply the approximated short time evolution steps onto a matrix product state (MPS), to obtain a representation of the time-evolved state. Alternatively, the time-evolution operator itself, e −iĤt , can be approximated by a matrix product operator (MPO), constructed also from the iteration of trotterized steps. Therefore, we can use various techniques for the classical simulation of the quantum methods above.
Finding the MPO representation of each term e −iĤtm as an MPO allows us to estimate tr(P δ (E)) as linear combination of the corresponding traces, which for MPOs can be computed very efficiently. This strategy will however fail as we approach the edge of the spectrum for large system sizes while considering small values of δ. The reason is the extremely imbalanced distribution of the DOS D(E), which becomes exponentially small when E is far from the center of the spectrum. For a traceless, local and bounded Hamiltonian, in the thermodynamic limit D(E) converges weakly to a Gaussian distribution with mean energy E = 0 and width proportional to √ N [29,30]: where d is the local Hilbert space dimension and σ 0 is some constant independent of the system size. Thus, for energies ≈ √ N , we expect D(E) to become exponentially smaller than its value at the center. Therefore given all tr e −iĤtm that are reasonably precise, the ratio in (7) could still not be properly achieved all through the whole spectrum: the applicable energy range will be proportional to the square root of the system size, and hence O(1/ √ N ) in energy density, a restriction we also observed in a previous work [9].
Fortunately, this difficulty can be overcome with the importance sampling method described in II B. To begin with, in this method we only need to evaluate the ratio (14) for states for which the probability factor φ|P δ (E) |φ is above some threshold. In particular, when E is away from the center of the spectrum, contributions from the exponentially large maximum of the DOS will be suppressed. Of course, one needs that the chosen basis {|φ } has enough states around the target energy E. Since for our numerical simulations we are free to choose any basis from which we can sample efficiently and whose states can be written as MPS, we can exploit this freedom to try to ensure this condition. A product basis minimizes the cost of the contractions and is often a good choice, since, as mentioned above, it covers an extensive window of the energy spectrum. If this is not the case, we can use any of the methods mentioned at the end of section II to find an MPS with small bond dimension close to the desired energy, and use it to construct a complete basis of the form (15). For the cases we consider in this paper, the computational basis is already an adequate choice, sufficient to produce accurate numerical results over the full spectrum, as we illustrate in the next section.
In the Monte Carlo simulation, the time evolution can be done either at the level of the states, i.e. directly evolving the sampled state as an MPS, or at the level of the operators, i.e. approximating the evolution unitaries for each state as MPOs, storing them in memory, and using them later to do contractions with states randomly sampled from the basis. The second option has the advantage of simulating the dynamics a single time, as the same operators can be reused when doing the sampling over different states. Hence it is faster, but it also consumes much more memory to store the required MPOs [42]. Also, the bond dimension needed to approximate an evolution operator as an MPO is significantly larger than the one used to approximate a time-evolved MPS for the same time. As a concrete example, for system size N = 80, we find bond dimension D MPS = 40 to be enough for evolving MPS, and D MPO ∼ 100 for storing the MPOs, for (δ, α) ∝ (1, √ N ) and 5 × 10 4 samples. The typical time scales taken are 1 week and 1 day, respectively with Intel® Xeon® Gold 6138 processor.

B. Filter parameters and the microcanonical limit
The cosine filter depends on two parameters: the width δ and the period of the filter α. As mentioned, δ determines the maximum time we need to evolve, t max = 2x/δ, while for a fixed δ, α determines the number of terms in the expansion R = xα/δ. Additionally, in the Monte Carlo algorithm we introduce a cutoff parameter and discard states for which the probability is found to be below this threshold.
In this section we discuss the significance of these parameters, as well as which choices ensure approaching the microcanonical limit. Note that the conclusions are valid for both the classical and the quantum version of the algorithms.

Filter width δ
A width δ and mean energy E 0 determine the properties of the filter ensembleP δ (E 0 ). But to estimate the corresponding energy distribution we need to take into account the density of states. Assuming a Gaussian DOS as in (16), the energy distribution of the filter ensemble will also be a Gaussian given by where γ = 1 + δ 2 /N σ 2 0 . We omitted an energyindependent factor in (17). It can be concluded that in the thermodynamic limit, the mean energy and width of the filtered ensemble are given by If we choose δ ∝ √ N , the mean energy of the ensemble is shifted with respect to the parameters of the filter, as explicitly shown in Fig. 4. A filter width that scales as This observation is especially relevant if we are interested in approaching the microcanonical limit: in general, assuming ETH, in the thermodynamic limit a microcanonical energy shell centered at E will yield the thermal values for intensive quantities at energy density E/N if the width ∆ satisfies ∆/N → 0. This condition is already satisfied for the filter ensemble with δ ∝ √ N , which means that the expectation values will converge to the thermal ones, only at shifted energies, according to the previous argument (see Fig. 4). We can similarly estimate the energy distribution of the pure state resulting from the application of the filter onto an individual state |ψ . There is actually a similar argument forP δ (E) |ψ if |ψ is a product state, as such states also have essentially Gaussian LDOS whose widths are proportional to √ N [43]. With a spectral decomposition |ψ = k c k |E k , where |E k are energy eigenstates, the filtered state results as where Γ = 1/ ψ|P δ (E) 2 |ψ is the normalization factor. By choosing the center of the filter at the mean energy of the state E = ψ|Ĥ|ψ =: E ψ , the average energy of the filtered state does not change. Assuming the LDOS of |ψ has a Gaussian form with width σ ψ √ N , where σ ψ is independent of system size, the energy variance of |ψ can be estimated through substituting the sum over eigenstates by an integral over energy values with Gaussian weights. We obtain Again, we may want to consider how this affects approaching the microcanonical limit as the width of the filter is decreased. A major difference in this respect between the filter ensemble and the filtered state is that the second contains coherent contributions from different energy eigenstates. Thus (5) includes contributions from off-diagonal matrix elements in the energy basis, which only become negligible when the width of the energy distribution decreases sufficiently fast with N . More concretely, from canonical typicality arguments we can expect that, for non-integrable systems, the expectation value of a local observable converges to the thermal value when the energy deviation of the state decreases as a polynomial of 1/N [34]. In [35] we observed a trend to convergence already with a slower decrease ∼ 1/ log(N ). For these scalings of the filter width, according to (20), the width of the filtered state will scale in the same way.

Period of cosine filter α
Different to the Gaussian one, the cosine filter (8) is periodic, but it remains a good approximation of the former when the argument is bounded within one period. More concretely, operatorsP δ (E) andF δ,α (E) are close to each other when Ĥ − E ≤ απ/2. At the same time, because the number of terms that need to be evaluated in the sum (9) is proportional to α, it is convenient to choose the smallest possible value that ensures the previous property. For a local Hamiltonian, a value α ∝ N is enough for the condition to hold for all values of E within the energy spectrum. If the operator acts only on a limited energy window, a smaller value of α can be chosen, as long as all relevant states are almost supported in [E − απ/2, E + απ/2]. This can be used, for instance, when the filter acts on a product state, whose energy distribution is approximately Gaussian, with support on an energy interval ∝ √ N . If additionally, the filter is centered near the mean energy of the state and δ = o( √ N ), it is enough to choose α ∝ √ N [10]. In practice, we find α = 3 max(σ φ √ N , δ) to work well for all system sizes.

Monte Carlo cutoff threshold
For the discussions in III B 2, it should be ensured that the samples in Monte Carlo simulations not stepping into other energy periods of the cosine filter when α ∝ √ N . In other words, the weights of the states whose mean energy are close the edges of [E − απ/2, E + απ/2] should be small enough. A cutoff threshold can be applied to the weights of samples to improve numerical stability in Monte Carlo simulations. To be more concrete, a proposed state |φ will be directly discarded (its probability assimilated to 0) if D δ,φ < . Besides making the numerics more stable, the presence of the cutoff prevents the peak of the DOS from shifting the ensemble when we target energies near the edges, because it restricts the visited energy range in more general cases, as we show next.
If we consider a product state basis, an individual state |φ will have mean energy E φ and width σ φ √ N . Again, assuming a Gaussian distribution, we can estimate its weight in the sum as When δ = Ω(

Choosing the parameters for microcanonical values
According to the discussions above, we can summarize in Table I some possible choices of parameters for the various algorithms, such that we obtain convergence to thermal values in the thermodynamic limit. The table shows the scaling with system size of α and δ, as well as the resulting cost (in terms of maximum evolution time and number of evolutions to run) and the energy range where the methods are applicable. The fastest method [(I) in the table] is directly computing (7) by taking the traces of the evolution operators approximated as MPOs, but, as discussed in section III A, it is only applicable in an energy interval of width proportional to √ N around the center of the spectrum. With Monte Carlo sampling (II), in contrast, it is possible to reach the whole spectrum. The filter width required to obtain convergence to thermal values in the thermodynamic limit should scale at most as O( √ N ). A larger width corresponds to a shorter evolution time t max , and hence a smaller bond dimension required for the MPS or MPO, but it also shows slower convergence to the thermal values as the system size is increased. We thus show two possible choices (IIa) and (IIb), both of which we explore numerically, using different approaches for time evolution, in IV A. In (IIa), a cutoff threshold is applied in the Monte Carlo simulations to avoid the energy shift due to DOS when δ ∝ √ N . Finally, when filtering a state (III), the time t max required to approach microcanonical values is polynomial in the system size, which means that with TN techniques we will be able to extract microcanonical values with this method only for small system sizes. For this last method, however, the achievable width can be optimized by applying the algorithm on a state with reduced energy width. We present a way to implement this improvement by using MPS obtained after a variational minimization of the variance.

IV. Results
To demonstrate and benchmark the various methods described above, we apply them to a quantum Ising chain with open boundary conditions, The model is integrable if either g = 0 or h = 0.
Here we choose a particular set of parameters (J, g, h) = (1, −1.05, 0.5) far from integrability [44]. In the thermodynamic limit, the corresponding energy density lies in the interval E/N ∈ [−1.33, 1.72]. For the observable, we focus on the average magnetization State filtering poly(1/N ) Ω( √ N ) poly(N ) poly(N ) TABLE I. Possible choices of the filter parameters (δ, α) that ensure approaching the microcanonical values in the thermodynamic limit for the various methods, and corresponding maximum evolution time and number of steps.

A. Filter ensemble
We start by illustrating the performance of the Monte Carlo algorithm to estimate expectation values in the filter ensemble (3) at all values of energy. Since the largest time we need to simulate is t max ∝ 1/δ, the classical simulation can efficiently treat widths δ = O(1/ log N ) and larger.
For the numerical benchmarking, we choose widths δ ∝ √ N and δ = const, which, according to the discussions above, are enough to approach the microcanonical values in the thermodynamical limit [see (II) in table I].
In the non-integrable model we consider, the values are thus expected to converge to the thermal ones. Thus, we can compare the results of the algorithm with the exact values in thermal equilibrium at the corresponding energies, which we can compute independently using standard TN techniques [6].
The calculations can be done using different options for the TN evolution (section III A). As long as the results are converged in bond dimension, both approximating the evolution operators as MPO or the individual evolved states as MPS are valid strategies, and we show results obtained with both of them. Figure 2 demonstrates the success of the method to find expectation values in the filter ensemble, for system sizes up to N = 80. In particular, for this plot, we chose to simulate and store the MPOs for all evolution operators before realizing the sampling over the computational basis. For the filter, we used filter parameters (δ, α) ∝ ( √ N , N ), which, as argued in section III B, in a generic case is enough for the observable to converge, in the thermodynamic limit, to the thermal expectation value if introducing a cutoff threshold (section III B 3). We choose to sample over the computational basis and the cutoff threshold = 10 −4 D δ,φ0 (E), where |φ 0 is the initial state in the Monte Carlo simulation, obtained by minimizing φ|(Ĥ − E) 2 |φ for |φ in the basis set.

MPO version of Monte Carlo simulation
As shown in the upper panel of fig. 2, the results clearly converge to the thermal value as the system size N is increased or δ is reduced. The inset shows explicitly this convergence for energy density E/N = 1.44, relatively close to the edge of the spectrum. Errors have two main sources, which are shown in this plot: the statistical error   from the Monte Carlo sampling (error bars) and the truncation error from the finite bond dimension of the MPO (shown as shadowed region). Only for the largest system size N = 80 and smallest width δ = 0.5 √ N we observe a small discrepancy, but compatible with our estimated errors from both sources.
The lower panels of fig. 2 show explicitly the convergence of the Monte Carlo sampling at the same energy density E/N = 1.44, for various system sizes and bond dimensions, and for two different values of the width. In all these cases we observe that, after 50000 steps, the results are practically converged, even though some fluctuations can be appreciated.

MPS version of Monte Carlo simulation
To illustrate the performance of the algorithm when individual states, rather than operators, are evolved, we choose a narrower filter width, which should lead to values closer to the thermal ones, while still being reachable by classical simulations. Notice, however, that similar values could have been obtained with the MPO option, at a different cost in memory and time.
In particular, considering a constant value δ = O(1), independent of system size, requires evolution until a constant time, which (for not too small values of δ) can be efficiently simulated using TN. Thus we choose parameters (δ, α) ∝ (1, √ N ) [(IIb) in Table I], and explore two values of the energy density, one near the center of the spectrum (E/N = 0.72) and one close to the edge (E/N = 1.44), and both within the reach of product states from the computational basis, which we take again as our sampling basis.
Results are shown in Fig. 3. In the upper panels, we plot the difference between (m z ) δ and the corresponding thermal value as a function of the system size, for two values of the bond dimension, with error bars indicating the statistical error. We find that, within error bars, the distance to the thermal value decreases as the system size grows, with a relative difference smaller than 0.5% for N = 80 and E/N = 1.44.
In the lower panels of Fig. 3, we again show explicitly the convergence of the Monte Carlo sampling. Note that here the same seed for randomization was used for different bond dimensions, so the fact that the solid and dashed lines (representing D = 40 and D = 60 respectively) are on top of each other indicates the convergence with regard to bond dimension already at D = 40.
As illustrated above, the combination of short-time dynamics simulation and sampling provides a powerful method to compute the expectation values in the filter ensemble, as long as the basis contains vectors with substantial weight in the energy region of interest. For the computational basis that we have used in the examples, mean energies lie in the interval E/N ∈ [−1, 1.5]. According to (21), the weights of basis states D δ,φ will decay exponentially with the system size for a fixed energy density E/N outside this interval. Fig. 2 shows that, for system size N = 80, the Monte Carlo sampling with computational basis remains valid at E/N = 1.68 (rightmost point) and E/N = −1.2 (leftmost point), much closer to the edges of the spectrum, so that we do not need to resort to the Pauli basis mentioned in II B. Using this 20 40  basis may however become necessary as we keep increasing the system size, or if we consider other models or higher dimensions. %it to finally fail as keeping increasing the system size, or when we deal with other models with larger gap between product states and the ground (or maximally excited) state. In these cases the Pauli basis mentioned in II B will be applicable.

Exploring the center of the spectrum without sampling
The fastest alternative to evaluate (7) with TN simulations is to directly evaluate numerator and denominator from traces of the evolution MPOs, without the sampling iteration [(I) in Table I]. As discussed above, this is only feasible in the central region of the spectrum, over a width ∝ √ N , before the density of states becomes exponentially small. If the filter width is not much smaller than this scale, the mean energy of the filter ensemble will be effectively shifted towards the maximum of the DOS, as explicitly computed in sec. III B 1. Figure 4 illustrates the behavior of this alternative for (δ, α) ∝ ( √ N , N ). The upper plot shows the results obtained for the magnetization m z for various system sizes and filter widths, as a function of the energy density corresponding to the center of the filter. Because of the shift discussed above, the results do not converge to the thermal ones (indicated by the solid line) at that energy,  but at a shifted value according to (18) (indicated by a dashed line for each δ).
We observe that, while in the central part of the spectrum better convergence is observed as δ decreases or N increases, near the edge of the spectrum the method fails to give the correct microcanonical values, especially for larger systems and smaller widths, as then it becomes sensitive to the exponentially smaller density of states. This is also visible in the lower panel of the figure, where we plot the difference between the computed values O δ (E) ad the thermal ones at the correspond-ing shifted energies. According to (17), the denominator of (7), tr P δ (E) should scale as exp(−E 2 /2γN σ 2 0 ), indicating the range of energy densities for which the method actually converges shrinks as 1/ √ N .

B. Filtered pure state
As discussed in section III B, we can apply the filter on a state to decrease its energy variance and, in the generic case, obtain convergence to the microcanonical properties. Also in this case, the largest time that needs to be simulated is t max = 2x/δ, which can be done efficiently by TN simulations if the width is at least O(1/ log N ). But, in contrast to the calculations for the filter ensemble, diagonal in the energy basis, a much smaller δ is required in this case to approach microcanonical values in the thermodynamic limit. More concretely, δ = O(poly(1/N )) should be enough, but this requires t max = Ω(poly(N )), and thus a bond dimension that increases exponentially with N . Additionally, when filtering a product state, the total number of terms to probe will be (2R + 1) 2 , where R = xα/δ, which grows at least as N 3 .
One way to mitigate the second problem is to apply the filter on a state with already reduced energy variance. This allows us to choose a smaller period α and correspondingly keep a more moderate value of R and to test the strategy for moderate sizes. We implement this strategy using as initial states MPS with a given bond dimension, found by variationally minimizing Ĥ − E 2 at the value of E we are interested in. To test this strategy, we have targeted a value E/N = 0.72 for system sizes 20 ≤ N ≤ 80, and obtained initial MPS with reduced widths from the minimization of (H − E) 2 with bond dimensions D 0 ∈ {1, 2, 5, 10}. For each one of this states, we compute the width σ D , and then apply a filter with parameters δ = σ D /2 √ N and α = 3σ D . We show the results in Fig. 5. To analyze the convergence towards the thermal value as the width decreases, we plot the relative error of the magnetization with respect to the thermal one as a function of 1/δ (left panel) and 1/(N 2 δ) (right panel) for each system size. The first case shows no clear scaling laws, which indicates that δ = O(1) is not enough for convergence in the thermodynamic limit. The right panel, instead, exhibits a trend to convergence for δ = O(1/N 2 ), even though numerically it becomes difficult to reach the lower right corner for large systems.

V. Summary and Discussion
We have presented a quantum-inspired classical method, based on a TNS simulation of the quantumassisted Monte Carlo algorithm proposed in [10], that allows us to compute microcanonical and diagonal values for quantum many-body systems. Our method estimates broadened spectral functions, which takes the form of the trace of an energy filter operator, or a product of the latter with an observable, via sampling over time-evolved product states. Because the longest required time is proportional to the inverse filter width, that is very short and easy to simulate with tensor networks, it allows us to find expectation values in a diagonal ensemble which would only be reached after a much longer time evolution from an initial product state. While filter widths O( √ N ) are enough to find the diagonal ensemble values of generic product states, we can also reach energy filters of constant widths, as these only require O(1) evolution times. These scalings are enough to obtain convergence to thermal equilibrium in the thermodynamic limit, in the generic case.
We have benchmarked the algorithm on the nonintegrable Ising chain, for sizes up to N = 80 sites (far beyond the reach of exact diagonalization), and we have checked different choices of the parameters that affect the efficiency and applicability of the method. In particular, we explicitly show diagonal expectation values for Gaussian ensembles of width O( √ N ), obtained with low computational cost over the whole range of energies, and observe their convergence towards the thermal equilibrium values. Reducing further the width, we obtain microcanonical expectation values with high precision.
We can also classically simulate the provably efficient quantum algorithm in [10], in which the filter is applied on a fixed initial state. If we want to use this method to explore the microcanonical properties, nevertheless, the filter width needs to decrease with the system size in order to guarantee convergence in the thermodynamic limit, which results in increasing times and an exponentially growing bond dimension. We have however optimized the procedure by choosing as initial states MPS with minimal variance. In this way, we can run the algorithm and observe convergence for reasonably sized systems.
The results shown in this paper demonstrate the potential of the algorithm. Accessing the microcanonical values would be helpful to investigate all sorts of out-ofequilibrium quantum many-body behavior, for instance many-body localization [45] and quantum scars [46], for large systems. Although we have only implemented it for a translationally invariant spin chain with short-range interactions, the method can be easily extended to systems with disorder, or long-range interactions, and also to bosonic or fermionic systems. We have recently learned that A. Schuckert et al. are using a related method to study long-range interacting models [47].
In principle, the same method can be also extended to higher dimensional systems, using PEPS (projected entangled pairs states) [48], which is a promising possibility, given the absence of numerical methods for extracting microcanonical expectation values in that case, beyond exact diagonalization. The numerical challenge is then higher, due to the larger computational cost of the corresponding algorithms, and determining the preferable implementation option should be analyzed. Also further extensions are possible that consider other filter functions or combinations of filters.