# Accuracy of the typicality approach using Chebyshev polynomials

Henrik Schlüter , Florian Gayk , Heinz-Jürgen Schmidt , Andreas Honecker and Jürgen Schnack

## Abstract

Trace estimators allow us to approximate thermodynamic equilibrium observables with astonishing accuracy. A prominent representative is the finite-temperature Lanczos method (FTLM) which relies on a Krylov space expansion of the exponential describing the Boltzmann weights. Here we report investigations of an alternative approach which employs Chebyshev polynomials. This method turns out to be also very accurate in general, but shows systematic inaccuracies at low temperatures that can be traced back to an improper behavior of the approximated density of states with and without smoothing kernel. Applications to archetypical quantum spin systems are discussed as examples.

## 1 Introduction

The (numerically) exact evaluation of thermodynamic quantum equilibrium observables is restricted to small systems due to the exponential growth of the Hilbert space for systems with finite-size single-site Hilbert spaces such as Heisenberg or Hubbard models. For quantum systems with unrestricted single-site spaces, the situation is even more severe. Only very few analytically solvable systems are known which creates a massive need for numerical (approximation) schemes. One rather successful means to approximate thermodynamic quantities rests on trace estimators which approximate a trace by an expectation value with respect to a random vector [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12]. These schemes, sometimes also called typicality or (microcanonical) thermal pure quantum states [13], [14], [15], [16], have been used very successfully, in particular, in the field of correlated electron systems, see e.g., [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33], [34], [35], [36], [37] but also in quantum chemistry [38, 39].

A prominent formulation of the method is the finite temperature Lanczos method (FTLM) [4, 21, 40], [41], [42] which employs a Krylov space expansion for exp { β H } (operators are marked by a tilde). It turns out that FTLM produces very accurate approximations when estimates are averaged over random vectors (order of ∼100, fewer for larger spaces); compare [10, 16, 31, 36, 43], [44], [45].

Despite this success, the authors of [8] suggest that an alternative approximation using an expansion of the density of states in terms of Chebyshev polynomials should be more accurate [8]. The major argument is that this expansion does not suffer from the loss of orthogonality during recursive state generation used in Krylov space methods. This property is certainly responsible for the high accuracy obtained in numerical unitary time evolution using a Chebyshev expansion, see e.g., [8, 46], [47], [48], [49], [50].

In the present paper, we therefore study several Heisenberg quantum spin systems and derive numerical as well as formal conclusions about the accuracy of the method. We can summarize that the approach via Chebyshev polynomials is indeed accurate, but not more accurate than FTLM [44]. On the contrary, under certain circumstances the employed kernel which smoothens (unphysical) oscillations of the approximated density of states introduces systematic inaccuracies. The same holds for the mapping of the energy spectrum onto the interval [−1 + ɛ/2, 1 − ɛ/2] to comply with the domain of definition of the polynomials.

The paper is organized as follows. In Section 2, we recapitulate the Chebyshev method. In Section 3, we present our numerical examples. The article closes with a discussion in Section 4.

## 2 Method

In this section, we briefly introduce the Chebyshev method and its parameters to be able to discuss the method’s accuracy. For a more detailed description of the algorithm we recommend [8].

In a quantum mechanical system with a discrete energy spectrum, the microcanonical density of states is defined as

(1) ρ ( E ) n δ ( E E n ) .

The canonical partition function Z(β) is determined by the integral over the density of states weighted with the Boltzmann factor:

(2) Z ( β ) = n e β E n = ρ ( E ) e β E d E

with β = 1 k B T . Correspondingly, the heat capacity is evaluated as

(3) C ( β ) k B = β 2 1 Z ( β ) ρ ( E ) e β E E 2 d E 1 Z ( β ) ρ ( E ) e β E E d E 2 .

For the susceptibility, we use the S z -symmetry of Heisenberg systems and decompose the density into contributions from all orthogonal subspaces with total magnetic quantum number M, i.e.,

(4) χ ( β ) ( g μ B ) 2 = β 1 Z ( β ) M M 2 ρ ( E , M ) e β E d E 1 Z ( β ) M M ρ ( E , M ) e β E d E 2 .

We further calculate only contributions for M ≥ 0, since the respective contributions for negative M are degenerate and can be added accordingly.

The idea of the Chebyshev algorithm is to expand the microcanonical density of states ρ(E) in terms of Chebyshev polynomials and then approximate the integral (2) by Gauss–Chebyshev integration. We would like to state already at this stage that some accuracy problems shown later in this article arise if the approximated density of states does not behave like a proper density, e.g., if it becomes negative.

Since the Chebyshev polynomials are restricted to the interval [−1, 1], a variable transformation of the Chebyshev polynomials to arbitrary intervals must be introduced as in [51, Section 1.3.2]. The transformation of the Hamiltonian results in

(5) H 1 m ( H c 1 )

with

(6) m = E max E min 2 ε , c = E max + E min 2

where, as suggested in [8], a parameter ɛ is introduced to prevent truncation of the approximated delta peaks corresponding to the extremal eigenvalues. The original energy interval is thus scaled to the interval [ 1 + ε 2 , 1 ε 2 ] .

The corresponding scaled density of states ρ ̄ ( x ) is then expanded in terms of Chebyshev polynomials C n (x):

(7) ρ ̄ ( x ) 1 π 1 x 2 μ 0 + 2 n = 1 N  deg μ n C n ( x ) .

It can be shown that the coefficients of the expansion are given by the traces

(8) μ n = Tr C n ( H ) .

These traces are approximated using the typicality approach, i.e.,

(9) μ n dim (H) R r = 1 R r | C n ( H ) | r r | r = Θ n ( R ) ,

with

(10) | r = ν r ν | ν

being a random vector with Gaussian distributed components r ν with respect to a chosen orthonormal basis {|ν⟩}. The relative error of an estimate Θ n (R) is proportional to 1 / R d i m ( H ) as shown in e.g. [8], where R is the number of random vectors and dim ( H ) the dimension of the Hilbert space. At this point it should be noted that these traces can also be evaluated to numerical accuracy using a complete basis. This possibility will be used to distinguish statistical and systematic deviations later on.

Due to the finite order of the expansion, so-called Gibbs’ oscillations can occur which causes the approximated density of the states to have negative values. If one wants to obtain a “physical” representation of the density, i.e., without negative values, one can modify the coefficients μ n by a kernel [8]. The kernel fixes this problem at the cost of introducing a systematic error which vanishes for N deg → ∞. In this paper, we restrict our discussion to the use of the Jackson kernel whose coefficients read

(11) g n = 1 N deg + 1 ( N deg n + 1 ) cos π n N deg + 1 + sin π n N deg + 1 cot π N deg + 1 N deg + 1 .

In figure captions or legends, we will write g n = JK when the kernel is applied, otherwise g n = 1.

For an arbitrary function f(x), the Gauss–Chebyshev integration gives rise to the approximation

(12) 1 1 f ( x ) 1 x 2 d x π N ̃ k = 1 N ̃ f ( x k ) , see  [ 52 ] ,

(13) x k = cos π ( k 1 2 ) N ̃ , k = 1 , , N ̃ .

This approximation is an exact identity if f(x) is a polynomial of order 2 N ̃ 1 or smaller [51]. In the case at hand, f(x) has to be chosen as

(14) f ( x ) g 0 μ 0 + 2 n = 1 N  deg g n μ n C n ( x ) e β ( m x + c ) ,

and is thus no polynomial of order smaller than 2 N ̃ 1 . However, the approximation through Gauss-Chebyshev integration is still a good choice for numerical purposes as it can be computed through a discrete cosine-transform (type III). Deploying the f(x) to Eq. (12) gives

(15) Z = ρ ( E ) e β E d E

(16) = 1 1 ρ ̄ ( x ) e β ( m x + c ) d x

(17) 1 N ̃ k = 1 N ̃ γ k e β ( m x k + c ) ,

where the values of the weights γ k read

(18) γ k π 1 x k 2 ρ ̄ M ( x k )

(19) = g 0 μ 0 + 2 n = 1 N deg g n μ n cos n π ( k 1 2 ) N ̃ .

If one chooses N ̃ N deg , the sum can be complemented to an upper limit of N ̃ with additional g n = 0 terms. The γ k can then be computed through a discrete cosine-transform (type III) of the coefficients g n μ n which allows a faster computation of the sum. The time needed scales with N ̃ ln N ̃ instead of N ̃ N deg [8].

If there are known symmetries, the scheme can be performed for each orthogonal subspace H Γ separately. The approximated partition function can then be written as

(20) Z ( β ) = 1 N ̃ Γ k = 1 N ̃ γ k Γ e β ( m Γ x k + c Γ ) .

All of the following systems possess S z symmetry which implies orthonormal subspaces corresponding to the quantum number of the total magnetization Γ = M. In our numerical examples, subspaces with dimension D < 15,000 are fully diagonalized for larger systems. For smaller systems, i.e., those containing only subspaces of dimension D < 15,000, only subspaces with dimension D < 1000, are fully diagonalized. In a real application, one would of course diagonalize all subspaces numerically exactly where this is possible.

To summarize, the Chebyshev method depends on several parameters that can have an effect on the accuracy of the results. These are the order of the expansion N deg, the scaling parameter ɛ, the number of random vectors R, the use of a kernel g n , and the number of supporting points N ̃ .

## 3 Numerical results

The Chebyshev algorithm uses random vectors for trace estimation, compare Eq. (9), therefore the results are expected to exhibit a statistical distribution. To assess this statistical behavior, we perform two kinds of studies: (A) We investigate a thermodynamic observable as a function of the number R of random vectors used for the trace estimator (9), and (B) we study the variance among P realizations per fixed parameter set (for some of the cases presented below).

By considering each realization O i (β) as a random measurement, a mean O ̄ ( β ) and a variance δO(β)2 can be defined:

(21) O ̄ ( β ) = 1 P i = 0 P 1 O i ( β ) ,

(22) δ O ( β ) 2 = O 2 ̄ ( β ) O ̄ ( β ) 2 .

If additionally an exact result O E is known, the systematic deviation

(23) Δ O ( β ) = | O ̄ ( β ) O  E ( β ) |

can be defined.

As specific quantum spin systems we investigate three archetypical systems that show fundamentally different behavior at low temperatures, namely a spin ladder that is gapped in the thermodynamic limit [53], a spin chain that is gapless in the thermodynamic limit [54], and a sawtooth chain in the vicinity of a quantum-critical point [55], [56], [57]; compare Figure 1.

Figure 1:

Systems investigated in Section 3.1, from top to bottom: Ladder, chain, sawtooth chain. Periodic boundary conditions will be applied.

In this subsection, the accuracy of the Chebyshev algorithm is investigated using a Heisenberg ladder for various numbers of spins N with spin quantum number s = 1/2 and periodic boundary conditions. The Hamiltonian reads

(24) H = J 1 i s i , 1 s i , 2 + J 2 i , j s i , j s i + 1 , j + g μ B B i , j s i , j z .

where the first subscript i ∈ {1, …, N/2} of the spin operators denotes the rung and the second subscript j ∈ {1, 2} denotes the leg of the spin. Thus, the exchange interaction J 1 connects the nearest neighbor spins on rungs, and J 2 does the same on legs. Both are chosen to be antiferromagnetic, J 1 = J 2 = 1.

Table 1 shows the standard configuration of parameters used in the following. The order of expansion N deg and the number of random states R are chosen for low computation times and sufficiently accurate results. Their influences on the accuracy are discussed in Section 3.1.2a and Section 3.1.1. As the parameter ɛ is introduced as a corrective variable, its influence will be shown separately in Section 3.1.2b and is omitted for the time being. The same argument holds for the kernel g n , shown in Section 3.1.2d. To make use of the discrete cosine-transform (type III), the number of points of integration N ̃ has to be greater than or equal to N deg. The equality is chosen as a starting point.

Table 1:

Standard configuration of parameters used in the following calculations.

R N deg ɛ N ̃ g n
200 100 0 100 1

#### 3.1.1 Statistical deviations

Since the approximation of the traces by using random vectors is the cause of the statistical variations in the result, the number of random vectors R is varied to investigate the latter. The values given in Table 1 are used as a standard configuration. In Figures 2 and 3, the heat capacity and susceptibility for various values of R are plotted next to the result determined by exact diagonalization. In addition, a result determined by the Chebyshev algorithm is shown where the traces μ n are computed numerically exactly instead of approximating them using random vectors.

Figure 2:

The heat capacity of the Heisenberg ladder with N = 16 spins s = 1/2 at zero field in linear and logarithmic plots computed using the Chebyshev algorithm in the standard parameter configuration, see Table 1, for various values of R (colored curves) and with exact diagonalization (ED). The purple curve shows a result of the Chebyshev algorithm where the calculation of the traces μ n is done using a complete basis.

Figure 3:

The differential magnetic susceptibility of the Heisenberg ladder with N = 16 spins s = 1/2 at zero field in linear and logarithmic plots computed using the Chebyshev algorithm in the standard parameter configuration, see Table 1, for various values of R (colored curves) and with exact diagonalization (ED). The purple curve shows a result of the Chebyshev algorithm where the calculation of the traces μ n is done using a complete basis.

One can see that both the heat capacity (Figure 2) as well as the differential susceptibility (Figure 3) match the respective curve derived from exact diagonalization very well for T/|J 1| > 10−2. However, for R = 10 a noticeable deviation in the maximum of the main peak of both observables can be seen. Additionally, all curves of the heat capacity show a “ghost dip” at low temperatures for T/|J i | ≈ 10−3 − 10−2. Nevertheless, for most purposes the achieved accuracy for the standard parameter configuration and R > 100 is more than sufficient at higher temperatures.

It is noticeable that in the region of the “ghost dip”, all approximate curves deviate from the exact solution independently of R, suggesting a small statistical but significant systematic error. This is confirmed by the fact that the curve determined with numerically exact traces shows this deviation as well.

#### 3.1.2 Systematic deviations

Next, we discuss how tuning the parameters N deg , ε , N ̃ , and g n affects systematic deviations. This is mostly done by observing the behavior of the “ghost dip” of the heat capacity under variation of each parameter.

1. The order of expansion can be increased to push the systematic deviations, i.e., the “ghost dip”, to lower temperatures. This is demonstrated in Figure 4. Computation time increases linearly with N deg.

2. The scaling parameter ɛ seems to have a negative rather than the intended positive effect on the heat capacity as shown in Figure 5. The best result is achieved for ɛ = 0. Conversely, when considering the results for the density of states of the largest subspace with M = 0 (see Figure 6), one can see that the parameter ɛ has the intended effect to prevent the peaks of the lowest and highest eigenvalues to be cut off. There are also situations where a ɛ ≈ 10−6 seemingly decreases the depth of the dip compared to ɛ = 0. Anyhow, such small values of ɛ are not sufficient to prevent the cut-off of the density of states, and the improvement was not significant when compared to statistical deviations.

3. The number of supporting points N ̃ is investigated in Figure 7. It is difficult to give universal recommendations regarding this parameter. Our experience is that it should be chosen equal to the order of expansion N deg. Other systems could show very different results, but it is important to note that a good choice of N ̃ scales with N deg.

4. The smoothing with the Jackson kernel causes the unphysical ghost dip to become a ghost peak, see Figure 8, but it does not seem to improve the result. It could even produce a negative effect as a dip (negative C) is more easily identified as an error than a peak. A good approach could be to always compare both the smoothened and the native result. This can be done without great computational effort.

Figure 4:

The heat capacity C/k B of the Heisenberg ladder with N = 16 spins s = 1/2 computed using the Chebyshev method in the standard parameter configuration, see Table 1, for various N deg (colored curves) compared to exact diagonalization (ED).

Figure 5:

The heat capacity C/k B of the Heisenberg ladder with N = 16 spins s = 1/2 computed using the Chebyshev algorithm in the standard parameter configuration, see Table 1, for various values of ɛ (colored curves) compared to exact diagonalization (ED).

Figure 6:

Detail of the scaled density of states ρ M (x) of the Heisenberg ladder with N = 16 spins s = 1/2 on the subspace with M = 0 computed using the Chebyshev algorithm in the standard parameter configuration, see Table 1, and the scaled ground state eigenvalue obtained by the Lanczos algorithm for various values of ɛ.

Figure 7:

The heat capacity C/k B of the Heisenberg ladder with N = 16 spins s = 1/2 computed using the Chebyshev method in the standard parameter configuration, see Table 1, for various N ̃ (colored) and with exact diagonalization (ED).

Figure 8:

The heat capacity C/k B of the Heisenberg ladder with N = 16 spins s = 1/2 computed using the Chebyshev algorithm (colored) in the standard parameter configuration, see Table 1, for various N deg with and without Jackson kernel compared to exact diagonalization (ED).

For larger systems a kernel can have an even stronger negative effect as is demonstrated in Figure 9 for N = 24 and s = 1/2. The Jackson kernel is significantly setting back the convergence of the expansion. For an order of N deg = 100, which produces very accurate approximations without kernel, the application of the kernel renders the result to become unusable, compare top of Figure 9. One needs to expand the polynomial to an order of N deg = 500 to counteract the inaccuracy introduced by the kernel, but even then the result with kernel is still not significantly better than the result without the kernel. The result without the kernel seems already sufficiently accurate for k B T/|J 1| < 10−2 and N deg = 100.

Figure 9:

The heat capacity C/k B of the Heisenberg ladder with N = 24 spins s = 1/2 computed using the Chebyshev method (CM) and exact diagonalization (ED). Displayed are the results with and without Jackson kernel (JK) for N deg = 100 and N deg = 500.

### 3.2 Heisenberg ring

Since the Heisenberg ladder is a gapped system, i.e., a spin system with a non-zero excitation energy between the ground state and the first excited state in the thermodynamic limit, we would also like to investigate a system that is gapless in the thermodynamic limit. The behavior of thermodynamic functions of the system at low temperatures highly depends on this excitation energy. One could argue that the deviation shown in the previous section are due to this dependency. Hence, in this section we will discuss an antiferromagnetic Heisenberg ring with s = 1/2 for which the excitation energy vanishes in the thermodynamic limit.

In Figure 10, the results for the Heisenberg ring with N = 24 spins are displayed. One can see that the deviations here are very similar to the ones for the Heisenberg ladder with the same system size and choice of parameters, compare Figure 9, even though here, they occur at slightly lower temperatures.

Figure 10:

The heat capacity C/k B of the Heisenberg ring with N = 24 spins s = 1/2 computed using the Chebyshev method (CM) and exact diagonalization (ED). Displayed are the results with and without Jackson kernel (JK) for N deg = 100 and N deg = 500.

To further investigate the influence of the gap’s size on the deviations in the results, the curves for the heat capacity per site for various numbers of spins are displayed in Figure 11. While the differences of the results between k B T/|J| = 2 × 10−2 and 2 × 10−1 are mostly due to finite-size effects, so even the exact curves would deviate from each other, the differences between k B T/|J| = 10−3 and 10−2 on the other hand are due to the method’s inability to reproduce the low temperature behavior of the Heisenberg ring for various system sizes. However, there is no definite trend of better results for larger systems recognizable, see the lower graph in Figure 11. Thus, the inaccuracies of the results do not directly depend on the gap size.

Figure 11:

The heat capacity per site C/(Nk B) of the Heisenberg ring with s = 1/2 computed using the Chebyshev method (CM) for N = 16, 20, 24, 28 and exact diagonalization (ED) for N = 24.

The antiferromagnetic Heisenberg ring with N = 10, s = 5/2 and the nearest neighbor interaction is an interesting example as well, as this system is realized as a magnetic molecule (abbreviated Fe10) called the “ferric wheel” [44] that can be accurately described by this model. In Figure 12, P = 100 estimates with the Chebyshev method using R = 1 random vectors and their mean are compared to an estimate using R = 100 random vectors. They are displayed with and without kernel.

Figure 12:

The heat capacity C/k B of the Heisenberg ring with N = 10 spins s = 5/2 computed using the Chebyshev method (CM) and exact diagonalization (ED). Displayed are the results of P = 100 realizations with R = 1 (light blue curves), their average (dark solid curve), the exact result (red solid curve), and one realization (P = 1) with R = 100 (dotted curve); N deg = 100.

One can see that without kernel the estimates with R = 1 are broadly scattered at low temperatures while their mean and the estimate with R = 100 random vectors are almost perfectly aligned with the exact result. When employing the kernel the estimates with R = 1 are distributed less broadly but their mean and the estimate with R = 100 deviate strongly from the result of the exact diagonalization. From the experience collected before, we assume that these deviations can be resolved by using higher orders of expansion N deg, see Figure 13 with a linear temperature axis. Note that even for N deg = 200 the result without kernel outperforms the one with kernel. Furthermore, the result with N deg = 100 without kernel is more accurate than the result with N deg = 200 and kernel.

Figure 13:

The heat capacity C/k B of the Heisenberg ring with N = 10 spins s = 5/2 computed using the Chebyshev method (CM) and exact diagonalization (ED). Displayed are the results with and without Jackson kernel (JK) for N deg = 100 and N deg = 200.

However, because of the narrower distribution of the estimates when using the kernel we want to investigate the statistical behavior of the results obtained with a higher order of expansions, see Figure 14 for the results for N deg = 200. One can see that the narrowing of the R = 1 estimates by using the kernel is less significant than in the N deg = 100 case but still non-negligible especially in the low temperature regime. This can be seen in the deviation of mean from the exact result which is greater in the case without kernel. But again the deviation is still there for the same and even slightly higher temperatures. The best result is obtained with the R = 100 estimates without kernel. In this case an order of expansion of N deg = 100 is sufficient.

Figure 14:

The heat capacity C/k B of the Heisenberg ring with N = 10 spins s = 5/2 computed using the Chebyshev method (CM) and exact diagonalization (ED). Displayed are the results of P = 100 realizations with R = 1 and one with R = 100; N deg = 200. Compare Figure 12.

Also for the differential susceptibility (not shown), we obtain that the mean of the R = 1 estimates deviates from the exact diagonalization result without kernel more strongly than with kernel. The R = 100 estimate on the other hand is almost accurate without kernel but deviates as strongly as the mean in the case with kernel.

### 3.3 Sawtooth chain

The sawtooth chain (also known as delta chain) is an example with a highly degenerate spectrum. The Hamiltonian reads

(25) H = J 1 i = 1 N s i s i + 1 + J 2 k N / 2 s 2 k 1 s 2 k + 1

with periodic boundary conditions, ferromagnetic nearest neighbor interaction J 1 < 0 and antiferromagnetic next-nearest neighbor interaction J 2 > 0. We select a case with |J 2/J 1| = 0.45 which is close to the quantum critical point (QCP) at |J 2/J 1| = 1/2 [55]. The typicality approach has shown to be very efficient for these systems in schemes such as FTLM [44]. This can be confirmed for the Chebyshev method as well. In Figure 15, the R = 1 estimates of the heat capacity are distributed very narrowly around their mean which itself is perfectly aligned with R = 100. While the result without kernel is also aligned with the FTLM estimate, the result with kernel shows significant deviations from the FTLM result.

Figure 15:

The heat capacity C/k B of the sawtooth chain with N = 24 spins s = 1/2 computed using the Chebyshev method (CM) the finite temperature Lanczos method (FTLM). Displayed are the results of P = 100 realizations with R = 1 and one with R = 100.

We again show another result for a higher order of expansion N deg = 500, the lower graph in Figure 16. The deviation due to the kernel can be minimized, but still does not fall below those of the results without kernel.

Figure 16:

The heat capacity C/k B of the delta chain with N = 24 spins s = 1/2 computed using the Chebyshev method (CM) and the finite temperature Lanczos method (FTLM). Displayed are the results with and without Jackson kernel (JK) for N deg = 100 and N deg = 500.

## 4 Discussion and conclusions

We have seen that the Chebyshev method achieves very accurate results when handled with care. There are many possible choices for the parameters introduced for this method. We tried to identify some “good” choices and some methods to optimize them.

In particular we found that the number of points of integration N ̃ should be chosen closely to the order of expansion which itself has to be chosen according to the dimension of the problem. In the cases investigated, N deg = 100–200 is a sufficient choice as well as R ≥ 100. The parameter ɛ had no positive effect at least not when trying to approximate thermodynamic functions. So it seems advisable to set it equal to zero.

The most interesting “parameter” was whether to smooth the result with the Jackson kernel or not. We found here as well that for the investigated systems there was no positive effect. However, there might be a different use of the kernel. When the expansion of the density of states is completed the kernel can be employed without great computational effort to check the approximation with and without kernel for differences. For good results with a high order of the expansion the kernel changes the results only where they were already wrong, e.g., where heat capacity is negative. Therefore, if the kernel does not change the result too much one can be reasonably sure that the choice of the order of expansion is sufficient.

Finally, we can summarize that the approach via Chebyshev polynomials is accurate, but does not show any advantage compared to FTLM [44].

Corresponding author: Jürgen Schnack, Fakultät für Physik, Universität Bielefeld, Postfach 100131, Bielefeld D-33501, Germany, E-mail:

Award Identifier / Grant number: 355031190 (FOR 2692)

Award Identifier / Grant number: 397300368 (SCHN 615/25-1)

Award Identifier / Grant number: 449703145 (SCHN 615/28-1)

## Acknowledgments

Computing time at the Leibniz Supercomputing Center in Garching is gratefully acknowledged.

1. Author contribution: All the authors have accepted responsibility for the entire content of this submitted manuscript and approved submission.

2. Research funding: This work was supported by the Deutsche Forschungsgemeinschaft DFG (355031190 (FOR 2692); 397300368 (SCHN 615/25-1); 449703145 (SCHN 615/28-1)).

3. Conflict of interest statement: The authors declare no conflicts of interest regarding this article.

## References

[1] J. Skilling, “Maximum entropy and Bayesian methods,” in The Eigenvalues of Mega-Dimensional Matrices, Dordrecht, Kluwer, 1988, pp. 455–466.10.1007/978-94-015-7860-8_48Search in Google Scholar

[2] M. Hutchinson, “A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines,” Commun. Stat. Simulat. Comput., vol. 18, p. 1059, 1989. https://doi.org/10.1080/03610918908812806.Search in Google Scholar

[3] D. A. Drabold and O. F. Sankey, “Maximum entropy approach for linear scaling in the electronic structure problem,” Phys. Rev. Lett., vol. 70, p. 3631, 1993. https://doi.org/10.1103/physrevlett.70.3631.Search in Google Scholar

[4] J. Jaklič and P. Prelovšek, “Lanczos method for the calculation of finite-temperature quantities in correlated systems,” Phys. Rev. B, vol. 49, p. 5065, 1994.10.1103/PhysRevB.49.5065Search in Google Scholar PubMed

[5] R. N. Silver and H. Röder, “Densities of states of mega-dimensional Hamiltonian matrices,” Int. J. Mod. Phys. C, vol. 5, p. 735, 1994. https://doi.org/10.1142/s0129183194000842.Search in Google Scholar

[6] G. H. Golub and U. von Matt, “Tikhonov regularization for large scale problems,” Stanford University, Tech. Rep., Technical report SCCM-97-03, 1997.Search in Google Scholar

[7] A. Hams and H. De Raedt, “Fast algorithm for finding the eigenvalue distribution of very large matrices,” Phys. Rev. E, vol. 62, p. 4365, 2000. https://doi.org/10.1103/physreve.62.4365.Search in Google Scholar PubMed

[8] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, “The kernel polynomial method,” Rev. Mod. Phys., vol. 78, p. 275, 2006. https://doi.org/10.1103/revmodphys.78.275.Search in Google Scholar

[9] H. Avron and S. Toledo, “Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix,” J. ACM, vol. 58, no. 8, p. 1, 2011. https://doi.org/10.1145/1944345.1944349.Search in Google Scholar

[10] F. Roosta-Khorasani and U. Ascher, “Improved bounds on sample size for implicit matrix trace estimators,” Found. Comput. Math., vol. 15, p. 1187, 2015. https://doi.org/10.1007/s10208-014-9220-1.Search in Google Scholar

[11] A. K. Saibaba, A. Alexanderian, and I. C. F. Ipsen, “Randomized matrix-free trace and log-determinant estimators,” Numer. Math., vol. 137, p. 353, 2017. https://doi.org/10.1007/s00211-017-0880-z.Search in Google Scholar

[12] A. Wietek, P. Corboz, S. Wessel, B. Normand, F. Mila, and A. Honecker, “Thermodynamic properties of the Shastry-Sutherland model throughout the dimer-product phase,” Phys. Rev. Res., vol. 1, p. 033038, 2019. https://doi.org/10.1103/physrevresearch.1.033038.Search in Google Scholar

[13] K. Inoue, Y. Maeda, H. Nakano, and Y. Fukumoto, “Canonical-ensemble calculations of the magnetic susceptibility for a spin-1/2 spherical Kagome cluster with Dzyaloshinskii-Moriya interactions by using microcanonical thermal pure quantum states,” IEEE Trans. Magn., vol. 55, p. 1, 2019. https://doi.org/10.1109/tmag.2018.2873212.Search in Google Scholar

[14] S. Sugiura and A. Shimizu, “Thermal pure quantum states at finite temperature,” Phys. Rev. Lett., vol. 108, p. 240401, 2012. https://doi.org/10.1103/physrevlett.108.240401.Search in Google Scholar

[15] S. Sugiura and A. Shimizu, “Canonical thermal pure quantum state,” Phys. Rev. Lett., vol. 111, p. 010401, 2013. https://doi.org/10.1103/PhysRevLett.111.010401.Search in Google Scholar PubMed

[16] S. Okamoto, G. Alvarez, E. Dagotto, and T. Tohyama, “Accuracy of the microcanonical Lanczos method to compute real-frequency dynamical spectral functions of quantum models at finite temperatures,” Phys. Rev. E, vol. 97, p. 043308, 2018. https://doi.org/10.1103/PhysRevE.97.043308.Search in Google Scholar PubMed

[17] R. Alben, M. Blume, H. Krakauer, and L. Schwartz, “Exact results for a three-dimensional alloy with site diagonal disorder: comparison with the coherent potential approximation,” Phys. Rev. B, vol. 12, p. 4090, 1975. https://doi.org/10.1103/physrevb.12.4090.Search in Google Scholar

[18] H. De Raedt and P. de Vries, “Simulation of two and three-dimensional disordered systems: Lifshitz tails and localization properties,” Z. Phys. B Condens. Matter, vol. 77, p. 243, 1989. https://doi.org/10.1007/bf01313668.Search in Google Scholar

[19] P. de Vries and H. De Raedt, “Solution of the time-dependent Schrödinger equation for two-dimensional spin-1/2 Heisenberg systems,” Phys. Rev. B, vol. 47, p. 7929, 1993. https://doi.org/10.1103/physrevb.47.7929.Search in Google Scholar PubMed

[20] E. Dagotto, “Correlated electrons in high-temperature superconductors,” Rev. Mod. Phys., vol. 66, p. 763, 1994. https://doi.org/10.1103/revmodphys.66.763.Search in Google Scholar

[21] M. Aichhorn, M. Daghofer, H. G. Evertz, and W. von der Linden, “Low-temperature Lanczos method for strongly correlated systems,” Phys. Rev. B, vol. 67, no. R, p. 161103, 2003. https://doi.org/10.1103/physrevb.67.161103.Search in Google Scholar

[22] I. Zerec, B. Schmidt, and P. Thalmeier, “Kondo lattice model studied with the finite temperature Lanczos method,” Phys. Rev. B, vol. 73, p. 245108, 2006. https://doi.org/10.1103/physrevb.73.245108.Search in Google Scholar

[23] J. Schnack and O. Wendland, “Properties of highly frustrated magnetic molecules studied by the finite-temperature Lanczos method,” Eur. Phys. J. B, vol. 78, p. 535, 2010. https://doi.org/10.1140/epjb/e2010-10713-8.Search in Google Scholar

[24] J. Ummethum, J. Schnack, and A. Laeuchli, “Large-scale numerical investigations of the antiferromagnetic Heisenberg icosidodecahedron,” J. Magn. Magn Mater., vol. 327, p. 103, 2013. https://doi.org/10.1016/j.jmmm.2012.09.037.Search in Google Scholar

[25] O. Hanebaum and J. Schnack, “Advanced finite-temperature Lanczos method for anisotropic spin systems,” Eur. Phys. J. B, vol. 87, p. 194, 2014. https://doi.org/10.1140/epjb/e2014-50360-5.Search in Google Scholar

[26] C. Psaroudaki, J. Herbrych, J. Karadamoglou, P. Prelovšek, X. Zotos, and N. Papanicolaou, “Effective s=12$s=\frac{1}{2}$ description of the s = 1 chain with strong easy-plane anisotropy,” Phys. Rev. B, vol. 89, p. 224418, 2014. https://doi.org/10.1103/physrevb.89.224418.Search in Google Scholar

[27] R. Steinigeweg, J. Gemmer, and W. Brenig, “Spin and energy currents in integrable and nonintegrable spin-12$\frac{1}{2}$ chains: a typicality approach to real-time autocorrelations,” Phys. Rev. B, vol. 91, p. 104404, 2015. https://doi.org/10.1103/physrevb.91.104404.Search in Google Scholar

[28] R. Steinigeweg, J. Herbrych, F. Pollmann, and W. Brenig, “Typicality approach to the optical conductivity in thermal and many-body localized phases,” Phys. Rev. B, vol. 94, no. R, p. 180401, 2016. https://doi.org/10.1103/physrevb.94.180401.Search in Google Scholar

[29] Y. Yamaji, T. Suzuki, T. Yamada, S.-i. Suga, N. Kawashima, and M. Imada, “Clues and criteria for designing a Kitaev spin liquid revealed by thermal and spin excitations of the honeycomb iridate na2iro3,” Phys. Rev. B, vol. 93, p. 174425, 2016. https://doi.org/10.1103/physrevb.93.174425.Search in Google Scholar

[30] B. Schmidt and P. Thalmeier, “Frustrated two dimensional quantum magnets,” Phys. Rep., vol. 703, p. 1, 2017. https://doi.org/10.1016/j.physrep.2017.06.004.Search in Google Scholar

[31] J. Schnack, J. Schulenburg, and J. Richter, “Magnetism of the N = 42 Kagome lattice antiferromagnet,” Phys. Rev. B, vol. 98, p. 094423, 2018. https://doi.org/10.1103/physrevb.98.094423.Search in Google Scholar

[32] P. Prelovšek and J. Kokalj, “Finite-temperature properties of the extended Heisenberg model on a triangular lattice,” Phys. Rev. B, vol. 98, p. 035107, 2018.10.1103/PhysRevB.98.035107Search in Google Scholar

[33] I. Rousochatzakis, S. Kourtis, J. Knolle, R. Moessner, and N. B. Perkins, “Quantum spin liquid at finite temperature: proximate dynamics and persistent typicality,” Phys. Rev. B, vol. 100, p. 045117, 2019. https://doi.org/10.1103/physrevb.100.045117.Search in Google Scholar

[34] J. Schnack, J. Schulenburg, A. Honecker, and J. Richter, “Magnon crystallization in the Kagome lattice antiferromagnet,” Phys. Rev. Lett., vol. 125, p. 117207, 2020. https://doi.org/10.1103/physrevlett.125.117207.Search in Google Scholar

[35] P. Laurell and S. Okamoto, “Dynamical and thermal magnetic properties of the Kitaev spin liquid candidate α-RuCl3,” npj Quantum Mater., vol. 5, p. 2, 2020. https://doi.org/10.1038/s41535-019-0203-y.Search in Google Scholar

[36] K. Morita and T. Tohyama, “Finite-temperature properties of the Kitaev-Heisenberg models on Kagome and triangular lattices studied by improved finite-temperature Lanczos methods,” Phys. Rev. Res., vol. 2, p. 013205, 2020. https://doi.org/10.1103/physrevresearch.2.013205.Search in Google Scholar

[37] A. Honecker, J. Richter, J. Schnack, and A. Wietek, “Loop-gas description of the localized-magnon states on the Kagome lattice with open boundary conditions,” Condens. Matter Phys., vol. 23, p. 43712, 2020. https://doi.org/10.5488/cmp.23.43712.Search in Google Scholar

[38] U. Manthe and F. Huarte-Larranaga, “Partition functions for reaction rate calculations: statistical sampling and mctdh propagation,” Chem. Phys. Lett., vol. 349, p. 321, 2001. https://doi.org/10.1016/s0009-2614(01)01207-6.Search in Google Scholar

[39] F. Huarte-Larranaga and U. Manthe, “Vibrational excitation in the transition state: the CH4 + H → CH3 + H2 reaction rate constant in an extended temperature interval,” J. Chem. Phys., vol. 116, p. 2863, 2002. https://doi.org/10.1063/1.1436307.Search in Google Scholar

[40] J. Jaklič and P. Prelovšek, “Finite-temperature properties of doped antiferromagnets,” Adv. Phys., vol. 49, p. 1, 2000.10.1080/000187300243381Search in Google Scholar

[41] P. Prelovšek and J. Bonča, “Strongly correlated systems, numerical methods,” in Ground State and Finite Temperature Lanczos Methods, Berlin, Heidelberg, Springer, 2013.10.1007/978-3-642-35106-8_1Search in Google Scholar

[[42] E. Pavarini, E. Koch, R. Scalettar, and R. M. Martin, Eds. “The physics of correlated insulators, metals, and superconductors,” in The Finite Temperature Lanczos Method and its Applications by P. Prelovšek, 2017. http://hdl.handle.net/2128/15283 [ISBN 978-3-95806-224-5].Search in Google Scholar

[43] L. Lin, Y. Saad, and C. Yang, “Approximating spectral densities of large matrices,” 2013, ArXiv E-Prints, arXiv:1308.5467 [math.NA].Search in Google Scholar

[44] J. Schnack, J. Richter, and R. Steinigeweg, “Accuracy of the finite-temperature Lanczos method compared to simple typicality-based estimates,” Phys. Rev. Res., vol. 2, p. 013186, 2020. https://doi.org/10.1103/physrevresearch.2.013186.Search in Google Scholar

[45] J. Schnack, J. Richter, T. Heitmann, J. Richter, and R. Steinigeweg, “Finite-size scaling of typicality-based estimates,” Z. Naturforsch. A, vol. 75, p. 465, 2020. https://doi.org/10.1515/zna-2020-0031.Search in Google Scholar

[46] T. Iitaka and T. Ebisuzaki, “Algorithm for linear response functions at finite temperatures: application to ESR spectrum of s=12$s=\frac{1}{2}$ antiferromagnet Cu Benzoate,” Phys. Rev. Lett., vol. 90, p. 047203, 2003. https://doi.org/10.1103/PhysRevLett.90.047203.Search in Google Scholar PubMed

[47] A. Holzner, A. Weichselbaum, I. P. McCulloch, U. Schollwöck, and J. von Delft, “Chebyshev matrix product state approach for spectral functions,” Phys. Rev. B, vol. 83, p. 195115, 2011. https://doi.org/10.1103/physrevb.83.195115.Search in Google Scholar

[48] A. C. Tiegel, S. R. Manmana, T. Pruschke, and A. Honecker, “Matrix product state formulation of frequency-space dynamics at finite temperatures,” Phys. Rev. B, vol. 90, p. 060406, 2014. https://doi.org/10.1103/physrevb.90.060406.Search in Google Scholar

[49] A. C. Tiegel, S. R. Manmana, T. Pruschke, and A. Honecker, “Erratum: matrix product state formulation of frequency-space dynamics at finite temperatures [Phys. Rev. B 90, 060406(R) (2014)],” Phys. Rev. B, vol. 94, p. 179908, 2016. https://doi.org/10.1103/physrevb.94.179908.Search in Google Scholar

[50] J. L. Lado and O. Zilberberg, “Topological spin excitations in Harper-Heisenberg spin chains,” Phys. Rev. Res., vol. 1, p. 033009, 2019. https://doi.org/10.1103/physrevresearch.1.033009.Search in Google Scholar

[51] J. Mason and D. Handscomb, Chebyshev Polynomials, Boca Raton, Florida, CRC Press, 2002.10.1201/9781420036114Search in Google Scholar

[52] I. Bronstein, K. Semendjaew, G. Grosche, V. Ziegler, and D. Ziegler, Springer-Handbook of Mathematics, E. Zeidler, Ed., Leipzig, Springer, 2007.Search in Google Scholar

[53] E. Dagotto and T. M. Rice, “Surprises on the way from one- to two-dimensional quantum magnets: the ladder materials,” Science, vol. 271, p. 618, 1996. https://doi.org/10.1126/science.271.5249.618.Search in Google Scholar

[54] H.-J. Mikeska and A. K. Kolezhuk, “One-dimensional magnetism,” in Quantum Magnetism, U. Schollwöck, J. Richter, D. J. J. Farnell, and R. F. Bishop, Eds., Berlin, Heidelberg, Springer, 2004, pp. 1–83.10.1007/BFb0119591Search in Google Scholar

[55] V. Y. Krivnov, D. V. Dmitriev, S. Nishimoto, S.-L. Drechsler, and J. Richter, “Delta chain with ferromagnetic and antiferromagnetic interactions at the critical point,” Phys. Rev. B, vol. 90, p. 014441, 2014. https://doi.org/10.1103/physrevb.90.014441.Search in Google Scholar

[56] D. V. Dmitriev, V. Y. Krivnov, J. Richter, and J. Schnack, “Thermodynamics of a delta chain with ferromagnetic and antiferromagnetic interactions,” Phys. Rev. B, vol. 99, p. 094410, 2019. https://doi.org/10.1103/physrevb.99.094410.Search in Google Scholar

[57] A. Baniodeh, N. Magnani, Y. Lan, et al.., “High spin cycles: topping the spin record for a single molecule verging on quantum criticality,” npj Quantum Mater., vol. 3, p. 10, 2018. https://doi.org/10.1038/s41535-018-0082-7.Search in Google Scholar

Revised: 2021-05-31
Accepted: 2021-06-01
Published Online: 2021-06-21
Published in Print: 2021-09-27