SciPost Submission Page
Realspace spectral simulation of quantum spin models: Application to generalized Kitaev models
by Francisco M. O. Brito, Aires Ferreira
This Submission thread is now published as
Submission summary
Authors (as registered SciPost users):  Francisco Brito 
Submission information  

Preprint Link:  https://arxiv.org/abs/2110.01494v3 (pdf) 
Date accepted:  20240122 
Date submitted:  20231227 13:48 
Submitted by:  Brito, Francisco 
Submitted to:  SciPost Physics Core 
Ontological classification  

Academic field:  Physics 
Specialties: 

Approaches:  Theoretical, Computational 
Abstract
The proliferation of quantum fluctuations and longrange entanglement presents an outstanding challenge for the numerical simulation of interacting spin systems with exotic ground states. Here, we present a toolset of Chebyshev polynomialbased 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 KitaevHeisenberg (KH) and the KitaevIsing (KI) models. We start by applying the Chebyshev toolset to compute nearestneighbor spin correlations, specific heat and entropy of the KH model on a 24spin 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 KI model. Our findings suggest that the efficiency, versatility and lowtemperature stability of the Chebyshev framework developed here could pave the way for previously unattainable studies of quantum spin models in two dimensions.
Author comments upon resubmission
Reply to Referee Report 1
We are pleased that the referee recognizes the improvements in the extended revision of the original manuscript and, in particular, appreciated our detailed analysis of the performance of different methods, the introduction of two new methods (FTCP and HLC), and the expansion on the physics side. We are also glad that the referee showed a particular interest in the energy cutoff dependence of TPQ. We believe that we have implemented the minor changes requested by the referee in this second resubmission:

The limitation on system size stems from the computer memory requirements for the storage of Dvectors, i.e. vectors with the dimension of the Hilbert space; this limitation applies to all the methods discussed in our work (Lanczos, TPQ and Chebyshevbased). Yet, in AIP Conf.Proc.1297:135,2010, the author explains how conservation laws can be used to improve the efficiency of Lanczos exact diagonalization algorithms, in particular from a computer memory standpoint. The arguments presented in the aforementioned reference apply to the Chebyshev approach as well (and to TPQ, for that matter). For example, by using translational symmetry, some configurations of the spins are seen to be equivalent to one another, apart from a phase factor. Implementing this symmetry explicitly in our code (via a reduced basis that reflects the block structure of the Hamiltonian) would enable the study of larger systems. Here, there is an important caveat: this technique does not apply to more general systems, for example with open boundary conditions and/or random couplings. A sentence was added to the paper explaining this technique in the second to last paragraph of Sec. 5: “As far as the system size is concerned…”

We improved the color bar in Fig. 12, so as to clarify the transition to white. We have also added a sentence in the caption explaining that white space indeed corresponds to a vanishing spin susceptibility, as correctly inferred by the referee.

We extended our discussion beyond the ground state, presenting data to back our claim in the last sentence of the first paragraph in Sec. 5. We target excited states (now defined in Fig. 3 in the revised version of the manuscript) and compute the nearest neighbor spinspin correlation for these states as well (see bottom panel of Fig. 4 in the revised version). Although there is already a difference in performance for the ground state, we find it to be more pronounced for excited states, i.e. CPGF is about an order of magnitude faster than MCLM when the excited states are targeted. We attribute this to the uniform convergence of the CPGF, which guarantees that the required Chebyshev iterations remain broadly the same, regardless of the target energy. Differences in convergence speed throughout the phase diagram are explained by the fact that the resolution must be adjusted in order to probe each excited state with sufficient resolution to compute the spin correlator accurately. In turn, finer resolutions require more polynomials.

We augmented our list of references with the suggestions of the referee (see page 3).
Reply to Referee Report 2
We kindly thank the referee for acknowledging our efforts and pointing out the strengths and weaknesses of the first resubmission. We were pleased that the referee: i) valued our pedagogical approach to unbiased EDlike numerical methods; ii) recognized the advantages in computational performance in our Chebyshevbased approach based on the novel methods introduced in this work; iii) recognized our original results for the signatures of quantum phase transitions in the KitaevIsing model obtained from spin dynamics using one of our newly introduced unbiased methods. Below, we address the requested changes, explaining how they were incorporated in our second resubmission, which we believe to overcome the weaknesses of the first resubmission.
1.We achieve remarkable lowtemperature stability with our FTCP approach because of the decomposition in a product of exponentials of Equation (34). The inverse temperature steps are kept small enough that the Chebyshev expansions of Equation (36) remain numerically stable, i.e. the argument of the rapidly increasing modified Bessel functions is controlled. In contrast, in the lowtemperature regime, the CTPQ expansion of Equation (31) involves rapidly increasing factorial functions of the inverse temperature and the number of iterations. At low temperatures, Equation (31) can no longer be easily evaluated with enough accuracy because the limits of machine precision are reached and eventually it generates overflows. This gradual loss of accuracy becomes visible for the CTPQ results as the temperature is lowered. Eventually, the specific heat diverges, rather than going to zero, as expected. This is inconsistent with the other methods and is merely a numerical artefact due to loss of precision. We added a sentence in Sec. 4.1.2., explaining why FTCP avoids the problem in evaluating the specific heat in contrast to CTPQ: “The operator breakup…”

We agree with the referee. These references now appear in the revised manuscript. The Reference Phys. Rev. B 92, 214431 (2015) is cited in Sec. 3.3.3., in the sentence “We could expand the operator…”. The Reference Phys. Rev. Lett. 121, 220601 (2018) is cited in Sec. 3.4: “On the other hand, the Chebyshev approximation of the time evolution operator  which is used e.g. in Ref. Phys. Rev. Lett. 121, 220601 (2018) in combination with CTPQ…”. Unlike Phys. Rev. B 92, 214431 (2015), in our finite temperature Chebyshev polynomial (FTCP) approach, we do not Chebyshevexpand the operator exp(beta H/2) directly. Instead, we first write exp(beta H/2) as a product over L inverse temperature steps in Eq. (34). Then, we Chebyshevexpand each term of the product. This procedure allows us to maintain numerical stability because the arguments of the fast growing modified Bessel functions of Eq. (35) are kept small. If this is not done, we run into overflows when probing low temperatures. In Phys. Rev. B 92, 214431 (2015), the temperature is kept high enough that this never becomes a problem. For the systems we study in our work, lowtemperature features are of particular interest, which motivated us to develop our FTCP method. We added a sentence in Sec. 3.3.3. explaining these differences with respect to the method used in Phys. Rev. B 92, 214431 (2015): “We could expand the operator exp(beta h / 2) in Chebyshev polynomials directly as done e.g. in Phys. Rev. B 92, 214431 (2015). However, such an expansion is vulnerable to numerical instabilities for large beta (low temperature) due to the rapid growth of the Bessel functions.” In Phys. Rev. Lett. 121, 220601 (2018), the authors use Chebyshev expansions to approximate the time evolution operator. The authors remark that in order to ensure convergence of the spectral functions in their approach, it is necessary to take a long enough time window. Additionally, it is important to consider the time step and the number of polynomials. In contrast, our approach does not rely on timedomain correlators. Instead, we use a Chebyshev expansion to compute the frequencydomain spectral functions directly. The main advantage of our approach is that convergence is set via the number of polynomials. There is no need to worry about a time window, or a time step. Moreover, the authors of Phys. Rev. Lett. 121, 220601 (2018) mention using up to 8000 time steps, with Chebyshev expansions up to order 500. This suggests that our method might be more efficient, possibly due to treating spectral functions directly in the frequency domain. Our approach involves only a single Chebyshev expansion with 3000 polynomials to resolve the spectral function, which is roughly equivalent to 6 time steps using 500 polynomials in Phys. Rev. Lett. 121, 220601 (2018).

We chose the ground state to benchmark CPGF simply to check consistency with Lanczos ED results. We agree with the referee that Lanczos is the method of choice to study the ground state. Thus, we have expanded our comparison between MCLM and CPGF beyond the ground state. We now target excited states (as explained in Figure 3 in the revised version of the manuscript). In particular, we compare the performance of MCLM and CPGF in a computation of the nearest neighbor spinspin correlation, similarly to what had been done for the ground state in the previous submission. MCLM and CPGF results are consistent. In the bottom panel of Figure 4 in the revised version, we show the CPGF results, which match MCLM results. Finally, we added a panel to Figure 14, showing our comparison between the computer time required to target the excited states with the two methods. Although CPGF was already faster for the ground state, the difference in performance is more dramatic for excited states (CPGF is about an order of magnitude faster than MCLM). We conclude that CPGF is the method of choice for targeting excited states because of its control over resolution and its advantageous convergence properties.
List of changes
The references suggested by the first referee were added in page 3 after the sentence: “The severity of the problem depends on the computational basis…”
A sentence was added in Sec. 3.3.3. to clarify an aspect of the FTCP method as requested by Referee 2: “We could expand the operator exp(beta h / 2) in Chebyshev polynomials directly as done e.g. in Phys. Rev. B 92, 214431 (2015). However, such an expansion is vulnerable to numerical instabilities for large beta due to the rapid growth of the Bessel functions.”
A sentence was slightly modified in Sec. 3.4.: “On the other hand, the Chebyshev approximation of the time evolution operator  which is used e.g. in Phys. Rev. Lett. 121, 220601 (2018) in combination with CTPQ  relies on…”. We also added the following footnote “In Phys. Rev. Lett. 121, 220601 (2018), CTPQ is used to generate an initial thermal state at t=0. Time evolution is carried out using the Chebyshev approximation of exp(i t H)”.
Figure 3 has been changed slightly so as to show the target energy used to compare MCLM and CPGF. The caption has also been changed to reflect this.
The bottom panel of Figure 4 now shows CPGF results targeting an excited state. MCLM results match CPGF and thus they are not shown so as not to overcrowd the figure. The caption and legend have also been updated to explain the extra data. A sentence has been added to comment on the excited state data: “Also shown in the bottom panel of Fig. 4…”
A sentence has been added in Sec. 4.1.1. after Figure 6, commenting on the advantages of using CPGF to target the excited states: “Nevertheless, when targeting the ground state, CPGF is 25% faster on average. For the excited states…”
A sentence was added in Sec. 4.1.2., reinforcing the numerical stability of FTCP in contrast to CTQP: “The operator breakup…”
The color bar in Figure 12 has been adapted. The white space indeed means that the spectral function vanishes, as the referee pointed out. This is now evident in the bottom of the color bar, where the blue color is smoothly changed to white near 0. The caption has also been extended to include an explanation of the color bar.
A sentence was added to the second to last paragraph of Sec. 5, explaining that larger systems can be tackled with EDlike methods by taking advantage of the symmetries of the problem at hand.
A new panel has been added to Figure 14 with a comparison of the computer time required to target excited states with MCLM and CPGF. The caption has been changed accordingly and the final paragraph of the appendix has also been expanded. We find that CPGF significantly outperforms MCLM. In fact, with CPGF, the overall computational cost of targeting the excited states is comparable to that of targeting the ground state. With MCLM, the cost is larger by about an order of magnitude in comparison to CPGF.
Published as SciPost Phys. Core 7, 006 (2024)