Real-space spectral simulation of quantum spin models: Application to generalized Kitaev models

The proliferation of quantum fluctuations and long-range entanglement presents an outstanding challenge for the numerical simulation of interacting spin systems with exotic ground states. Here, we present a toolset of Chebyshev polynomial-based iterative methods that provides a unified framework to study the thermodynamical properties, critical behavior and dynamics of frustrated quantum spin models with controlled accuracy. Similar to previous applications of the Chebyshev spectral methods to condensed matter systems, the algorithmic complexity scales linearly with the Hilbert space dimension and the Chebyshev truncation order. Using this approach, we study two paradigmatic quantum spin models on the honeycomb lattice: the Kitaev-Heisenberg (K-H) and the Kitaev-Ising (K-I) models. We start by applying the Chebyshev toolset to compute nearest-neighbor spin correlations, specific heat and entropy of the K-H model on a 24-spin cluster. Our results are benchmarked against exact diagonalization and a popular iterative method based on thermal pure quantum states. The transitions between a variety of magnetic phases, namely ferromagnetic, N\'eel, zigzag and stripy antiferromagnetic and quantum spin liquid phases are obtained accurately and efficiently. We also accurately obtain the temperature dependence of the spin correlations, over more than three decades in temperature, by means of a finite temperature Chebyshev polynomial method introduced here. Finally, we report novel dynamical signatures of the quantum phase transitions in the K-I model. Our findings suggest that the efficiency, versatility and low-temperature stability of the Chebyshev framework developed here could pave the way for previously unattainable studies of quantum spin models in two dimensions.


Introduction
In strongly correlated materials, the interplay between different types of interactions can engender rich T = 0 phase diagrams characterized by a series of transitions between paramagnetic, magnetically ordered and quantum spin liquid phases.These quantum phase transitions are driven by one or more parameters in the Hamiltonian, such as an external magnetic field or a spin exchange coupling constant.When a drastic change in the ground state occurs as one of these parameters is varied, there is an accompanying change in the thermodynamic properties.This change manifests itself in the form of critical behavior of quantities such as the structure factor and the susceptibility, which scale with the parameter that drives the transition [1,2].
The study of the competition between quantum fluctuations and interactions at the heart of quantum phase transitions often calls for a numerical approach [3][4][5][6][7][8][9][10][11].Exact solutions are only known for a handful of cases (a well-known example is the isotropic Heisenberg chain [12]).Moreover, in low dimensions, the presence of strong quantum fluctuations limits the applicability of mean-field approaches.Exact diagonalization (ED) methods, such as those based on the Lanczos algorithm [13], are the next step beyond exact analytical solutions.There are several examples of the application of this method to interacting spin systems, for example Refs.[3,11,14].Unfortunately, these are limited to relatively small system sizes, even when algorithms are optimized to reflect symmetries of the model.The culprit is the exponential scal-ing of the computational cost with the system size, which is particularly severe in dimensions greater than one.Beyond the Lanczos algorithm and its more recent variants [15][16][17], several attempts to go beyond the limitations of full exact diagonalization have been made.Potent numerical techniques have been deployed with varying degrees of success, including series expansions [18][19][20][21][22][23][24], quantum Monte Carlo (QMC) [7,8,25,26], density matrix renormalization group (DMRG) [27][28][29][30], tensor-network approaches (such as iPEPS) [10,[31][32][33] and thermal pure quantum (TPQ) states [32,[34][35][36].Efficient numerical schemes amenable to large-scale computations share a key feature: they aim at reconstructing expectation values of quantum observables without having to fully diagonalize the Hamiltonian.The resulting computational cost depends crucially on how the expectation values of the observables are evaluated.Here, two relevant aspects are at play.The first has to do with how the corresponding operators are reconstructed.The second relates to the process by which one obtains the expectation value.Usually this process is a stochastic one, unless there is prior knowledge about some of the system's features, in which case a variational approach can be viable [25,37,38].
In principle, QMC methods can be used to probe large systems in any number of dimensions, while remaining numerically exact [7,8].However, they often encounter the so-called sign problem [26,39].This is a situation where the variance of the estimators of quantities of interest increases exponentially due to quantum statistics.The severity of the problem depends on the computational basis used to tackle the specific model [40][41][42][43][44][45][46].Generally, the sign problem tends to be more acute in frustrated systems [47,48], hampering the use of QMC to extract quantities of interest, such as correlation functions.The sign problem and the limited range of models that QMC is able to access emphasize the need for a general purpose method that can be used more broadly as an alternative to Lanczos ED and QMC.
In this work, we apply a toolset of spectral methods to compute both static (e.g.spin correlations) and dynamic (e.g.spin susceptibility) quantum observables in paradigmatic frustrated systems with competing interactions.Spectral methods based upon Chebyshev expansions have recently proven useful in different contexts [49][50][51][52][53][54][55][56][57][58][59][60], and here we will be interested in extending this approach to models of interacting quantum spins.We focus on generalized honeycomb Kitaev models, i.e. systems that combine Kitaev interactions with other types of magnetic exchange.Specifically, we study the Kitaev-Heisenberg (K-H) model and the Kitaev-Ising (K-I) model [3][4][5]61].The Kitaev model on the honeycomb lattice is one of the rare examples of an exactly solvable microscopic model [62] showing exchange frustration, i.e. nearest neighbor (NN) interactions that cannot be simultaneously minimized.This is similar to geometric frustration which notably occurs in the case of the antiferromagnetic Ising model on the triangular lattice [63].In the Kitaev case, frustration is created by bond-directional interactions which give way to a fractionalized excitation spectrum of Majorana fermions.This has attracted great attention because it opens up the possibility of synthesizing spin liquid materials with exotic topological orders [64,65].Notably, in honeycomb iridates, which are transition metal oxides with partially filled d-shells, a subtle interplay of spin-orbit coupling and electronic correlations produces the type of bond-directional interactions that appear in the Kitaev model.Thus, these materials make good spin liquid candidates.In these spin-orbit assisted Mott insulators, the Kitaev exchange interaction is thought to be responsible for the emergence of a spin liquid phase [3,[65][66][67][68].
The Kitaev exchange interaction also plays a key role in the modelling of other compounds such as the van der Waals ruthenate α-RuCl 3 .A recent study proposes a minimal microscopic 2D spin model for α-RuCl 3 [69].The model at play is an extension of the Kitaev model that considers both Kitaev and Heisenberg exchange interactions, third neighbor exchange and Γ interactions (i.e., terms that couple different spin components for each nearest neighbor bond).The authors treat this generalized K-H model using a a mean-field random-phase approximation, aiming at extracting quantities such as the dynamical structure factor [69].The use of this mean-field approach is justified by comparing its results with those of exact diagonalization [14].However, exact diagonalization is limited to relatively small system sizes (for example, in Ref. [14], a 24-site cluster is used).Thus, this approach runs into the risk of overlooking large scale properties of the model.In fact, in Refs.[70,71], the authors detect finite-size effects when computing the dynamical spin structure factor using QMC simulations of the Kitaev model in the presence of disorder -which deems the study of large scale properties crucial -even for clusters of 288 sites. 1 These developments illustrate the need to develop accurate, general purpose computational methods that scale favorably with the system size.
The Chebyshev polynomial methods we shall introduce below have a comparable complexity to Lanczos-based methods.They both share the advantage of being free of the sign problem, even though they can't reach the same system sizes as QMC.Yet, Chebyshev expansions have a few key features that can make them more advantageous than their Lanczos counterparts, namely superior robustness and accuracy.For example, the finite temperature Lanczos method [15] has a low-temperature counterpart [16] that was developed to tackle loss of accuracy due to statistical convergence issues at low-temperature.In this work, we develop a seamless Chebyshev approach that is accurate and does not require a low-temperature counterpart.Moreover, as we will show below, it has advantages even over TPQ [34,35].Another application is to study dynamics (e.g., spectral functions), where it proves more flexible and efficient than its Lanczos-based counterpart [17].
This paper is organized as follows.Section 2 provides a bird's-eye view of the Chebyshev scheme in condensed matter physics.In Sec. 3 we give details of the techniques used throughout this study: i) in the microcanonical ensemble, we use the microcanonical Lanczos (MCLM) method (Sec.3.2.1), the microcanonical (MTPQ) variant of the TPQ (Sec.3.2.2),and the iterative Chebyshev polynomial Green's function (CPGF) method (Sec.3.2.3);ii) in the canonical ensemble, we use the finite temperature Lanczos (FTLM) method (Sec.3.3.1),the canonical (CTPQ) variant of the TPQ (Sec.3.3.2),and the newly developed finite temperature Chebyshev polynomial (FTCP) method that we introduce in this work (Sec.3.3.3);iii) for dynamical studies, we use the Lanczos method using the continued fraction approach and a hybrid Lanczos-Chebyshev approach, also introduced in this work (Sec.3.4.2).Section 4 compares the convergence properties and performance of all these methods in detail by applying them to study the K-H and K-I models on the honeycomb lattice: i) we start by checking consistency at zero temperature by computing the ground state energy and nearest-neighbor spin-spin correlation function of the K-H model with all methods and comparing them with previously known results that we reproduced using the ED Lanczos technique; ii) for the spinspin correlation, we present a detailed analysis of the dependence of our estimators on the relevant parameters: the number of initial random states, truncation order and, in the case of the CPGF, also the energy resolution; iii) we study the temperature dependence of the nearestneighbor spin-spin correlation, specific heat and entropy of the K-H model and show that our results match the tensor network results of Ref. [10]; iv) we bench-mark our implementation of the hybrid Lanczos-Chebyshev method and compute the dynamical spin susceptibility for the K-I model, finding dynamical signatures of the quantum phase transitions described in Ref. [61].Finally, in Sec. 5, we point out the pros and cons of each method.In particular, we summarize how this work highlights the efficiency of the CPGF, FTCP and the hybrid Lanczos-Chebyshev methods and suggest potentially interesting applications of these methods.We also discuss the dynamical signatures of the quantum phase transitions in the K-I model that we found using the hybrid Lanczos-Chebyshev approach.Figure 1: Microcanonical CPGF approach, the rescaled energy spectrum is probed using a coarse-grained average of energy states within a specified energy resolution η.

Chebyshev spectral methods: Rationale
Spectral methods are an increasingly popular tool for the simulation of condensed matter systems that fulfill the requirement of general applicability [49,52,[56][57][58][59][60][72][73][74][75][76][77].These methods rely on the iterative reconstruction of the target functions (e.g.static or dynamic correlation functions), generally in terms of Chebyshev polynomial expansions due to their favorable convergence properties [78].The iterative scheme is stable and can be made as accurate as required within a specified parameter.For example, if one is interested in computing expectations of quantum observables in the microcanonical ensemble, this parameter is the energy resolution.The canonical ensemble analogue of this parameter is the temperature.Let us take the instructive case of the microcanonical ensemble.The spectral approach uses a coarse-grained description of energy states to provide estimates for quantum observables in large systems (see Fig. 1).This is to be contrasted with exact diagonalization, which relies on the knowledge of individual states and thus is limited to small systems.Furthermore, spectral methods can be combined with stochastic techniques for computation of traces to further reduce the computational cost and are amenable to parallelization as we will see briefly.
The CPGF approach we exploit in this work to compute microcanonical averages [51,79] has proven effective in dealing with tight-binding models, allowing unparalleled large-scale simulations with billions of atomic orbitals [51,55].Motivated by these developments, the main aim of this work is to introduce a finite-temperature spectral framework that can capture the physics of two-dimensional quantum spin models over a wide range of temperatures (of particular interest will be to probe the low-temperature behavior of spin liquids).
Spectral methods leverage efficient stochastic estimators for expectation values that use random vectors 2 to evaluate traces of operators [80].This technique, dubbed stochastic trace evaluation (STE), is ubiquitous in the study of condensed phases and is used in ED methods, such as those based on the Lanczos algorithm and TPQ states [34,35], and in the kernel polynomial method [49].The rationale in the STE is to approximate the trace of an operator by an average of expectation values using N rd.vec.random vectors, |φ The relative error scales favorably with the Hilbert space dimension, D (in fact the relative 2 A random vector is defined as |φ 0 〉 = D i=1 ξ i |i〉, with {|i〉} an arbitrary basis and ξ i ∈ random variables that satisfy ξ i = 0, ξ i ξ j = δ i j and ξ * i ξ j = δ i j (here the bar denotes statistical average).The r-th realization, |φ 0 〉 corresponds to a specific set of coefficients, {ξ r i , i = 1, 2, . . ., D.} and different sets are assumed to be uncorrelated.
error is proportional to 1/ D for typical sparse operators) and, for a fixed system size, can be made as small as desired by increasing N rd.vec. .Moreover, stochastic trace estimators are free from the sign problem.
A crucial feature of Chebyshev expansions is that they offer uniform convergence [78] (and in the case of CPGF this translates into an energy resolution that can be specified exactly [51]).This is appealing for studies of phase transitions, particularly when one wishes to characterize the critical behavior of thermodynamic functions, among others.There are two ways to define a resolution.One of them uses a kernel that modifies the coefficients of the Chebyshev expansion.This modification smears out so called Gibbs oscillations which occur upon truncation of orthogonal polynomial series [49].A resolution may then be defined as the spread of the kernel in the x y-plane, 3 and generally depends on the truncation order and energy.
Here, instead, we use a Green's function-based method that was proposed independently in the works of Ferreira and Mucciolo [51] and Braun and Schmitteckert [79].This approach, coined CPGF [51], has two main features: (i) it is based on a stable, asymptotically exact expansion of lattice Green's functions in Chebyshev polynomials; and (ii) the energy resolution is specified from the outset, in the form of a simple imaginary self-energy.
By extending these ideas to quantum spin models, we propose a Chebyshev-based method for the computation of quantum expectations in the canonical ensemble, where the temperature plays the role of a resolution.This method -which we shall detail below -bypasses potential low-temperature convergence issues by using an adaptive temperature step, while maintaining rigorous control over convergence.Alternatively, TPQ-based methods can be used to approximate either microcanonical [34] or canonical [35] averages by successive application of the Hamiltonian operator onto an initial random state.In TPQ, the number of iterations is proportional to a quantity that plays the role of an effective temperature.Broadly speaking, this effective temperature acts similarly to a resolution that becomes finer as more iterations are completed.Yet, this annealing scheme is susceptible to slow convergence, particularly in the vicinity of critical points as shown later.
In the following sections, we will compare Lanczos, TPQ, and Chebyshev-based approaches since they all scale linearly with the dimension of the Hilbert space D.Moreover, all methods scale linearly with the number of polynomials required for spectral convergence (or iterations in the case of TPQ) N poly/it.and with the number of realizations of the initial random state required for statistical convergence N rd.vec. .

Methodology
The methods described here involve two main steps.First, a numerically exact or approximate spectral representation of the target state is obtained by means of algorithms with polynomial computational complexity.This is achieved by recursive application of the Hamiltonian, Ĥ, to an initial random state, |φ 0 〉.Some examples of typical target states are the ground state, a microcanonical state (with energy restricted to an energy shell) or a canonical (finite temperature) state.After the target state is converged to the desired precision, physical observables may be computed by means of the STE technique introduced earlier in Sec. 2, i.e. by averaging the expectation value 〈φ 0 | Ô|φ 0 〉 over an ensemble of random vectors.
These techniques are general in scope and, as discussed below, when combined provide a powerful means of accessing excited states, reconstructing Green's functions, computing average and local density of states, and can be easily extended to the study of quantum dynamics, either in the time domain, by exploiting a spectral approximation of the time evolution operator, or in the frequency domain, via the resolvent operator.

Lanczos exact diagonalization
While the iterative methods discussed here generate different polynomials of the Hamiltonian during the recursion, they have a crucial aspect in common.After M iterations, they create a state in the so called Krylov subspace, defined as follows: K M ( ĥ, |φ 0 〉) ≡ span{|φ 0 〉 , ĥ |φ 0 〉 , ĥ2 |φ 0 〉 , . . .ĥM |φ 0 〉} , where ĥ = Ĥ/N is the Hamiltonian normalized to the number of lattice sites (referred to as Hamiltonian density throughout) and |φ 0 〉 is a normalized initial random state.The Lanczos method [13] converges quickly to the ground state and low-lying excitations.It consists of iteratively generating a set of orthonormal states {|φ j 〉, j = 0, 1, . . ., M } spanning the Krylov space.Let α j = 〈φ j | ĥ|φ j 〉.Then, α 0 is used to generate an unnormalized orthogonal state: which can then be normalized to obtain the second Lanczos state: Subsequent Lanczos states are generated using the recursion and normalization scheme: with β j = 〈Φ j |Φ j 〉.Notice that acting with 〈Φ 1 | upon Eq. ( 2), yields 〈φ 1 | ĥ|φ 0 〉 = β 1 .Other nonzero matrix elements of the Hamiltonian are obtained by acting with either 〈φ j−1 |, 〈φ j |, or 〈φ j+1 | on the recursion of Eq. ( 4): Thus, the representation of the Hamiltonian in the Lanczos basis is a tridiagonal matrix, which is exact when M coincides with the size of the Hilbert space, D. A low-energy approximation of the Hamiltonian is obtained by truncating the tridiagonal matrix at M ≪ D: In our Lanczos implementation, T M is diagonalized using the method of Multiple Relatively Robust Representations (MR) [81], implemented e.g. in LAPACK [82][83][84].MR was chosen to maximize efficiency because it has O(M 2 ) computational complexity and allows one to specify a range of desired eigenpairs, rather than computing all eigenpairs.This is useful because we are only interested in the lowest eigenvalues of T M , {ϵ j=0,1,...,λ } (with λ ≪ M ), which accurately approximate the low-lying eigenvalues of ĥ.The corresponding eigenstates, {|ψ j 〉} are obtained by transforming to the original basis using the eigenvectors of T M , v j = (v j0 , v j1 , . . ., v j M ): The dominant memory cost of the methods discussed throughout is incurred via the storage of vectors of dimension D (D-vectors).This is because the Hamiltonian is never stored in memory, e.g. as a sparse matrix.Instead, the matrix-vector multiplications encoding the action of the Hamiltonian on a state are carried out "on-the-fly", based on the bit representation of spin states explained in Ref. [85].For completeness, we provide a brief description.Each of the D = 2 N states in a basis of product states of individual spins is encoded by an integer between 0 and D−1, represented by a set of bits.Mapping the lattice sites to the N bit positions and the individual spin states to the value of the bit (either 0 or 1), the Hamiltonian acts on a basis state in one of the two ways.Either the state: i) gets multiplied by a constant that depends on the value of two bits at different positions; or ii) it gets converted to a state encoded by a different integer, obtained by flipping only two bits, and then multiplied by a constant that depends on the values of the two flipped bits.Once a model is translated into these simple ruleswhich are stored at virtually no memory cost -any matrix-vector multiplication boils down to applying the rules to basis states.In particular, the Lanczos recursion requires only two D-vectors (|φ i 〉, |φ i−1 〉) to be stored in memory in each step, i.Consequently, constructing the corresponding eigenstates, |ψ j 〉 entails a second Lanczos recursion in order to regenerate the Lanczos vectors, while accumulating the weighted sum of Eq. ( 7).This can only be done once the eigenvectors v j are obtained at the end of the first recursion.Suppose we are interested in constructing one of the low-lying states, |ψ j 〉.Then, we require an additional vector to be stored in memory so as to accumulate the weighted sum of Eq. ( 7) during the second recursion, implying that the memory cost is dominated by three D-vectors.
Once a low-lying eigenstate, |ψ j 〉, is found, the static expectation value of a quantum observable in that state can be evaluated.If multiple low-lying states are desired, it is still possible to preserve the 3-vector memory cost by carrying out multiple Lanczos recursions to evaluate the relevant expectations for different low-lying states.In contrast, constructing the whole set of eigenstates, {|ψ j 〉} during the second recursion requires as many extra vectors as desired eigenstates to be stored in memory.While the ground state and low-lying excitations are important, there are problems that require knowledge of higher-excited states or even the whole spectrum.Below, we compare different approaches to go beyond low-lying states.

Microcanonical ensemble
In the standard Lanczos algorithm, the core of the spectrum of ĥ is inaccessible.Loss of orthogonality due to finite machine precision impedes convergence beyond low-lying excitations.To complicate matters further, reorthogonalization schemes are computationally expensive [85].An alternative approach is to set a target energy, construct a quasi-eigenstate corresponding to that energy and compute observables using the obtained microcanonical state.

Microcanonical Lanczos
In the microcanonical Lanczos method (MCLM) [17,86], excited states in the core of the spectrum -and thus inaccessible to the "ground state" Lanczos method above -are probed by setting a target energy density, ϵ and finding the lowest-lying eigenpair of v = ( ĥ − ϵ) 2 . ( The lowest eigenvalue found by performing a Lanczos recursion with v approaches 0 and, using its corresponding Lanczos vectors, one can construct the quasi-eigenstate |ψ ϵ 〉.The MCLM converges slower than the standard Lanczos approach in most applications.For example, the microcanonical variant was found to require O(10 3 ) iterations to construct quasieigenstates in the spin-1/2 Heisenberg chain [86].This is to be contrasted to the standard Lanczos algorithm, which typically requires O(10 2 ) iterations to retrive a low-lying state [17].The energy uncertainty reads As a rule of thumb, in Ref. [17], the authors state that in order to resolve the desired energy level with small energy spread σ ϵ /W < 10 −3 , where W is the spectrum width, M ′ ∼ 10 3 iterations are typically needed.Then, the quasi-eigenstate can be used to compute observables reliably.The computational complexity is dominated by two main components.One of them comes from the diagonalizations of the tridiagonal matrices at each Lanczos iteration.Since we use MR [81,83,84] for these diagonalizations, the number of floating point operations for this part scale as The other component comes from the M ′ matrix-vector multiplications in the Lanczos recursion, each carried out "on-the-fly", incurring a cost O(zD log 2 D), where z is the coordination number of the lattice.Thus, the computational effort from matrix-vector multiplications scales as O(zM ′ D log 2 D), where D = 2 N for spin-1/2 systems.
For N ≳ 20, the ratio between the two costs is z D log 2 D/M ′ 2 ≳ 1, and the computational complexity is dominated by the cost of matrix-vector multiplication.However, if the more standard implicit QR method is used for diagonalization instead of MR, the complexity of the diagonalization increases to O(M ′ 4 ).The relevant ratio of computational costs becomes z D log 2 D/M ′ 3 , which only becomes significantly larger than 1 for N ≳ 30.
As a final note on memory cost, we remark that the ĥ2 -term in Eq. ( 8) requires an additional vector to be stored in memory compared with the "ground state" Lanczos, increasing the number of stored D-vectors to four.

Thermal pure quantum states
In this subsection, we follow closely the work of Sugiura and Shimizu [34].The rationale of the TPQ method is to find a pure state that faithfully captures the equilibrium properties of a quantum system at finite temperature as accurately as possible using microcanonical TPQ states with well defined energy, constructed as follows.First, one generates a random state For simplicity, {|i〉} is usually taken as the set of product states of individual spins.The distribution of energy in |φ 0 〉 is proportional to the density of states where s(u; N ) is the entropy density, which converges to a function of the energy density, s(u; ∞), in the thermodynamic limit [34].
The basic procedure is an iterative one similar to minimization annealing schemes.The goal is to modify the distribution of energy in the random state so that it becomes sharply peaked at the desired energy density, ϵ.This is achieved by operating with a suitable polynomial of the Hamiltonian density onto |φ 0 〉 iteratively.Take a constant ϵ upper ∼ O(1), such that ϵ upper ≥ ϵ M , where ϵ M is the maximum eigenvalue of the Hamiltonian density.Thus, ϵ upper is an upper bound on the spectrum of ĥ.Then, start from |φ 0 〉 and iteratively compute, respectively the energy density, u k and the (normalized) new state, |φ k+1 〉 at iteration k: iteratively for k = 1, 2, . . ., N it , the maximum number of iterations.Later on, we will see that N it plays a role analogous to the inverse of the resolution in the Chebyshev expansion.Since the microcanonical states are now generated directly, as opposed to being reconstructed via a Lanczos recursion, the memory cost is now dominated by two D-vectors rather than four.
The first energy density corresponds to the infinite temperature state at β = 0. Thus, g(u; N ) has its maximum at u = u 0 .Then, the energy density decreases gradually towards the ground state energy, ϵ m , as k is increased: One stops iterating at k = N it , when u k gets close enough to the ground state energy density, ϵ m .The obtained TPQ states in the sequence |φ 0 〉, |φ 1 〉, . . ., |φ N it 〉 correspond to decreasing thermal energy densities Thus, an estimate of the equilibrium average value of an arbitrary observable Â is obtained as 〈 Â〉 k = 〈φ k | Â|φ k 〉 as a function of u k .Notably, in the large system limit, the effective temperature associated with each TPQ iteration, β k , accurately reproduces the true thermodynamic temperature of a state with energy density u k .In fact, it is possible to approximate the thermodynamic temperature with an error of order O(1/N 2 ) [34].
Finally, the static expectation value 〈 Â〉 k obtained for each realization of the random coefficients {ξ i } depends exponentially less on the number of sites, N , as the latter is increased, due to self averaging properties.Hence, accurate results are often obtained with a few or even a single random vector realization [34,35].

Chebyshev polynomial Green's function
This method consists of numerically evaluating the lattice resolvent operator Ĝ via an exact expansion in terms of Chebyshev polynomials of the Hamiltonian density.Here, z = ϵ + iη is a complex energy variable.A key aspect is that the Green's function is reconstructed with uniform energy resolution over the entire energy range.For numerical stability, the resolution parameter should satisfy η = Im z ≳ δϵ, where δϵ is the mean level spacing.To expand Eq. ( 14) in Chebyshev polynomials, we consider the following linear transformation of the Hamiltonian and the energy variables: h = ( ĥ − b)/a and z = ε + i η, where ε = (ϵ − b)/a, η = η/a, and where ϵ M and ϵ m are the extremal eigenvalues and f ≃ 1.001 is a safety factor to ensure that the spectrum of the reconstructed operator falls inside the Chebyshev domain of convergence at each iteration step.As customary, we work with Chebyshev polynomials of the first kind {T n (x) = cos(n arccos x) , n = 0, 1, 2 . . .} due to their favorable convergence properties [78].
Typical target functions of energy, including density of states and static expectations values, are evaluated by making use of the Chebyshev polynomial expansion of the imaginary part of the rescaled Green's function [51] with The operators T n ( h) of Eq. ( 16) are constructed using the operator versions of the Cheby-shev recursion relation The series is truncated when the desired accuracy is achieved for a given choice of resolution.The (N poly + 1)-th order approximation of the lattice Green's function is therefore ĜN For most cases, N poly = c η−1 , with c = O(1) sufficing to achieve machine precision [55].
The spectral operator within the CPGF approach can be defined as follows By applying this operator to the r-th realization of the random state |φ (r) 0 〉, we obtain |ε, η〉 r , a quasi-eigenstate with rescaled energy ε (within the rescaled resolution η).To compute static expectation values, we start by defining Provided that the resolution η is adequate (η → δϵ ⇔ η → δϵ/a, where δϵ is the mean level spacing), one obtains an accurate estimate of the expectation value of Â for a given energy ϵ, A(ϵ) by averaging over realizations of the initial random state and using Eq. ( 21):

Canonical ensemble
While microcanonical methods are useful, one might also be interested in evaluating observables using canonical states.In practice, one may decide which ensemble is more convenient to perform a given calculation because the principle of ensemble equivalence guarantees that results are consistent across statistical ensembles.For example, canonical methods have the advantage of allowing direct specification of temperature as an input, so they may be preferable to study temperature dependence of systems in thermal equilibrium.Moreover, for finite systems, calculations done with the microcanonical ensemble tend to show significant statistical fluctuations [17].As temperature increases, higher energy states in the spectrum become increasingly important for determining the properties of the system and these statistical fluctuations are smeared out.The most interesting features of the systems we tackle in this work appear at low temperature, so the canonical methods detailed below are particularly useful in this context.

Finite temperature Lanczos method
The finite temperature Lanczos method (FTLM) has been introduced in Ref. [15] and discussed in depth in Ref. [17].The basic idea is to generate a set of eigenpairs {ϵ j,r , |ψ (r) j 〉} using M FT Lanczos steps and starting from different realizations of the initial random state.Throughout the recursion, we need only store two sets of overlaps: Q r, j ≡ 〈φ 0 〉.The memory cost is still dominated by the D-vectors: two for the recursion and one to store the initial random state so as to allow the computation of the overlaps, totaling three D-vectors.The STE estimator of the canonical average, 〈A〉 (β, N ), defined as in the FTLM [15,17] is obtained as Since the statistical fluctuations of this estimator increase significantly as the temperature is decreased, the need for an optimized low-temperature Lanczos method (LTLM) arose and this method has been introduced in Ref. [16].Apart from the Q r, j overlaps defined above, the computation of the STE estimator for this method requires the storage of O(M 2 FT ) additional overlaps: l 〉, with l, j = 0, 1, . . ., M FT .In terms of these overlaps and using the symmetric form in the right hand side of Eq. ( 23), the final estimator becomes Notice that the LTLM requires a double sum with O(M 2 FT ) terms.It becomes increasingly more expensive to compute these overlaps as the temperature is decreased and more Lanczos iterations M FT are required.However, the estimator of Eq. ( 25) has the advantage of reaching smoothly the zero temperature limit -as opposed to that of Eq. ( 24) -which is the reason behind its advantageous statistical convergence properties.Below, we introduce other alternative methods to FTLM since LTLM has an inherently high cost.

Canonical thermal pure quantum states
The microcanonical TPQ state we outlined previously is specified by the independent variables (u, N ).It can be shown [35] that its (unnormalised) canonical counterpart |β, N 〉 -specified by the inverse temperature, β instead of u -is obtained as follows: A simple analytic transformation reminiscent of the principle of ensemble equivalence allows one to cast canonical TPQ states in terms of their microcanonical counterparts.This correspondence is obtained as follows.First, we assume that the minimum and maximum eigenvalues of ĥ, respectively ϵ m and ϵ M , are known.These can be obtained numerically, for example with Lanczos.Let us now define an unnormalised microcanonical TPQ state for a given realization of the initial random state: where W = ϵ M − ϵ m is the bandwidth.Here, dividing by W ensures some degree of numerical stability since the operator inside parentheses is then bounded.Multiplying and dividing Eq. ( 26) by e N βϵ M /2 and Taylor expanding the exponential, one finds: In the canonical TPQ formulation (CTPQ), the STE estimator is then given by where the subscript r means that the canonical states of Eq. ( 28) have been constructed using the r-th realization of the initial random state.Naively, one might expect that computing this expectation would involve performing a double sum over the iterations and storing O(N 2 it ) overlaps, yielding a cost comparable to LTLM.This is because Eq.( 28) implies that If the observable of interest is a constant of motion (i.e.[ Â, ĥ] = 0), Eq.( 30) simplifies significantly.Let A k,r = 〈k| Â|k〉 r , A ′ k,r = 〈k| Â|k + 1〉 r .Then, we have In such cases, we only need to store 4N it ≪ D overlaps for each random vector: Finally, for each inverse temperature, the STE expectation of Eq.( 29) can be reconstructed using the stored overlaps: While this derivation is only strictly valid when [ Â, ĥ] = 0, in Ref. [35], the authors show that it holds remarkably well in general.We shall confirm this in practice below.
Memory-wise, CTPQ is similar to its microcanonical counterpart, requiring two D-vectors.Unfortunately, numerical instabilities build up rapidly as k increases.Lower-temperature properties are thus challenging to probe.

Finite temperature Chebyshev polynomial approach
So far, we have reviewed the generalization of the ideas behind TPQ and Lanczos to the study of canonical expectations.In what follows, we introduce the finite temperature Chebyshev polynomial (FTCP) method, a new approach that we developed to extend the CPGF method to the canonical ensemble description of interacting quantum systems.
In Ref. [51], the g-coefficients of Eq. ( 16) are obtained by exploiting the operator version of the Jacobi-Anger identity where J n (z) is the Bessel function of order n.We follow a similar route, and seek a Chebyshev expansion of the operator e −β h/2 .Suppose we are interested in low-temperature behavior, i.e. a high inverse temperature, β max .We could expand the operator e −β max h/2 in Chebyshev polynomials directly as done e.g. in Ref. [87].However, such an expansion is vulnerable to numerical instabilities for large β max due to the rapid growth of the Bessel functions.Using where L k=1 δβ k = β max /2, we can decompose e −β max h/2 into a string of L operators.This expansion enables us to bypass the divergent behavior of J n (z) in Eq. ( 33) for large negative imaginary arguments z = −iβ max .Moreover, the inverse temperature steps, δβ k need not be uniform, but may instead vary for each operator, e −δβ k h, in our string of L operators.This opens the door to the use of an adaptive temperature step.
Using the modified Bessel functions -which obey the relation x ∈ -we can use the Jacobi-Anger identity to cast each operator in our string of operators as a Chebyshev series: Applying Eq. ( 35) to a random state |φ (r) 0 〉 produces a sequence of approximate finite temperature states (with inverse temperatures β l = 2 l k=1 δβ k with l = 1, . . ., L), with accuracy controlled by the truncation order.In practice, a Chebyshev truncation order N poly,k ∼ O(10) ensures convergence for a typical inverse temperature step of δβ ≲ 10 2 (in rescaled units).The l-th finite temperature state reads where, as customary, the Chebyshev vectors are generated starting from the r-th realization of the initial random state.We note that, at all steps, the arguments of the fast growing modified Bessel functions need to be kept small in order to avoid numerical instabilities.
As the l-th operator in the string of operators is applied to |φ (r) l−1 〉, the canonical average of a quantum observable, Â, can be evaluated for the l-th inverse temperature.In order to compute a thermal average, it suffices to notice that 〈φ (r) Then, using the RHS of Eq. ( 23), we obtain the STE expectation with the FTCP method: , where 2 The adaptive inverse temperature step used in this work allows us to maximize efficiency by focusing the computational effort at low temperature, where a finer temperature grid (and thus a larger spacing δβ k ) is required to capture the key features of the systems at play.The reconstruction of 〈A〉 STE (β l , N ) for a discrete set of L inverse temperatures, {β l , l = 1, 2, . . ., L} involves storing 2L overlaps, 〈φ l 〉 for each random vector realization.The total number of Chebyshev iterations is thus N Cheb = L l=1 N poly, l ∼ 10 3 .As shown shortly in Sec. 4, the FTCP favorable convergence properties will allow us to reach very low temperatures that are hard to access with FTLM and CTPQ.Finally, FTCP has the same memory requirement of three D-vectors as CPGF: two for the Chebyshev recursion and one to cumulatively generate the finite temperature state at each inverse temperature.

Dynamical properties
The prototype simulation aimed at studying the dynamics of a quantum system starts from a well defined initial state |Ψ(t = 0)〉, such as the ground state of the model Hamiltonian at play, |GS〉, which can be obtained, e.g. using the Lanczos method.This initial state is then evolved using the time evolution operator.Successive small time steps are taken in order to maintain enough numerical accuracy, while keeping track of the evolution of quantities of interest, such as time-domain correlators of the type where Â, B are two generic quantum observables in the Heisenberg picture.Both Lanczos and a Chebyshev-based approaches exist to approximate the time evolution operator.Within the Lanczos approach, the time evolution operator for a short time step δt is approximated as where {ϵ j } and {|ψ j 〉} are sets of energies and corresponding eigenstates obtained using M t Lanczos steps and starting the Lanczos procedure from a previously computed state |Ψ(t ′ )〉.
On the other hand, the Chebyshev approximation of the time evolution operator -which is used e.g. in Ref. [88] in combination with CTPQ 4 -relies on Eq.( 33): Both methods require O(10) iterations for a standard time step δt ≈ W −1 [17].Yet, the Chebyshev approach has an important advantage.Unlike Lanczos, where a tridiagonal matrix has to be diagonalized at each time step to generate the coefficients of the Lanczos expansion, the coefficients in the Chebyshev expansion in Eq. ( 41) can be easily and efficiently evaluated using freely available numerical libraries.
The efficiency of the Chebyshev approximation of the time evolution operator suggests that a Chebyshev approach can also be advantageous when studying zero-temperature spectral functions, C B Â(ω), obtained by Fourier transforming the time-domain correlators of Eq. ( 39): where we recall that ϵ m is the ground state energy density, obtained e.g. with Lanczos.

Dynamical autocorrelation response functions with Lanczos
In the particular case where B = Â † , the spectral function C Â † Â is referred to as the autocorrelation response function for observable Â and is defined as follows (with z = ω + ϵ m + iη): Once the ground state, |GS〉, is reconstructed with Lanczos, the response function above can be computed by performing an additional Lanczos recursion (where the number of iterations needed for satisfactory convergence is typically M ∼ 10 3 ) with the initial state Similarly to Sec. 3.1, this recursion also generates a (much larger) tridiagonal matrix whose entries can be used to compute the response function [17] with no need to compute the eigenpairs {ϵ j , ṽj , j = 0, 1, . . ., M }.The resolvent (z − ĥ) −1 can be approximated using a continued fraction, thus giving the "Lanczos" response function where the continued fraction is terminated with β M +1 = 0.This procedure is signficantly more expensive computationally than simply approximating the ground state with Lanczos as described in Sec.3.1.This is due to the accumulated cost of matrix-vector multiplications as more iterations are completed, which is O(z M D log 2 D).Unlike Sec.3.2.1, the accumulated cost of diagonalizing the tridiagonal matrix at each iteration using MR [81,83,84] only applies to the first recursion.This diagonalization cost is O(M 3 ) ∼ 10 6 , which is very small compared to the cost of matrix-vector multiplications, e.g.O(z M D log 2 D) ∼ 10 12 for N = 24.

Hybrid Lanczos-Chebyshev method for spectral functions
Here, we use the Chebyshev expansion of the resolvent operator of Eq. ( 19) to compute spectral functions directly in the frequency domain.This approach is inspired by a similar technique described in Ref. [49], where a kernel polynomial approximation based on Chebyshev polynomials is used.This approach was further exploited in Refs.[50,74], where Chebyshev expansions were combined with Matrix Product States (MPS) and DMRG to investigate onedimensional strongly correlated systems.Yet, this technique has so far relied on the use of a kernel convolutions to damp Gibbs oscillations in the Chebyshev expansion.Here, we combine the numerically exact Chebyshev expansion of Eq. ( 19), which avoids the use of a kernel, with Lanczos.The key advantage of this approach is the rigorous control over resolution, a feature that is shared with the CPGF method that was described above in Sec.3.2.3.In principle, the ideas of the method described below could be combined with MPS and DMRG as well, but such a task is outside the scope of this work.
Similarly to Lanczos, we start the procedure with the state | φ0 〉 -obtained from two prior Lanczos recursions -and, instead of a third Lanczos recursion, we carry out a Chebyshev recursion to generate the polynomials T n ( h) using Eq. ( 18), while storing the moments The autocorrelation response function is then obtained as follows: with z = (ω + ϵ m + iη − b)/a, where a, b are defined in Eq. ( 15) and we recall that ϵ m is the ground state energy density.This procedure has significant advantages.The first is that two moments can be obtained per matrix-vector multiplication, which implies that the the numerical effort is halved (other parameters being fixed).Thus, the computational effort is proportional to the number of iterations, Ñit = Ñpoly /2, where Ñpoly is the number of Chebyshev moments needed for convergence.This is derived by applying the product identity 2T m (x)T n (x) = T m+n (x) + T m−n (x) to the overlaps 〈 φn | φn 〉 and 〈 φn+1 | φn 〉, where For a given iteration, given a new | φn 〉, two moments can now be computed: Two remarks are now in order: • Ñit is typically of the same order of magnitude as M , which guarantees that CPGF is at least as fast as Lanczos.In practice, we observe faster performance with CPGF.We attribute this to the possibility of better parallelization with CPGF because it requires half as many vector update loops.These loops are needed in order to carry out the Lanczos and CPGF recursions with only two vectors of dimension D stored in memory.They incur a cost that, whilst not dominating over that of matrix-vector multiplications, still compares closely.In fact, the complexity of each of these loop is proportional to D. With Lanczos, the two steps of Eq. ( 4) involve two of these loops that need to be executed one after the other, in opposition to CPGF, which needs only a single loop vector update.
• CPGF can easily be modified to compute more general spectral functions (for which we might have B ̸ = Â † ) without a significant additional memory or computer time cost.The 3-vector memory cost is preserved because the vector used to generate |GS〉 during the second Lanczos recursion is not used in CPGF once the initial state, | φ0 〉 is generated.Thus, this vector can be used to store |ϕ〉 ≡ B † |GS〉/ 〈GS| B B † |GS〉, which in turn can be used to compute the modified moments, In contrast, the continued fraction Lanczos approach does not work in the case B ̸ = Â † .One must then resort to Eq. ( 40) to directly study the behavior of the time-domain correlator.This The bond-directional character of the Kitaev interaction implies that to each bond corresponds a different type of interaction.Similarly to the case of the Ising model on the triangular lattice, where we have geometrical frustration, here we have exchange frustration due to the nature of the interaction and it is not possible to find a spin configuration that simultaneously minimizes the energy on all bonds.
leads to short time expansions with M t Lanczos vectors and the initial state The eigenvectors of the tridiagonal matrix of Eq. ( 45), {ṽ j } give 〈 ψj | φ0 〉 = ṽj0 .However, the overlaps of the type 〈GS| B| ψj 〉 must be evaluated explicitly using the vector B|GS〉, which now has to be stored in memory separately, thus adding to the memory cost: Moreover, we must update the initial state of the Lanczos expansion at each time interval, |Ψ(t)〉 using short time Lanczos expansions.Then, we re-compute the eigenvectors of a new tridiagonal matrix and re-evaluate 〈GS| B| ψj 〉 for each time step.This process becomes computationally expensive very quickly since we may require a large number of time steps to capture important features of G B Â(t ).On the other hand, the CPGF treats the cases B ̸ = Â † and B = Â † on equal footing.Therefore, the CPGF is a general purpose approach, which accesses spectral functions for the case B ̸ = Â † using the same methodology and with the same computational complexity and memory requirements as the case B = Â † .

Applications
In this section we apply the methods described in Sec. 3 to two generalized Kitaev models on the honeycomb lattice for a 24-spin hexagonal cluster with periodic boundary conditions (see left panel of Fig. 2): the K-H and the K-I models.

Kitaev-Heisenberg model
The K-H model combines Kitaev and Heisenberg interactions.The Hamiltonian can be cast as a sum over NN bonds 〈i, j〉 γ on the honeycomb lattice (the superscript refers to the type of bond, see right panel of Fig. 2): with K = sin ϕ, J = cos ϕ, and where γ = x, y, z is one of the three types of bond on the honeycomb lattice.We use the conventions of Ref. [5]: ϕ ∈ [0, 2π] parameterizes the strength of each term and ensures that all possible ratios of the Heisenberg and Kitaev interactions are considered, and an overall energy scale is defined and set to unity throughout: The Kitaev interaction is bond-directional, i.e. for each distinct type of bond γ = x, y, z, there is a correspondent interaction (respectively S x S x , S y S y , S z S z ), as shown on the right panel of Fig. 2. The calculations presented throughout the next subsections serve two purposes: to prove that the classes of methods considered in this work are consistent and correctly capture the physics of frustrated quantum magnets and, more importantly, to show that the Chebyshev polynomial-based approaches offer significant advantages in terms of performance, stability and generality, especially when combined with Lanczos algorithms.

Zero-temperature consistency and microcanonical approaches
Firstly, we reproduced the results of Refs.[3][4][5] using ED (via the "ground state" Lanczos algorithm), as used in the original references.Then, we recovered these results using both microcanonical and canonical approaches based on Lanczos, TPQ and Chebyshev recursions.
In order to bench-mark the microcanonical approaches -MCLM, MTPQ and CPGF -we set the ground state energy obtained from Lanczos as the target energy.Crucially, MTPQ and CPGF (and their canonical counterparts) require an estimate of the maximum eigenvalue as well.Although we could have applied Lanczos to the operator − ĥ to obtain it, the Hamiltonian of Eq.( 54) happens to have a symmetry that can be used to avoid this additional computation: − ĥ(ϕ) = ĥ((ϕ + π) mod (2π)).Thus, the maximum eigenvalues can be obtained by simply reorganizing the minimum eigenvalues as a function of ϕ and switching their sign (see Fig. 3).These are the minimum and maximum energies that are then used as inputs to the TPQ and Chebyshev methods.We remark that even though the ground state energy can be accurately estimated solely using MTPQ, the method requires an upper bound on the maximum eigenvalue (that one would normally obtain from Lanczos anyway).Moreover, MTPQ has much slower convergence to the ground state than Lanczos, as we will also show later, so it is simply more efficient to use Lanczos to obtain extremal eigenvalues.Still, MTPQ provides us with a useful consistency check, while being less memory-intensive -requiring 2 rather than 3 vectors of size D -and giving access to good approximations of finite temperature states.
The results we present throughout this section are for the ground state energy density and NN spin-spin correlation of the K-H model.The NN spin-spin correlation is used for our benchmark for two reasons.Firstly, in Ref. [5], the authors show that longer-range correlations vanish in the vicinity of the spin liquid phases.Given that we are particularly interested in this region of the phase diagram, it is reasonable to focus on NN correlations.Secondly, the step-like behavior of the NN spin-spin correlation coincides with quantum critical points [5].Moreover, the behavior of the NN correlation as a function of the temperature is intimately connected to peaks in the specific heat [10] -that we also compute using Lanczos, TPQ and Chebyshev -and that are particularly relevant experimentally [89].
For each value of ϕ, we computed the minimum eigenvalue, or ground-state energy of the Hamiltonian, ϵ GS , using the Lanczos ED technique.In the top panel of Fig. 4, we show these results, along with the second derivative −d 2 ϵ GS /dϕ 2 , which accurately detects quantum phase transitions between magnetically ordered phases (Néel, stripy and zigzag antiferromagnets and a ferromagnet) and two spin liquid phases, which we refer to as "antiferromagnetic" (ϕ ∼ π/2) and "ferromagnetic" (ϕ ∼ 3π/2) QSL phases, or AFQSL and FQSL in short.The phase boundaries we obtained are summarized in Table 1 and agree well with Ref. [5].The Lanczos algorithm shows small statistical fluctuations, enabling one to use a single realization of the initial random state across the entire phase diagram.As mentioned in Sec.3.2.2, a single random vector suffices to achieve good accuracy in the MTPQ approach as well due to the self-averaging properties of the TPQ estimator.In fact, error bars obtained with more realizations are negligibly small and thus are not shown in our plots.We now move gears to the spin correlator.Figure 4 shows excellent agreement between the spin-spin correlation computed with Lanczos and MTPQ.As mentioned earlier, the step discontinuities in the correlator signal the quantum phase transitions of the K-H model [3][4][5].MTPQ achieves its maximum effective resolution at low-temperatures (i.e.large N it.) where Fig. 4 shows that it accurately approximates the ground state [34][35][36]90].These methods are designed to reconstruct the ground state in a recursive fashion.Thus, it is perhaps not too surprising that the energy and N-N spin correlation can be both computed with great accuracy provided enough iterations are completed.Also shown in the bottom panel of Fig. 4 is the behavior of the spin correlation for the target energies larger than the ground state by 10% of the spectrum width, obtained with CPGF (the MCLM results coincide, and thus are omitted to avoid unnecessary clutter on the figure).These excited states are still close enough to the ground state that the spin correlation preserves its general shape, albeit with a considerable broadening of its sharper features.The correlator is also computed with CPGF for states with target energy ϵ target .The symbol 〈i, j〉 denotes an average over all nearest neighbors 〈i, j〉.
The top panel in Fig. 5 summarizes a careful convergence study aimed at understanding how the optimal number of iterations, N * it , 5 depends on the microscopic details of the model.Our simulations show that convergence is relatively fast in the purely Heisenberg limits (ϕ = 0, π) and becomes slower as we approach the quantum phase transitions and, in particular the Kitaev limits (ϕ = π/2, 3π/2), in both Lanczos and MTPQ.Note that in spite of the strong dependence of N * it with ϕ, Lanczos converges considerably faster throughout the K-H parameter space, requiring about 2 orders of magnitude less iterations for convergence than MTPQ.
Next, we address convergence around the quantum critical points near the Kitaev limits in more detail; see Fig. 5 (bottom panel).The convergence of the Lanczos and TPQ estimators for the spin-spin correlation is found to slow down notoriously as the Kitaev term on the Hamiltonian becomes dominant (|K|/|J| ≫ 1) and the ground state energy per site ϵ GS (ϕ) and the NN spin-spin correlation 〈S i S j 〉 〈i, j〉 GS (ϕ) start to differ significantly.As a result of the slower convergence away from the purely Heisenberg points, the computational effort grows substantially, notably so close to the Kitaev points.For example, for ϕ ≃ 0.506π ⇐⇒ |K|/|J| ≃ 53, MTPQ requires around 4 × 10 4 iterations for optimal convergence as defined earlier to be achieved.
Figure 6 shows the energy density and spin correlation function of the 24-site cluster across the phase diagram calculated with CPGF.We compare the latter with MCLM, which also has the ability to probe an input target energy.We also establish the accuracy of the STE of Eq. ( 22) using CPGF.A detailed analysis shown in Appendix A confirms that the standard deviation of the estimate for the N-N spin-spin correlation function scales as expected, i.e. as the inverse square root of the number of realizations of the initial random state.
These results also show that finer CPGF resolutions are needed to probe the phase diagram around the Kitaev points (ϕ = π/2, 3π/2).On the contrary, near the Heisenberg limits (ϕ = 0, π), convergence occurs with comparatively coarse resolution.A useful feature of the CPGF method is that the optimal number of Chebyshev iterations, N * poly = N * poly (η), follows a predictable pattern: it is roughly proportional to the spectrum width and inversely proportional to the required resolution [51,55].The results reported in Appendix B confirm this expected behavior.We find that the convergence speed of MCLM and CPGF compares very differently throughout the phase diagram.MCLM seems to show unpredictable behavior as some points of the phase diagram require significantly more computional effort than others.CPGF behaves more intuitively, requiring a greater effort for points close to the phase transi- tions, which need finer resolutions and, consequently more iterations and thus more computer time.Although the CPGF approach is significantly faster in some regions of the K-H parameter space, the MCLM approach is faster close to the transitions (see Appendix B).Nevertheless, when targeting the ground state, CPGF is 25% faster on average.For the excited states, the difference in performance is much more pronounced and the results of Appendix B show that CPGF is the method of choice due to its faster overall convergence (about an order of magnitude less CPU time required).
Next, we bench-mark the MCLM and CPGF methods around the Kitaev limits in detail, i.e. for spreads of ϕ-values in the intervals [0.4800, 0.5200]π and [1.5160, 1.5760]π.To assess convergence in the K-H model, we carefully track the evolution of the correlation function as more iterations are completed with each of the two methods.
In Fig. 7, we show the dependence of the NN spin correlation on the number of iterations in the MCLM and CPGF approaches.Here, we focused on a range of resolutions in the interval [2.5 × 10 −5 , 5 × 10 −4 ], but only plotted the curves corresponding to the energy resolutions that yield convergence.In this case, we consider that a given resolution yields convergence when the energy density and spin correlation no longer change appreciably as finer resolutions are considered.To ensure that convergence has indeed occurred, we compute the difference between the energy density for two consecutive resolutions and ensure that the difference is smaller than the resolution itself.As it turns out, convergence is achieved for N poly ∼ 10 2 −10 3 depending on the value of ϕ. Figure 7 confirms that enhanced resolutions are crucial in the vicinity of quantum critical points.Moreover, it also shows that MCLM and CPGF have complementary convergence behavior.The top left panel shows that MCLM converges quicker around  0.500π, but dramatically slows down around 0.480π and 0.520π.The top right panel shows that CPGF has more consistent and faster convergence, with only the points closest to the transition (0.494π, 0.506π) requiring more iterations due to the need for a finer resolution.The bottom panels exhibit a similar tendency.The points near 1.576π display faster convergence with CPGF than with MCLM, while the curves in the center of the two plots show that convergence is faster with MCLM due to the increased resolution that is needed with CPGF.Here, we emphasize that we found each iteration to be faster with CPGF due to a lower number of necessary operations than with MCLM.Thus, CPGF is still faster even in situations where the number of iterations required for convergence is similar, e.g. for the points near 1.516π in the bottom panels.The main contribution to the extra time per iteration in MCLM is the action of ĥ2 upon a state at each iteration.Since the Hamiltonian is generated "on-the-fly", one must act with ĥ twice in order to obtain the action of ĥ2 .Thus, each iteration incurs an extra cost due to the additional matrix-vector multiplication, which is the most computationally expensive operation in these methods.

Canonical approaches: Finite temperature
In what follows, we discuss the performance of the finite-temperature approaches introduced in Sec.3.3.We also include the microcanonical variant of the TPQ method in this analysis and confirm that the temperature corresponding to each energy density can be accurately estimated as explained in Ref. [34].For concreteness, we focus on the Kitaev point at ϕ = 1.5π and restrict the temperature to the range T = [0.004,24] in units of the Kitaev coupling, which suffices to capture the relevant features of the model.We set out to obtain the NN spin correlation, specific heat and entropy with all methods, with a careful convergence analysis.Thus, we compute the minimum (i.e., optimal) number of iterations such that the target functions are reliably captured within the desired accuracy at all temperatures.Naturally, the relevant convergence parameters need to be adjusted separately for each method due to their different characteristics.These are summarized in Table 2.It is important to note that the results presented below not only agree with each other, but also with QMC and exponential tensor renormalization group (XTRG) studies of Ref. [10], thus supporting the validity of our implementation for all methods.The specific heat is computed as follows: c = N β 2 (〈 ĥ2 〉 − 〈 ĥ〉 2 ).Similarly to Ref. [91], the entropy density is computed by integrating c/T .We perform the integral numerically using Simpson's rule.At first sight, Table 2 seems to indicate that FTLM is the most efficient method.However, notice that in Fig. 8, FTLM shows large low-temperature fluctuations, reaching about 26 times the fluctuations of FTCP (calculated from root-mean-square deviations).The inset of Fig. 8 shows that these fluctuations are larger for FTLM than for the other methods throughout the temperature range.Here, we note that the same initial random vectors are used in all methods so that statistical fluctuations can be directly compared.Moreover, the total computational cost is dictated not only by the number of iterations, but also by the number of matrix-vector multiplications per iteration.This number is higher for FTLM than FTCP, leading to about twice the average computer time per iteration (2.05 seconds with FTLM compared to 0.98 seconds with FTCP).
Figure 9 shows a very important shortcoming of the CTPQ method.Unlike the NN spin correlation, the specific heat has an important feature that cannot be reproduced with CTPQ, Figure 8: Finite temperature N-N spin correlator computed with MTPQ, FTCP, CTPQ and FTLM for the 24-site cluster at the Kitaev limit ϕ = 3π/2.We used 50 random vector realizations in all cases.Error bars are negligibly small, except for FTLM.Inset: Standard deviation of the N-N spin correlator.CTPQ is not shown because the fluctuations are comparable to MTPQ.The error bars are negligibly small in most of the temperature range, except for temperatures below that of the low-temperature peak.The error bars are of comparable size for every method (except for CTPQ, for which their meaning is ill-defined due to numerical instability).Model parameters as in Fig. 8.
namely the decay to zero of the low-temperature peak as T → 0. The limitation is related to the need to consider more terms in the summation of Eq. ( 31) so as to achieve convergence.We find that when we considered 500 terms, we reached the limit of machine precision (80bit floating-point on an Intel Core i5 processor) before reaching enough accuracy.The CTPQ results gradually lose accuracy as temperature decreases and at some point between 10 −1 and 10 −2 , they are no longer reliable.The entropy also shows signs of the numerical instability of CTPQ (see inset of Fig. 9).FTCP avoids this instability via the numerical stabilization procedure described in Sec.3.3.3.The operator break-up into the product of Eq. ( 34) ensures that each Chebyshev expansion in Eq. ( 36) is stable.The arguments of the fast growing modified Bessel functions are guaranteed to remain in check because they can be controlled via the inverse temperature step.This is to be contrasted with Eq. (31), where the also fast growing functions of the inverse temperature and the number of iterations are not controlled, leading to the numerical instability visible in Fig. 9.
Having discussed the limitations of FTLM (large low-temperature fluctuations) and CTPQ (numerical instability), we now move on to MTPQ.An often ommited detail in the literature is that the upper bound on the maximum eigenvalue of the Hamiltonian density, e M , given as an input to MTPQ determines: i) the temperature range that is covered, and ii) how much accuracy is achieved for a given temperature.The strong influence of e M can be traced back to the evolution of the TPQ energy density distribution with the number of iterations [34].We found that in order to cover the desired temperature range with enough accuracy in MTPQ, we had to increase the upper bound to 35 times the maximum eigenvalue of the Hamiltonian density (computed with Lanczos).We consider that enough accuracy has been reached when the specific heat and entropy computed with MTPQ match the FTCP and FTLM results.
The results of Fig. 9 indicate that the disparity between MTPQ with e upper = e M and the other methods is solely due to insufficient accuracy at high temperature.While MTPQ captures the high-temperature peak of the specific heat even with e upper = e M , the amplitude of this peak is severely underestimated.As more iterations are completed, the accuracy of the method increases and the low-temperature peak matches the results of other methods better, but still not perfectly.Moreover, even though the MTPQ result with e upper = e M captures the general behavior of the entropy density, there is not enough accuracy because of accumulated error from the lack of accuracy at high temperature and the reduced temperature range (note that we compute the entropy as an integral over temperature).The other methods match perfectly and show comparable fluctuations.Since there are only small error bars of approximately the same size for all methods at low temperature, the disparity can only be due to insufficient accuracy for MTPQ with e upper = e M .By probing higher values of e upper , we find that is necessary to input e upper ≈ 35e M for the high and low-temperature peaks to reproduce the correct behavior.When we consider this upper bound, MTPQ requires about twice the total computer time compared to FTCP in order to achieve comparable results (14268 versus 7376 seconds per core).Even though the number of iterations is nearly 3 times higher in MTPQ compared to FTCP, the former requires 40 % less matrix-vector operations per iteration.Still, FTCP remains about two times faster due to significantly faster convergence.
Finally, notice that for both the specific heat and the entropy in Fig. 9, FTLM does not show the low-temperature fluctuations of Fig. 8.We attribute this to the fact that Lanczos methods are designed to quickly achieve an approximation of the Hamiltonian in subspace restricted to the ground state and low-lying excitations.The specific heat and the entropy are calculated solely in terms of averages of the Hamiltonian density, ĥ and ĥ2 .Thus, convergence is better for these quantities than for more general observables, such as the NN spin correlation, which are not as closely associated with the Hamiltonian.

Kitaev-Ising model
The K-I model combines Kitaev and Ising interactions.The model Hamiltonian is: where we use the parametrization of Ref. [61].We allow not only for the isotropic case, but also a particular type of anisotropy, for which we have: where α ∈ [0, 1.5] is a parameter and the energy scale is set by In Ref. [61], the authors study this model for the 24-spin hexagonal cluster of the left panel of Fig. 2 using Lanczos at T = 0. We start by reproducing some of their results so as to validate our implementation.Then, we present a new study displaying dynamical signatures of the phase transitions described in Ref. [61].These signatures are found in the σ z spin susceptiblity, which we obtain using both the Lanczos and the hybrid Lanczos-Chebyshev (HLC) approach that we introduced in Sec.3.4.2.

Lanczos bench-mark
Figure 10 summarizes our Lanczos results and reproduces those of Ref. [61].The top panel shows two steps in the first derivative of the ground state energy density for fixed α = 0.7 and varying J I ∈ [10 −4 , 10].Correspondingly, the center panel displays two peaks in the second derivative.These steps/peaks indicate two quantum phase transitions.From left to right, the first quantum critical point separates a quantum spin liquid (blue) from a nematic (green) phase and the second one separates the latter from a ferromagnetic phase (red).Similarly to how the steps in the spin-spin correlation signal transitions for the K-H model (see bottom panel of Fig. 4), the mean squared magnetization also signals transitions in the K-I model, with its value approaching saturation very quickly as we enter the ferromagnetic phase.

Finite temperature: Comparing thermal pure quantum states and Chebyshev
In section 4.1.2,we bench-marked the FTCP approach by computing the temperature dependence of the specific heat, entropy density and NN spin correlation.In principle, MTPQ could These results match those of Ref. [61] and illustrate the transitions between ferromagnetic (red), nematic (green) and quantum spin liquid (blue) phases found in Ref. [61].Here, a single realization of the initial random state was found to be sufficient due to negligible statistical fluctuations.be the most viable competitor of FTCP (despite being a microcanonical method) because it avoids the shortcomings of FTLM and CTPQ.Yet, we found that MTPQ required about twice the computer time of FTCP in the context of the K-H model.Here, we further bench-mark FTCP by considering the K-I model.We also seek to clarify whether FTCP outperforms MTPQ in terms of computer time for a different model, which would suggest that the advantages of FTCP are not problem-dependent.We start by recovering the MTPQ results of Ref. [61] using our own implementation.Then, we repeat the calculation using FTCP.These results -shown in Fig. 11 -are for two specific points of the phase diagram: (α, J I ) = {(0.7,0.001), (0.7, 0.03)}.Going back to Fig. 10, we can see that these two points are located within the Kitaev quantum spin liquid and nematic regions of the phase diagram, respectively.
In Ref. [61], the authors remark that the well known results for the pure Kitaev model that we recovered in Figs. 8 and 9 remain qualitatively valid even when J I ̸ = 0, and even within the nematic phase (see right panels of Fig. 11 with results for (α, J I ) = (0.7, 0.03)).The specific heat has a two-peak structure, corresponding to a two-step release of entropy.At high temperature (T ∼ 0.5), the Majorana fermions c (defined in [61]) release their entropy (0.5k B ln 2).The other half of the entropy is released by Z 2 fluxes at low temperature (T ≲ 10 −2 ) [8].The high temperature crossover coincides with an enhancement of the expectation of the kinetic energy of the Majorana fermions c, defined in terms of the spin operators as follows: The features mentioned above are all apparent in Fig. 11 and our MTPQ and FTCP results match nearly perfectly.The striking difference between the two approaches is that once again, FTCP cuts the computer time in approximately half.To be more precise, for J I = 0.03, FTCP is around 2.1 times faster and for J I = 0.001, it is around 1.8 times faster.Table 3 summarizes these differences in computer time.
Figures 11 a),b) confirm the two-peak behavior of the specific heat.In the right panel (J I = 0.03), statistical fluctuations are more apparent and even doubling the number of random vectors compared with the case J I = 0.001 (going from 50 to 100), we find that statistical fluctuations remain higher.This is illustrated in the low temperature behavior of the standard deviation of the specific heat estimator, which is shown on the insets.MTPQ and FTCP show identical statistical properties, with these standard deviations matching remarkably well.Notice that the optimal value for the upper bounds on the maximum eigenvalue of Hamiltonian used in MTPQ were different in each case (40e M for J I = 0.001 and 30e M for J I = 0.03).These were chosen so as to ensure enough accuracy throughout the chosen temperature ranges (T ∈ [0.002, 24] for J I = 0.001 and T ∈ [0.004, 24] for J I = 0.03, in units of K x + K y + K z ).Figures 11 c),d) show the two-step release of entropy.Compared with the pure Kitaev case of Fig. 9, the left panel (c) shows a much more pronounced plateau-like behavior between T = 10 −1 and T = 10 −2 , ending in an abrupt decrease of entropy.In contrast, the right panel (d) shows no plateau at all, with a gentler decrease in entropy between T = 10 −1 and T = 10 −2 .This is a manifestation of the intrinsic differences between the two liquid phases (Kitaev QSL and nematic).Finally, figures 11 e),f) illustrate the high temperature enhancement of the kinetic energy of the Majorana fermions c, a behavior that is shared between the two phases.Here, statistical fluctuations are very small for both MTPQ and FTCP, with negligible error bars.

Dynamics: Hybrid Lanczos-Chebyshev approach
Lastly, we present novel results that elaborate on the picture of the K-I system that was outlined in Ref. [61].We find that the signatures of the quantum phase transitions are present not only in static quantities, such as the energy and squared magnetization, but also in the dynamical spin susceptibility.This spectral function is obtained by considering the relevant observable in Eq. ( 48) to be the Fourier-transformed spin operator, i.e.Â = r e −iq•r Ŝz r / N , where r is a position on the lattice and q is the wave vector.
In Fig. 12, we show the variation of the q = (0, 0) dynamical spin susceptibility, S Γ (ω), with the model parameter J I .These results are obtained with the hybrid Lanczos-Chebyshev Figure 12: Dynamical spin susceptibility of the K-I system for varying J I and α = 0.7, normalized to its maximum value.The phase transitions are marked as gray dashed lines.The results obtained with the Lanczos and HLC methods are identical.The white space corresponds to a vanishing S Γ (ω), as shown on the bottom of the color bar.
(HLC) method.The initial Lanczos run is stopped when the variation between the ground state energy density computed for consecutive iterations is less than 10 −9 in units of K x + K y + K z .We compute 3000 Chebyshev moments, which is enough to achieve convergence for all the values of η considered.Specifically, for each J I , the resolution parameter is fixed to 0.1% of the spectrum width, which translates to values in the interval η ∈ [0.0011, 0.0036](K x + K y + K z ).In terms of statistical sampling, we find that lower values of J I require more random vectors for the peaks of the dynamic susceptibility to be resolved satisfactorily.Thus, for J I < 0.0057, we use 16, rather than the 4 initial random vectors that we use for J I > 0.0057.Our interpretation of these results is in line with a similar reasoning for the K-H model presented in Ref. [5], albeit with a crucial difference due to the specifics of the liquid-to-liquid transition.In the ferromagnetic limit (J I ≳ 10 −1 ), the ω = 0 component dominates because of the strong ferromagnetic correlations.As J I is lowered and the transition to the nematic phase occurs, the gapless magnon mode gradually turns into a gapped mode.Moreover, there is a proliferation of sharp well-defined excitations in the nematic phase, which abruptly collapse onto a smaller set of modes as the transition to the Kitaev phase occurs (J I ∼ 10 −2 ).This rapid change in the spin susceptibility is consistent with the T = 0 first order topological phase transition described in Ref. [61].The gap is found to peak for the Kitaev liquid, at which point the lowest-ω mode occurs for a larger ω than in the nematic phase.
We close this section with a note on computational resources and memory cost.Throughout the paper, our calculations are done using Intel Xeon Gold 6138 processors running at 2 GHz and each simulation requires between 0.5 and 1 GB of memory.All Chebyshev methods showed improved performance.In the case of the dynamics studies of this section, we found that under the exact same circumstances, our HLC method is 33 % faster than the continued fraction Lanczos approach.This estimate was obtained as follows.We reproduced a previously known result for the dynamical spin susceptibility from Ref. [5] using the same number of Lanczos iterations/polynomials and averaging over ∼ 400 simulations, each running with parallelization enabled with 16 cores using the aforementioned processors.We used our own implementation of Lanczos in both cases, so our results are implementation-independent.
The main conclusion of this section is that our newly introduced Chebyshev approach is not only remarkably flexible, in the sense that it allows the study of generic (not necessarily autocorrelation) spectral functions, but also faster overall than the traditional Lanczos approach.

Concluding remarks
We studied the Kitaev-Heisenberg (K-H) and Kitaev-Ising (K-I) models on the honeycomb lattice for a 24-spin hexagonal cluster with periodic boundary conditions using three distinct approaches: Lanczos, TPQ and Chebyshev.This work is mainly divided in three parts.In the first part, we started by reproducing the results of Refs.[3][4][5] for the K-H model using Lanczos ED.Then, we recovered those results using microcanonical variants of the Lanczos (MCLM) and thermal pure quantum state (TPQ) methods and the Chebyshev polynomial Green's function (CPGF) method, independently of each other.For these three methods, we carefully examined the spectral and statistical convergence properties.While Lanczos is found to be ideal to approximate the ground state, we find that CPGF is the most efficient method capable of probing an arbitrary target energy with well controlled accuracy, proving to be faster than both MTPQ and MCLM on average throughout the phase diagram of the K-H model.
In the second part, we computed the temperature dependence of the N-N spin correlation, specific heat and entropy density.The aim of this study was to compare three methods based on Lanczos, TPQ and Chebyshev ideas, respectively: the finite temperature Lanczos method (FTLM), the canonical variant of TPQ (CTPQ) and the finite temperature Chebyshev Polyno-mial method (FTCP), introduced in this paper.The MTPQ method is also considered because it is capable of estimating the temperature corresponding to each energy density remarkably accurately.Our implementations are bench-marked against the exponential tensor renormalization group results of Ref. [10].We find our newly introduced FTCP method to be the most efficient and versatile of the three, namely showing a two-fold increase in speed compared with TPQ, while also avoiding the large low-temperature statistical fluctuations of FTLM.
The third and last part of this work started with the reproduction of Lanczos ED and TPQ results for the K-I model [61].This bench-mark allowed us to validate our implementation and carry out a novel dynamical study for the K-I model, where we used a hybrid Lanczos-Chebyshev method.The latter was also shown to be more flexible and about 33% faster than Lanczos on average.Our detailed calculation of the dynamical spin susceptibility identifies signatures of the quantum phase transitions in the K-I model.This third part of our work also provides further evidence for the efficiency of FTCP, confirming the two-fold speed-up with respect to TPQ in finite temperature calculations for the K-I model.
In what follows, we summarize the key aspects that support the conclusions above for each part of our work.
All microcanonical methods show low statistical fluctuations and considering even just a single realization of the initial random state seems to suffice to achieve negligible deviations from the exact diagonalization results throughout the phase diagram of the K-H model.Unlike the number of realizations, the resolution plays a central role when comparing the performance of the three microcanonical approaches.In CPGF, finer resolutions always require more polynomials to achieve convergence.There is no one-to-one correspondence between the CPGF resolution and the effective resolutions of MTPQ and MCLM (which vary as more iterations are completed).Nonetheless, we managed to compare the three methods.TPQ showed an erratic convergence speed, with particularly significant slow-down for ϕ = 0.506π in the K-H model, at which point around 4 × 10 4 polynomials/iterations are needed to achieve convergence and thus match the Lanczos ED results satisfactorily.On the other hand, CPGF requires very fine resolution to recover the ED results near the quantum phase transitions, thus converging relatively slower (but still faster than MTPQ) at these points of the phase diagram.Yet, away from quantum critical points, comparatively coarse resolutions are enough to reproduce the ED results.For most points of the phase diagram, CPGF has fast convergence and a relatively coarse input resolution suffices to match the results of ED.Even though the convergence behavior of MCLM throughout the phase diagram is not as predictable as CPGF, the convergence speed is comparable to that of CPGF and both typically converge faster than MTPQ.Overall, we find that CPGF requires less computer time and has well controlled accuracy through the resolution and number of polynomials, thus having a slight edge over MCLM.In spite of not being ideal for probing target energies with well controlled accuracy, we still find that MTPQ is very useful because of its ability to estimate the temperature corresponding to each energy density over the course of the iteration.This means that MTPQ is a viable method to carry out studies that would otherwise only be possible using canonical methods.
Regarding finite-temperature studies, we find shortcomings in both MTPQ and its canonical counterpart, CTPQ.Both seem advantageous at first sight due to their lower memory cost.However, in the case of MTPQ, this implies a trade-off that we show to greatly increase the computer time.On the other hand, in the case of CTPQ, a numerical instability limits its ability to probe very low temperatures.We show that the main competitor of TPQ, the Lanczos-based FTLM, has comparatively larger statistical fluctuations when studying the low-temperature behavior of the NN spin correlations.These fluctuations, along with a high number of matrixvector multiplications per iteration in FTLM, limit the method's efficiency.Finally, we show that the newly introduced FTCP method circumvents the shortcomings of the other methods.Its statistical fluctuations are smaller than FTLM and comparable to the TPQ methods.Whilst the FTCP memory cost is the same as FTLM (but slightly higher than TPQ), the trade-off is that the method is more efficient, i.e. twice as fast as MTPQ.This can be a crucial advantage in practical applications.Moreover, when using FTCP, accuracy can be controlled at each iteration, unlike in MTPQ, where the only way to guarantee sufficient accuracy is to increase the upper bound on the maximum eigenvalue of the spectrum by trial-and-error and carry out more costly, longer simulations.This is a demanding process, where one considers that accuracy is sufficient when no changes are detected in the relevant quantities for the desired temperature range as the upper bound is increased and the simulations are repeated.FTCP provides us with the the option of choosing the number of polynomials for convergence in each Chebyshev expansion throughout the iterative process, thereby ensuring that accuracy is maintained in the whole temperature range, without dramatically increasing the computational cost.
Our results show clear trade offs that must be taken into account when choosing which method to use.For example, MTPQ is designed to achieve maximum accuracy for the ground state.However, it cannot isolate excited states nor can it ensure uniform accuracy throughout large temperature windows.Concomitantly, the additional control afforded by the CPGF approach -which can access excited states directly -could be useful for studying non equilibrium systems, such as those studied in Ref. [9].Another example is that of canonical methods.While we find Lanczos to be efficient for the computation of observables closely related to the Hamiltonian, it has large statistical fluctuations for more generic observables that might be of interest, such as the NN spin correlation in the K-H model.
A particularly powerful competing method is DMRG, originally devised to investigate 1D interacting systems [27][28][29].DMRG aims to systematically truncate the exponentially large Hilbert space basis.The basis is rotated in the process in order to improve the accuracy of the truncation.This rotation is achieved via a series of global rotations generated by sweeping the lattice and thus focusing on a few sites at a time.The wide applicability of the DMRG procedure means that it can be used as a general purpose, sign problem free method.Moreover, it yields results that are competitive with QMC.However, the application of DMRG to 2D systems remains challenging [30].The Chebyshev-based methods used throughout this paper are a potential alternative to DMRG because, unlike the latter, they pose no restrictions on boundary conditions and their accuracy can be precisely controlled by ensuring statistical convergence and, in the case of CPGF, by adjusting the spectral resolution.As far as the system size is concerned, Ref. [85] details the use of conservation laws to improve the efficiency of ED methods, particularly from the computer memory point of view.For example, translational symmetry implies that some configurations of the spins are equivalent up to a phase factor.This is a consequence of the block structure of the Hamiltonian, which gives forth to a reduced basis approach, enabling the study of larger system sizes.The only caveat is that systems with open boundary conditions and/or random couplings cannot be tackled with this approach.
To sum up, our results show that Chebyshev methods are more versatile and efficient than their Lanczos and TPQ counterparts, unless one is interested in properties that are well described using solely the ground state, or a small set of low-lying excitations.While in that case, Lanczos is still the method of choice, Chebyshev methods have significant advantages in various other scenarios, namely for the study of: properties that depend on an arbitrary target energy; finite temperature behavior of observables of interest that cannot be expressed in terms of the Hamiltonian and low-order polynomials of the latter; and dynamical quantities, such as spectral functions.

A Stochastic trace evaluation
Statistical convergence is obtained when the error bars become acceptably small.This information is encoded in the scaling of the standard deviation with the number of used initial random states.In Fig. 13, we confirm that we obtain the expected scaling (σ ∝ N −1/2 rd.vec. ) [49], that is our error bars can made as small as required by simply averaging over more realizations of the initial random state.This calculation was carried out with CPGF for a point in the phase diagram of the Kitaev-Heisenberg model with the parametrization of Ref. [3].

B Spectral convergence and computational effort of CPGF
As mentioned in Sec. 4, the spectral convergence of CPGF follows a predictable pattern: the optimal number of polynomials needed for convergence, N * poly , is inversely proportional to the resolution, η.Our results shown in the top panel of Fig. 14 confirm this behavior.In the CPGF, we used a stricter definition of convergence than before: A variation of less than 10 −9 in the energy density between three consecutive iterations.This is necessary because convergence occurs in a "damped oscillatory" manner in CPGF (as we have shown in Fig. 7).
When targeting the ground state, the convergence of the CPGF method slows down close to critical points, despite showing comparably faster convergence in other parts of the phase diagram (center panel of Fig. 14).Surprisingly, MCLM behaves in a complementary fashion: convergence tends to become faster near a phase transition.Furthermore, each iteration is faster overall with CPGF, so even when in cases where both require comparable numbers of iterations for convergence, CPGF tends to be faster.Thus, the CPGF is faster than MCLM at reproducing the complete phase diagram.
The cost of targeting the excited states with CPGF is comparable to targeting the ground state.In contrast, MCLM requires significantly more iterations for convergence, resulting in a .N * poly is approximately proportional to the spectrum width and inversely proportional to the required resolution.Thus, the number of iterations required for convergence can be estimated in advance, unlike in other approaches.Center panel: Computer time required for convergence using MCLM and CPGF to target the ground state.Bottom panel: Similar to the center panel, but now targeting excited states with energy ϵ target .Here, we allowed a maximum of 10000 iterations with MCLM.In comparison, CPGF never required more than 3195 polynomials for convergence.

Figure 2 :
Figure2: Left panel -24-spin hexagonal cluster on the honeycomb lattice.The red, green and purple bonds illustrate the periodic boundary conditions.Right panel -Nearest neighbor bonds colored in red, green and purple (respectively γ = x, y, z in our cartoon).The bond-directional character of the Kitaev interaction implies that to each bond corresponds a different type of interaction.Similarly to the case of the Ising model on the triangular lattice, where we have geometrical frustration, here we have exchange frustration due to the nature of the interaction and it is not possible to find a spin configuration that simultaneously minimizes the energy on all bonds.

Figure 3 :
Figure 3: Minimum (ϵ m ) and maximum (ϵ M ) eigenvalues of the K-H Hamiltonian on a 24-spin hexagonal cluster on the honeycomb lattice, computed with the ground state variant of Lanczos ED.The energy is normalized to the number of spins and the shaded green area represents the spectrum width W = ϵ M − ϵ m for each value of ϕ. ϵ target = ϵ m + 0.1W are target energies used later to compare MCLM and CPGF.

Figure 4 :
Figure4: Top panel -Ground state energy density obtained Lanczos (solid blue line) and MTPQ (dashed orange line).The solid red line is minus the second derivative of the Lanczos ED curve obtained via finite differences.Its peaks signal the quantum critical points labeled "QCP".Bottom panel -Ground state nearest neighbor spin-spin correlation computed with Lanczos ED compared with MTPQ.The correlator is also computed with CPGF for states with target energy ϵ target .The symbol 〈i, j〉 denotes an average over all nearest neighbors 〈i, j〉.

Figure 5 :
Figure 5: Convergence in computations of ground state N-N spin-spin correlations.Top panel -Dependence of N * it with ϕ in Lanczos and MTPQ -(left and right vertical axes, respectively).Bottom panel -Convergence of Lanczos (left) and MTPQ (right) around the Kitaev limits: AFQSL on top and FQSL on the bottom .

Figure 6 :
Figure 6: Comparative study between CPGF and MCLM.Top panel -Ground state energy density obtained with MCLM (solid blue line) and CPGF (dashed lines) with varying resolution as a function of ϕ.The two methods show excellent agreement provided that the CPGF resolution is sufficiently fine.Bottom panel -Ground state N-N spin-spin correlation as a function of ϕ.We used a single realization of the initial random state for both methods.The results shown here also match the Lanczos and MTPQ results shown in Fig. 4.

Figure 7 :
Figure 7: Convergence of the ground state nearest neighbor spin-spin correlation computed with MCLM (left) and CPGF (right) around the Kitaev limits.

s/kB ln 2 Figure 9 :
Figure9: Specific heat (entropy density on the inset) computed with MTPQ, FTCP, CTPQ FTLM.We use 50 realizations of the initial random state for all methods.The error bars are negligibly small in most of the temperature range, except for temperatures below that of the low-temperature peak.The error bars are of comparable size for every method (except for CTPQ, for which their meaning is ill-defined due to numerical instability).Model parameters as in Fig.8.

Figure 10 :
Figure 10: Top and center panels -Negative first (top) and second (center) derivatives of the Lanczos ground state energy density.Bottom panel -Mean square of the magnetization, 〈m 2〉 ≡ 〈[(1/N ) i σ z i ] 2 〉.These results match those of Ref.[61] and illustrate the transitions between ferromagnetic (red), nematic (green) and quantum spin liquid (blue) phases found in Ref.[61].Here, a single realization of the initial random state was found to be sufficient due to negligible statistical fluctuations.

Figure 11 :
Figure 11: MTPQ and FTCP results for the K-I model on a 24-spin cluster with α = 0.7 and J I = 0.001 (left panels) and J I = 0.03 (right panels).We used 50 and 100 random vector realizations for J I = 0.001 and J I = 0.03, respectively.Panels a) and b) show the specific heat, with a zoom-in on the standard deviations -from which the error bars are derived -at low temperature; panels c) and d) show the entropy density, with the grey lines marking s = 0.5k B ln 2; panels e),f) show the finite temperature expectation of the kinetic energy of the Majorana fermions c.

Funding
information F. B. is supported by a DTP studentship funded by the Engineering and Physical Sciences Research Council.A.F. acknowledges the support from the Royal Society through a Royal Society University Research Fellowship.

Figure 13 :
Figure13: Behavior of the standard deviation of the NN spin correlation with the number of realizations of the initial random state for the CPGF method.Here, we consider a specific point of the phase diagram (α = 0.818 in the parametrization of Ref.[3]).We obtain the expected scaling: σ ∝ N −1/2 rd.vec. .

Figure 14 :
Figure 14: Top panel: Number of iterations required for convergence with CPGF times the required resolution (left vertical axis) and spectrum width (right vertical axis).N *poly is approximately proportional to the spectrum width and inversely proportional to the required resolution.Thus, the number of iterations required for convergence can be estimated in advance, unlike in other approaches.Center panel: Computer time required for convergence using MCLM and CPGF to target the ground state.Bottom panel: Similar to the center panel, but now targeting excited states with energy ϵ target .Here, we allowed a maximum of 10000 iterations with MCLM.In comparison, CPGF never required more than 3195 polynomials for convergence.

Table 1 :
Phase boundaries of the K-H model on a 24-spin hexagonal cluster obtained with Lanczos.

Table 2 :
Number of iterations (or Chebyshev polynomials) required for convergence for each method and CPU time per core per iteration.

Table 3 :
Average CPU times, t CPU , for MTPQ and FTCP calculations for the K-I model on a 24-spin hexagonal cluster with α = 0.7 and J I = 0.001, 0.03.