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

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 re-submission:

1. The limitation on system size stems from the computer memory requirements for the storage of D-vectors, i.e. vectors with the dimension of the Hilbert space; this limitation applies to all the methods discussed in our work (Lanczos, TPQ and Chebyshev-based). 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…”

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

3. 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 spin-spin 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.

4. 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 re-submission. We were pleased that the referee: i) valued our pedagogical approach to unbiased ED-like numerical methods; ii) recognized the advantages in computational performance in our Chebyshev-based approach based on the novel methods introduced in this work; iii) recognized our original results for the signatures of quantum phase transitions in the Kitaev-Ising 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 low-temperature 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 low-temperature 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 break-up…”

1. 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 Chebyshev-expand 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 Chebyshev-expand 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, low-temperature 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 time-domain correlators. Instead, we use a Chebyshev expansion to compute the frequency-domain 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).

2. We chose the ground state to bench-mark 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 spin-spin 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.

-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 break-up…”
-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 ED-like 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.