# Accuracy of the typicality approach using Chebyshev polynomials

H. Schlüter, F. Gayk, H. -J. Schmidt, A. Honecker, J. Schnack

► To cite this version:

H. Schlüter, F. Gayk, H. -J. Schmidt, A. Honecker, J. Schnack. Accuracy of the typicality approach using Chebyshev polynomials. Zeitschrift für Naturforschung A, 2021, 76 (9), pp.823-834. 10.1515/zna-2021-0116. hal-03210736

HAL Id: hal-03210736

https://hal.science/hal-03210736v1

Submitted on 3 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

# Accuracy of the typicality approach using Chebyshev polynomials

https://doi.org/10.1515/zna-2021-0116

Received April 28, 2021; accepted June 1, 2021;

published online June 21, 2021

**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.

**Keywords:** Chebyshev expansion; spin systems; trace estimators; typicality.

## 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–12]. These schemes, sometimes also called *typicality* or (*microcanonical*) *thermal pure quantum states* [13–16], have been used very successfully, in particular, in the field of correlated electron systems, see e.g., [16–37] but also in quantum chemistry [38, 39].

A prominent formulation of the method is the finite temperature Lanczos method (FTLM) [4, 21, 40–42] which employs a Krylov space expansion for $exp\{-\beta 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–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–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 + \epsilon/2, 1 - \epsilon/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

<figure><img src="image_3.png" /></figure>

<footnote>*Corresponding author: Jürgen Schnack, Fakultät für Physik, Universität Bielefeld, Postfach 100131, Bielefeld D-33501, Germany, E-mail: jschnack@uni-bielefeld.de. https://orcid.org/0000-0003-0702-2723</footnote>

<footnote>Henrik Schlüter and Florian Gayk, Fakultät für Physik, Universität Bielefeld, Postfach 100131, Bielefeld D-33501, Germany, E-mail: hschlueter@physik.uni-bielefeld.de (H. Schlüter), fgayk@physik.uni-bielefeld.de (F. Gayk)</footnote>

<footnote>Heinz-Jürgen Schmidt, Fachbereich Physik, Universität Osnabrück, Osnabrück D-49069, Germany, E-mail: hschmidt@uos.de. https://orcid.org/0000-0001-8087-3204</footnote>

<footnote>Andreas Honecker, Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089, CY Cergy Paris Université, Cergy-Pontoise Cedex F-95302, France, E-mail: andreas.honecker@cyu.fr. https://orcid.org/0000-0001-6383-3200</footnote>
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

$$ \rho(E) := \sum_n \delta(E - E_n). \qquad (1) $$

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

$$ Z(\beta) = \sum_n e^{-\beta E_n} = \int_{-\infty}^{\infty} \rho(E) e^{-\beta E} dE \qquad (2) $$

with $\beta = \frac{1}{k_B T}$. Correspondingly, the heat capacity is evaluated as

$$ \frac{C(\beta)}{k_B} = \beta^2 \left[ \frac{1}{Z(\beta)} \int_{-\infty}^{\infty} \rho(E) e^{-\beta E} E^2 dE - \left( \frac{1}{Z(\beta)} \int_{-\infty}^{\infty} \rho(E) e^{-\beta E} EdE \right)^2 \right]. \quad (3) $$

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.,

$$ \frac{\chi(\beta)}{(g\mu_B)^2} = \beta \left[ \frac{1}{Z(\beta)} \sum_M M^2 \int_{-\infty}^{\infty} \rho(E, M) e^{-\beta E} dE - \left( \frac{1}{Z(\beta)} \sum_M M \int_{-\infty}^{\infty} \rho(E, M) e^{-\beta E} dE \right)^2 \right]. \quad (4) $$

We further calculate only contributions for $M \ge 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 $\rho(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

$$ H' := \frac{1}{m}(H - c \cdot 1) \qquad (5) $$

with

$$ m = \frac{E_{\max} - E_{\min}}{2 - \varepsilon}, \quad c = \frac{E_{\max} + E_{\min}}{2} \qquad (6) $$

where, as suggested in [8], a parameter $\varepsilon$ 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 + \frac{\varepsilon}{2}, 1 - \frac{\varepsilon}{2}]$.

The corresponding scaled density of states $\bar{\rho}(x)$ is then expanded in terms of Chebyshev polynomials $C_n(x)$:

$$ \bar{\rho}(x) \approx \frac{1}{\pi \sqrt{1-x^2}} \left[ \mu_0 + 2 \sum_{n=1}^{N_{\deg}} \mu_n C_n(x) \right]. \qquad (7) $$

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

$$ \mu_n = \mathrm{Tr} [C_n(H')] \qquad (8) $$

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

$$ \mu_n \approx \frac{\dim(H)}{R} \sum_{r=1}^{R} \frac{\langle r | C_n(H') | r \rangle}{\langle r | r \rangle} = \Theta_n(R), \qquad (9) $$

with

$$ |r\rangle = \sum_{\nu} r_{\nu} |\nu\rangle \qquad (10) $$

being a random vector with Gaussian distributed components $r_\nu$ with respect to a chosen orthonormal basis $\{|\nu\rangle\}$. The relative error of an estimate $\Theta_n(R)$ is proportional to $1/\sqrt{R \dim(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 $\mu_n$ by a kernel [8]. The kernel fixes this problem at the cost of introducing a systematic error which vanishes for $N_{\deg} \to \infty$. In this paper, we restrict our discussion to the use of the Jackson kernel whose coefficients read
$$g_n = \frac{1}{\tilde{N}_{\text{deg}} + 1} \left[ (N_{\text{deg}} - n + 1) \cos\left(\frac{\pi n}{N_{\text{deg}} + 1}\right) + \sin\left(\frac{\pi n}{N_{\text{deg}} + 1}\right) \cot\left(\frac{\pi}{N_{\text{deg}} + 1} N_{\text{deg}} + 1\right) \right]. \quad (11)$$

In figure captions or legends, we will write $g_n = \text{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

$$\int_{-1}^{1} \frac{f(x)}{\sqrt{1-x^2}} dx \approx \frac{\pi}{\tilde{N}} \sum_{k=1}^{\tilde{N}} f(x_k), \quad \text{see [52],} \quad (12)$$

where the supporting points read

$$x_k = \cos\left(\frac{\pi(k - \frac{1}{2})}{\tilde{N}}\right), \quad k = 1, \dots, \tilde{N}. \quad (13)$$

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

$$f(x) := \left[ g_0 \mu_0 + 2 \sum_{n=1}^{N_{\text{deg}}} g_n \mu_n C_n(x) \right] e^{\beta(mx+c)}, \quad (14)$$

and is thus no polynomial of order smaller than $2\tilde{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

$$Z = \int_{-\infty}^{\infty} \rho(E) e^{-\beta E} dE \quad (15)$$

$$= \int_{-1}^{1} \bar{\rho}(x) e^{\beta(mx+c)} dx \quad (16)$$

$$\approx \frac{1}{\tilde{N}} \sum_{k=1}^{\tilde{N}} \gamma_k e^{-\beta(mx_k+c)}, \quad (17)$$

where the values of the weights $\gamma_k$ read

$$\gamma_k := \pi \sqrt{1 - x_k^2} \bar{\rho}_M(x_k) \quad (18)$$

$$= g_0 \mu_0 + 2 \sum_{n=1}^{N_{\text{deg}}} g_n \mu_n \cos\left(\frac{n\pi(k - \frac{1}{2})}{\tilde{N}}\right). \quad (19)$$

If one chooses $\tilde{N} \ge N_{\text{deg}}$, the sum can be complemented to an upper limit of $\tilde{N}$ with additional $g_n = 0$ terms. The $\gamma_k$ can then be computed through a discrete cosine-transform (type III) of the coefficients $g_n \mu_n$ which allows

a faster computation of the sum. The time needed scales with $\tilde{N} \ln \tilde{N}$ instead of $\tilde{N} N_{\text{deg}}$ [8].

If there are known symmetries, the scheme can be performed for each orthogonal subspace $\mathcal{H}_\Gamma$ separately. The approximated partition function can then be written as

$$Z(\beta) = \frac{1}{\tilde{N}} \sum_{\Gamma} \sum_{k=1}^{\tilde{N}} \gamma_k^{\Gamma} e^{-\beta(m_{\Gamma}x_k + c_{\Gamma})}. \quad (20)$$

All of the following systems possess $\tilde{S}^z$ symmetry which implies orthonormal subspaces corresponding to the quantum number of the total magnetization $\Gamma = 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_{\text{deg}}$, the scaling parameter $\epsilon$, the number of random vectors $R$, the use of a kernel $g_n$, and the number of supporting points $\tilde{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(\beta)$ as a random measurement, a mean $\bar{O}(\beta)$ and a variance $\delta O(\beta)^2$ can be defined:

$$\bar{O}(\beta) = \frac{1}{P} \sum_{i=0}^{P-1} O_i(\beta), \quad (21)$$

$$\delta O(\beta)^2 = \bar{O}^2(\beta) - \bar{O}(\beta)^2. \quad (22)$$

If additionally an exact result $O_E$ is known, the systematic deviation

$$\Delta O(\beta) = |\bar{O}(\beta) - O_E(\beta)| \quad (23)$$

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–57]; compare Figure 1.

## 3.1 Heisenberg ladder

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

$$ \begin{aligned} \hat{H} ={}& J_1 \sum_i \vec{s}_{\sim i,1} \cdot \vec{s}_{\sim i,2} + J_2 \sum_{i,j} \vec{s}_{\sim i,j} \cdot \vec{s}_{\sim i+1,j} \\ & + g \mu_B B \sum_{i,j} s_{\sim i,j}^z . \end{aligned} \quad (24) $$

where the first subscript $i \in \{1, \dots, N/2\}$ of the spin operators denotes the rung and the second subscript $j \in \{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_{\text{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 $\epsilon$ 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.

<figure><img src="image_4.png" /><figcaption>**Figure 1:** Systems investigated in Section 3.1, from top to bottom: Ladder, chain, sawtooth chain. Periodic boundary conditions will be applied.</figcaption></figure>

<table><caption>**Table 1:** Standard configuration of parameters used in the following calculations.</caption><thead><tr><th>R</th><th>N<sub>deg</sub></th><th>ε</th><th>N&#x0304;</th><th>g<sub>n</sub></th></tr></thead><tbody><tr><td>200</td><td>100</td><td>0</td><td>100</td><td>1</td></tr></tbody></table>

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_{\text{deg}}$. The equality is chosen as a starting point.

### 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

<figure><img src="image_5.png" /><figcaption>**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 $\mu_n$ is done using a complete basis.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**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 $\mu_n$ is done using a complete basis.</figcaption></figure>

result determined by exact diagonalization. In addition, a result determined by the Chebyshev algorithm is shown where the traces $\mu_n$ are computed numerically exactly instead of approximating them using random vectors.

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| \approx 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}$, $\epsilon$, $\tilde{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.

(a) The order of expansion can be increased to push the systematic deviations, i.e., the “ghost dip”, to

<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**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 $\epsilon$ (colored curves) compared to exact diagonalization (ED).</figcaption></figure>
lower temperatures. This is demonstrated in Figure 4. Computation time increases linearly with $N_{deg}$.

(b) The scaling parameter $\varepsilon$ 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 $\varepsilon = 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 $\varepsilon$ has the intended effect to prevent the peaks of the lowest and highest eigenvalues to be cut off. There are also situations where a $\varepsilon \approx 10^{-6}$ seemingly decreases the depth of the dip compared to $\varepsilon = 0$. Anyhow, such small values of $\varepsilon$ are not sufficient to prevent the cut-off of the density of states, and the improvement was not significant when compared to statistical deviations.

<figure><img src="image_9.png" /><figcaption>**Figure 6:** Detail of the scaled density of states $\rho_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 $\varepsilon$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**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 $\tilde{N}$ (colored) and with exact diagonalization (ED).</figcaption></figure>

(c) The number of supporting points $\tilde{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

<figure><img src="image_11.png" /><figcaption>**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).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**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$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**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$.</figcaption></figure>

results, but it is important to note that a good choice of $\tilde{N}$ scales with $N_{deg}$.

(d) 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.

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$.

## 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.

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

<figure><img src="image_14.png" /><figcaption>**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$.</figcaption></figure>
$k_B T / |J| = 2 \times 10^{-2}$ and $2 \times 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.

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 Fe$_{10}$) 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.

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_{\text{deg}}$, see Figure 13 with a linear temperature axis. Note that even for $N_{\text{deg}} = 200$ the result without kernel outperforms the one with kernel. Furthermore, the result with $N_{\text{deg}} = 100$ without kernel is more accurate than the result with $N_{\text{deg}} = 200$ and kernel.

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_{\text{deg}} = 200$. One can see that the narrowing of the $R = 1$ estimates by using the kernel is less significant than in the $N_{\text{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

<figure><img src="image_15.png" /><figcaption>**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_{\text{deg}} = 100$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**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_{\text{deg}} = 100$ and $N_{\text{deg}} = 200$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**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$.</figcaption></figure>

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.

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

$$ \sim H = J_1 \sum_{i=1}^{N} \vec{s}_i \cdot \vec{s}_{i+1} + J_2 \sum_{k}^{N/2} \vec{s}_{2k-1} \cdot \vec{s}_{2k+1} \quad (25) $$

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.

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.

## 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
<figure><img src="image_19.png" /><figcaption>**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 ($J$) for $N_{deg} = 100$ and $N_{deg} = 500$.</figcaption></figure>

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 $\tilde{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 \ge 100$. The parameter $\varepsilon$ 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].

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

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

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

**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.

[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.

[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.

[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.

[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.

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

[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.

[8] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, “The kernel polynomial method,” *Rev. Mod. Phys.*, vol. 78, p. 275, 2006.

[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.

[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.

[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.

[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.

[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.

[14] S. Sugiura and A. Shimizu, "Thermal pure quantum states at finite temperature," *Phys. Rev. Lett.*, vol. 108, p. 240401, 2012.

[15] S. Sugiura and A. Shimizu, "Canonical thermal pure quantum state," *Phys. Rev. Lett.*, vol. 111, p. 010401, 2013.

[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.

[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.

[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.

[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.

[20] E. Dagotto, "Correlated electrons in high-temperature superconductors," *Rev. Mod. Phys.*, vol. 66, p. 763, 1994.

[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.

[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.

[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.

[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.

[25] O. Hanebaum and J. Schnack, "Advanced finite-temperature Lanczos method for anisotropic spin systems," *Eur. Phys. J. B*, vol. 87, p. 194, 2014.

[26] C. Psaroudaki, J. Herbrych, J. Karadamoglou, P. Prelovšek, X. Zotos, and N. Papanicolaou, "Effective $s = \frac{1}{2}$ description of the $s = 1$ chain with strong easy-plane anisotropy," *Phys. Rev. B*, vol. 89, p. 224418, 2014.

[27] R. Steinigeweg, J. Gemmer, and W. Brenig, "Spin and energy currents in integrable and nonintegrable spin-$\frac{1}{2}$ chains: a typicality approach to real-time autocorrelations," *Phys. Rev. B*, vol. 91, p. 104404, 2015.

[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.

[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 $na_2iro_3$," *Phys. Rev. B*, vol. 93, p. 174425, 2016.

[30] B. Schmidt and P. Thalmeier, "Frustrated two dimensional quantum magnets," *Phys. Rep.*, vol. 703, p. 1, 2017.

[31] J. Schnack, J. Schulenburg, and J. Richter, "Magnetism of the $N=42$ Kagome lattice antiferromagnet," *Phys. Rev. B*, vol. 98, p. 094423, 2018.

[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.

[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.

[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.

[35] P. Laurell and S. Okamoto, "Dynamical and thermal magnetic properties of the Kitaev spin liquid candidate $\alpha$-RuCl$_3$," *npj Quantum Mater.*, vol. 5, p. 2, 2020.

[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.

[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.

[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.

[39] F. Huarte-Larranaga and U. Manthe, "Vibrational excitation in the transition state: the $CH_4 + H \rightarrow CH_3 + H_2$ reaction rate constant in an extended temperature interval," *J. Chem. Phys.*, vol. 116, p. 2863, 2002.

[40] J. Jaklíč and P. Prelovšek, "Finite-temperature properties of doped antiferromagnets," *Adv. Phys.*, vol. 49, p. 1, 2000.

[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.

[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].

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

[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.

[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.
[46] T. Iitaka and T. Ebisuzaki, "Algorithm for linear response functions at finite temperatures: application to ESR spectrum of $s = \frac{1}{2}$ antiferromagnet Cu Benzoate," *Phys. Rev. Lett.*, vol. 90, p. 047203, 2003.

[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.

[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.

[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.

[50] J. L. Lado and O. Zilberberg, "Topological spin excitations in Harper-Heisenberg spin chains," *Phys. Rev. Res.*, vol. 1, p. 033009, 2019.

[51] J. Mason and D. Handscomb, *Chebyshev Polynomials*, Boca Raton, Florida, CRC Press, 2002.

[52] I. Bronstein, K. Semendjaew, G. Grosche, V. Ziegler, and

D. Ziegler, *Springer-Handbook of Mathematics*, E. Zeidler, Ed., Leipzig, Springer, 2007.

[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.

[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.

[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.

[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.

[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.