Reorthonormalization of Chebyshev matrix product states for dynamical correlation functions
Abstract
The Chebyshev expansion offers a numerically efficient and easyimplement algorithm for evaluating dynamic correlation functions using matrix product states (MPS). In this approach, each recursively generated Chebyshev vector is approximately represented by an MPS. However, the recurrence relations of Chebyshev polynomials are broken by the approximation, leading to an error which is accumulated with the increase of the order of expansion. Here we propose a reorthonormalization approach to remove this error introduced in the loss of orthogonality of the Chebyshev polynomials. Our approach, as illustrated by comparison with the exact results for the onedimensional XY and Heisenberg models, improves significantly the accuracy in the calculation of dynamical correlation functions.
I Introduction
Dynamical correlation functions, such as optical conductivity or singleparticle spectral function, are measurement quantities and are of central importance in condensed matter physics. However, calculation of these quantities is a challenging theoretical problem. During the last two decades, several approaches have been proposed to calculate spectral functions using the density matrix renormalization group (DMRG) White (1992) or other methods related to tensor network states Östlund and Rommer (1995); Verstraete and Cirac (2004); Xie et al. (2014). In 1987, Gagliano and Baliseiro Gagliano and Balseiro (1987) (1987) proposed a continuedfraction method (also called the Lanczos vector method) to evaluate dynamic correlation functions based on the Lanczos diagonalization. Their method was first adopted in the DMRG calculation by Hallberg Hallberg (1995) in 1995. While this method requires only modest numerical resources, reliable results can be obtained only limited to low frequencies. In 1989, a correction vector method was proposed by Soos and Ramasesha Soos and Ramasesha (1989) to improve the continuedfraction method. This correction vector method was adopted in the DMRG calculation first by Ramasesha et al. Pati et al. (1997) (1997) and later by Kühner and White Kühner and White (1999) (1999). The key idea is to take the Green’s function at a particular frequency as a target state and calculate it by solving a set of large but sparse linear equations. This method can generate very accurate result at any given frequency, but its computational cost is very high. In 2002, Jeckelmann Jeckelmann (2002) improved the correction vector method and showed that it is much more efficient and accurate to determine the correction vector by minimizing a cost function. Recently, Dargel et al. proposed an adaptive Lanczos vector method Dargel et al. (2011) (2011) and an MPSbased Lanczos method Dargel et al. (2012) (2012) that invokes the idea of orthogonalization to improve the accuracy of dynamic correlation functions. More recently, Nocera and Alvarez Nocera and Alvarez (2016) (2016) proposed a Krylovspace approach to replace the conjugate gradient method in the calculation of correction vectors.
Dynamic quantities can be also calculated from the Fourier transformation of the timedependent Green’s function. However, to obtain good frequency resolution, one has to calculate the correlation functions over a long time interval, which is limited by either a loss of accuracy due to the approximation used or by finitesize effects such as reflections from open ends. In 2002, Cazalilla and Marston evaluated the timedependent correlation functions for a onedimensional (1D) system under an applied bias using the DMRG Cazalilla and Marston (2002). In 2003, Luo et al. Luo et al. (2003) proposed a packkeeping approach to optimize the basis states that are retained with the time evolution. In 2004, Vidal introduced the socall timeevolving blockdecimation (TEBD) method to efficiently simulate timeevolution of 1D systems with shortrange interactions using MPS Vidal (2004). Based on this approach, an adaptive timedependent DMRG method was explored by Daley et al. Daley et al. (2004) and by White and Feiguin White and Feiguin (2004).
Longtime dynamics, however, remains challenging due to the linear growth of the entanglement entropy as the state evolves in time Calabrese and Cardy (2005), which implies that an exponentially growing bond dimension is required in order to keep the error fixed. In 2009, Bañuls et al. Bañuls et al. (2009) proposed a folding scheme to present more efficiently the entanglement structure of the tensornetwork state in the time direction, providing an accurate method for evaluating longtime dynamics of the ground state in the thermodynamic limit. A quantum transfermatrix DMRG extension of this method to finite temperatures was proposed by Huang et al. Huang et al. (2014). Generalized DMRG approaches have also been extended to simulate time evolution at finite temperatures by introducing auxiliary degrees of freedom to purify the thermal statistical operators Barthel et al. (2009); Karrasch et al. (2012).
Recently, Holzner et al. Holzner et al. (2011) (2011) proposed a Chebyshev matrixproductstates (CheMPS) method to calculate dynamical spectral functions for operators and at zero temperature
(1) 
where is the Hamiltonian, is the ground state energy and is the corresponding eigenstate. The main idea is to use the Chebyshev polynomials to expand the delta function and to use a set of MPS to represent the Chebyshev vectors generated in the Chebyshev expansion. This method offers precise and convenient control of the accuracy and resolution. It yields results with accuracies comparable to those of the correctionvector DMRG, but at dramatically reduced numerical cost.
The CheMPS provides a balanced scheme between cost and accuracy. It resolves spectral functions well if the spectral weight is concentrated or limited to low energy part, but cannot resolve the highenergy spectrum accurately, similar as in the continuedfraction DMRG. Since each Chebyshev vector is approximated by an MPS at each step of the Chebyshev expansion, the recurrence relation between different Chebyshev vectors is satisfied only approximately. Thus the resolution may not be truly improved with the increase of the series number of Chebyshev expansion.
In this paper, we introduce a reorthonormalization scheme to improve the accuracy in the calculation of spectral function based on the CheMPS. (Similar approach was used to improve the LanczosMPS calculation Dargel et al. (2012)). In the discussion below, for simplicity, we refer our method as ReCheMPS. In this method, the Chebyshev expansion is applied twice. We first carry out a CheMPS calculation to obtain a set of Chebyshev vectors represented using MPS, and then reorthonormalize these MPS to obtain a set of manybody basis states Huang et al. (2016). Within the truncated Hilbert space spanned by this set of basis states, we rediagonalize the Hamiltonian and evaluate the spectral function, using the Chebyshev expansion again. This method, as shown below, improves the accuracy of CheMPS significantly, with very limited increase of computational cost in comparison with CheMPS. Furthermore, we propose a simple approach to extrapolate the finitesize results to the thermodynamic limit using the eigenenergies of the truncated Hamiltonian and the corresponding spectral weights. By comparison with the exact results in the thermodynamic limit, we find that this approach offers an efficient way to perform the finitesize scaling in the entire energy range.
This paper is organized as follows. In Sec. II, after a brief review of the CheMPS in Sec. II.1, the ReCheMPS is introduced in Sec. II.2. We propose a finite size scaling approach in Sec II.3, based on the eigenspectra and the spectral weights obtained from the effective Hamiltonian in the orgonalized but truncated basis space. Then we calculate the magnetic structure factor using the ReCheMPS for the 1D XY and Heisenberg models. The results are presented in Sec. III.1 for the XY model and in Sec. III.2 for the Heisenberg model, and compared with those obtained by the CheMPS. The finite size scaling result is showen in Sec. III.3 and compared with the exact results obtained in the thermodynamic limit. Finally, Sec. IV gives a summary.
Ii Method
ii.1 The CheMPS method
To carry out a CheMPS calculation Holzner et al. (2011), one needs to first calculate the wave function of the ground state and the corresponding energy using the DMRG or variational MPS methods. In order to determine the spectral function defined in Eq. (1) in a frequency interval , we rescale and shift the frequency and the Hamiltonian such that the nonzero parts of spectral function are completely concentrated on the frequency interval , where is usually taken to be slightly smaller than 1 for numerical stability Holzner et al. (2011). Then the spectral function becomes
(2) 
where , and
(3) 
The function in Eq. (2) can be expanded by an infinite Chebyshev polynomial series Weiße et al. (2006). In practical calculation, this expansion is approximated by the first terms of the Chebyshev polynomials Holzner et al. (2011), i.e.
(4) 
where and are the damping factors to suppress the Gibbs oscillations introduced by the truncation to the infinite Chebyshev series. There are several different choices for these damping factors Weiße et al. (2006), here we use the Jackson damping defined by
(5)  
In Eq. (4), are the Chebyshev polynomials defined by the recurrence relation
(6) 
with and . can be also represented as
(7) 
Substituting Eq. (4) into Eq. (2), we obtain
(8) 
where are the Chebyshev moments defined by
(9) 
To obtain the Chebyshev moments , we first calculate the Chebyshev vectors
(10) 
From the recurrence relation Eq. (6), we find that , , and
(11) 
These Chebyshev vectors span a Krylov basis space.
In the CheMPS calculation, the ground state and all the Chebyshev vectors are approximately represented by MPS. In this case, the above equations are valid just approximately. The MPS representation of the first two Chebyshev vectors are determined by variationally minimizing the functions
(12)  
(13) 
Similarly, other Chebyshev vectors are determined by minimizing the cost function
(14) 
From these Chebyshev vectors, the Chebyshev moments are determined by the formula
(15) 
Given a Chebyshev expansion number , it is known that the CheMPS has a spectral resolution of order in the interval Holzner et al. (2011). When the bandwidth is broadened or the spectrum changes drastically, one has to increase to increase the resolution.
ii.2 The ReCheMPS Method
Two approximations are used in the CheMPS. The first is to truncate the infinite Chebyshev series and approximate it by a finite polynomial. The error such introduced can be reduced by increasing the number of Chebyshev polynominals . The second is to approximate each Chebyshev vector by an MPS. This is a more severe approximation because it breaks the orthonormal condition of Chebyshev series, and the error is accumulated in the calculation of Chebyshew vectors with Eq. (14).
Here we propose a reorthonormalization scheme to solve the second problem. We start from the Chebyshev vectors obtained using the CheMPS, and then carry out the following reorthonormalization calculation:

Reorthonormalize the Chebyshev vectors using the GramSchmidt orthogonalization process
(16) so that
(17) Here with are the new orthonormalized basis states. are the normalization constants.

Represent the Hamiltonian and other operators in the basis space spanned by ,
(18) 
Redo the Chebyshev expansion and generate new Chebyshev vectors using in this new basis space.
In this new scheme, the Chebyshev expansions are applied twice. In the first expansion, approximate Chebyshev vectors represented by MPS are evaluated by variationally minimizing the cost functions defined in Eqs. (12  14). These MPS are not orthonormalized due to the truncation in the bond dimension.
In the second expansion, the Hamiltonian is approximately replaced by , which is an matrix defined in the space spanned by , and Chebyshev vectors that satisfy the standard recurrence relations of Chebyshev polynominals are generated using . Here is the number of Chebyshev vectors used for approximating the function using Eq. (4) for the Hamiltonian . It can be set either equal to or as an independent control parameter. The Chebyshev expansion Eq. (4) becomes exact in the limit .
In the calculation of ReCheMPS, only the inner products between any two MPS already generated in the CheMPS need to be calculated. Its cost in both computational time and memory is quite small in comparison with the calculation of CheMPS. Thus, as shown below, our reorthogonalization scheme can improve significantly the accuracy of the CheMPS, but with very limited increase in the computational time.
The Hamiltonian has a finite dimension. We can diagonalize it to obtain its spectrum
(19) 
where and are the eigenpairs of , which are arranged in an ascending order, i.e. if . In this truncated Hilbert space, the spectral function is given by
(20) 
In the second Chebyshev expansion, the above deltafunctions are blurred by Chebyshev polynomials. In this round of expansion, there is no technical difficulty to take arbitrary many Chebyshev vectors in the expansion. Therefore, we can set equal or not equal to , depending on the resolution required. Eq. (20) is recovered in the limit .
Furthermore, as the spectrum of is known, we can always set the bandwidth of , represented by , exactly equal to the difference between the largest and smallest eigenvalues of , i.e. , and find all the Chebyshev vectors without taking any approximation in the second Chebyshev expansion. Clearly, , thus if the same number of Chebyshev vectors is used as in CheMPS, i.e. setting , the energy resolution, which is qualitatively determined by the ratio , should be higher in ReCheMPS.
ii.3 Finitesize scaling
The CheMPS offers a convenient way to control the accuracy and resolution of the spectral function by simply adjusting the expansion order . This is an advantage that can be used to analyze the finitesize effects Holzner et al. (2011). For example, one can mimic the infinite lattice limit () by choosing a small enough or more precisely a large to smear out the finitesize subpeaks. However, in practical calculation, it is difficult to adjust to reduce the finite size effect for the systems with different lattice lengths.
An alternative approach is to take a large limit and assume that all spectral peaks can be reliably resolved. By fitting the peaks by Gaussian functions, one can obtain the central position and the weight for each peak. The spectral function at each peak position is then defined by the ratio between the peak weight and the half distance between the two neighboring peaks. This approach does yield quite impressive results for the Heisenberg model Holzner et al. (2011). However, it relies on the accurate calculation of spectral peaks, which is in fact difficult due to the MPS approximation of Chebyshev vectors. Technically, it is also difficult to use a large to do the Chebyshev expansion when the lattice size or the bond dimension becomes large.
Here we propose a novel approach to perform the finitesize scaling. In Sec. II.2, we show that by diagonalizing the Hamiltonian in the truncated basis space, we can approximately represent the spectral function as
(21) 
where and are the eigenpairs of , defined in Eq. (19), and
(22) 
is the spectral weight corresponding to .
Let us denote as the midpoint between and , i.e.
(23) 
The average spectral weight within the interval between and () is given by
(24) 
where
(25) 
is the energy difference between and , and
(26) 
is the average energy of these two points. Here should be taken such that the interval between and contains just one main peak.
Iii Results
In this section, we present the results obtained by the ReCheMPS and compare them with those obtained by the CheMPS for both the XY and the antiferromagnetic Heisenberg models in one dimension. These models are respectively defined by the Hamiltonians
(27) 
and
(28) 
where is the spin operator and is the lattice length. For both cases, open boundary conditions are assumed.
We have calculated the spinspin correlation function defined by
(29) 
where is the spin operator defined in momentum space
(30) 
and is the quasimomentum.
In the calculation of this correlation function, the energy spectrum is bounded by the band width , which equals and for the 1D XY and Heisenberg models, respectively. However, in the CheMPS calculation, this energy bound is broken by the approximation used in the MPS representation of Chebyshev vectors, and the effective band width can be larger than the exact one. To avoid this overflow of energy, we set equal to a value which is greater than the exact one so that the energy spectrum can be fully covered in the first Chebyshev expansion. In this work, we take for both models. By doing this, we do not need to employ the energy truncation method used in Ref. [Holzner et al., 2011].
During the second Chebyshev expansion, as already mentioned in Sec. II.2, we take the bandwidth exactly equal to the difference between the largest and the smallest eigenenergies, i.e. , of the approximate Hamiltonian . For all the calculation we have done, we find that and varies slightly with . The results presented below are obtained by taking if not particularly specified.
iii.1 The XYmodel
The 1D XY model can be mapped onto a free spinless fermion model through the JordanWigner transformation Lieb et al. (1961). The spinspin correlation function contributes mainly by the particlehole excitations of spinless fermions, whose energy is bounded respectively from below and above by
(31) 
As , the spinspin correlation function can be exactly represented as
(32) 
where . In the limit , it becomes
(33) 
Figure 1(a) shows as a function of frequency for the 1D XY model obtained with ReCheMPS on the lattice. As expected, the spectral resolution is improved with the increase of . The peak positions of , as illustrated in Fig. 1(b), converge rapidly to the exact energies, , again with the increase of .
Figure 2 compares the result obtained by ReCheMPS with that by CheMPS. The improvement of ReCheMPS over CheMPS is quite significant. By keeping the same number of Chebyshev vectors, the spectral resolution of ReCheMPS is clearly better than that of CheMPS. Besides, the ReCheMPS can also improve the accuracy of the peak positions, in comparison with the exact ones, especially for the intermediate and high energy peaks.
Moreover, this improvement becomes more pronounced when the system size becomes larger. Fig. 3 shows how the spectral function changes with for the system of . By comparison with the result shown in Fig. 2, we find that a larger is needed in order to obtain the same accuracy in the spectral function when the lattice size is increased, while the result obtained by the CheMPS does not show significant improvement with the increase of D because of the violation of the orthonormal condition of Chebyshev series. The peak positions we obtain with the ReCheMPS agree well with the exact ones, especially for the case .
With the further increase of the lattice size, the computational cost (including both computation time and memory) increases rapidly with increasing , and it is difficult to use a large to do the Chebyshev expansion. This limits the resolution of the CheMPS calculation. However, in the ReCheMPS calculation, we can use a large to increase the resolution. Fig. 4 shows, as an example, how the resolution changes with the increase of in ReCheMPS for the XY model with and . When , the resolutions of both ReCheMPS and CheMPS are very poor. However, when we increase to 500, the resolution is greatly improved and the ReCheMPS produces all the low energy peaks accurately.
iii.2 The Heisenberg model
For the spin1/2 Heisenberg model, the asymptotic behavior of the dynamical spin structure factor is known from the Bethe ansatz solution Caux and Hagemans (2006). It diverges at in the limit and
(34) 
Since the spectra weight of the Heisenberg model diverges in the zero energy limit, we find that the low energy spectrum can be accurately calculated with relatively small in a large lattice system.
Figure 5(a) compares the results obtained using CheMPS and ReCheMPS with for the Heisenberg model of . The spectral peak obtained with the ReCheMPS is higher because the ReCheMPS has higher resolution than the CheMPS. Moreover, the CheMPS peaks drop quickly with the increase of energy and become almost invisible at the high energy end. But the high energy peaks obtained with the ReCheMPS can be seen in the whole energy range.
The peak energies obtained with these two methods agree well with each other in the lowenergy limit. But the difference becomes larger when the energy is increased. Similar as for the XY model, we find that the low energy spectral peaks obtained with ReCheMPS converge quickly with the increase of . The accuracy in the highenergy spectral function, as shown in Fig. 5(b),can be improved by increasing from 64 to 256. As expected, the spectral weight drops quickly with increasing energy. The energy dependence of the peak height, which can be more clearly seen from the loglog plot of shown in the inset of Fig. 5(b), follows the trend that is expected from the Bathe Ansatz solution Eq. (34) obtained in the thermodynamic limit.
We also calculate the spin structure factor for the 1D Heisenberg model in the whole momentum space. An intensity plot of is shown in Fig. 6. can be measured by inelastic neutron scattering spectroscopy. Our ReCheMPS calculation produces correctly the behavior of spinon continuum spectra. The lower and upper bounds of the continuum agree well the exact results, i.e.
(35) 
iii.3 Finitesize scaling
Here we use the method introduced in Sec. II.3 to perform a finite size scaling analysis for the spectra obtained using the effective Hamiltonian . In practical calculation of the spectral function using Eq. (24), the energy interval should be taken such that it is roughly equal to the energy difference between two neighboring main peaks in the spectral function. For example, for the spectral shown in Fig. 5(b), we can clearly identify ten peaks in the energy interval from 0 to 1. In some cases, it is sufficient just to take . But in the following cases, one should enlarge the value of so that contains more than one energy level: (1) two or more energy levels are nearly degenerate, such that their energy difference is significantly smaller than the main peak distance; and (2) add one more energy level to does not change significantly the total spectral weight within this enlarged energy interval. Practically, one can start with , and then add one or more energy level in. If the change in the total spectral weight within the energy interval before and after adding the energy levels is smaller than a threshold, one can then augment to contain all those energy levels, provided that is not significantly larger than the peak distance. The threshold can be set, for example, from 0.1 to a few percent of the total spectral weight depending on the accuracy of calculated data.
We have calculated the dynamical spinspin correlation functions using Eq. (24) for the XY and Heisenberg models on finite lattice systems. The results are depicted in Figs. 7 and 8 for these two models, and compared with the exact results obtained in the limit . For both models, we find that the results we obtain on the finite size lattices using the approach introduced above agree well in the entire energy range with the exact ones obtained in the thermodynamic limit. Thus our approach provides an accurate and reliable tool to perform the finite size scaling.
Iv Summary
In summary, we propose an orthogonalization scheme to improve the accuracy in the calculation of dynamical spectra functions using the MPS in the Chebyshev expansion. We first perform a CheMPS calculation to obtain a set of Chebyshev vectors represented by MPS. These recursively generated Chebyshev vectors provide a set of truncated basis states to resolve the spectral function over the entire spectral width uniformly. However, these vectors do not truly satisfy the recurrence relations of the Chebyshev polynomials due to the approximation used in the CheMPS calculation. Our approach is to reorthonormalize these vectors to obtain a set of orthonormalized basis states. The matrix elements of the Hamiltonian are evaluated in this truncated basis space, which defines an effective Hamiltonian. We then calculate the spectral function using this effective Hamiltonian. Again the Chebyshev expansion is used to resolve the spectral peaks in a finite lattice system. Since the effective Hamiltonian is diagonalizable and the Chebyshev polynomials can be determined to an arbitrarily high order, this scheme provides a highresolution solution to the calculation of dynamic correlation functions. Using this approach, we have calculated the spectral functions using our method for the 1D XY and Heisenberg models. By comparison with the results obtained using CheMPS, we find that our method improves greatly both the accuracy and the resolution of dynamical spectra functions.
We can eliminate the broadening effect by taking a discrete representation of the spectral function using directly the eigenenergies of the effective Hamiltonian and the corresponding spectral weight. The value of spectral function is determined by the ratio between the spectral weight and the corresponding energy width for each eigenenergy. This can reduce the finite size effect and provide a reliable and accurate approach to scale the finitelattice results to an infinite lattice. Our results obtained with this approach, as shown in Figs. 7 and 8, agree very well with the exact results obtained in the thermodynamic limit.
Acknowledgments
This work was supported by the National R&D Program of China (Grant No. 2017YFA0302901) and the National Natural Science Foundation of China (Grants No. 11190024 and No. 11474331).
References
 S. R. White, “Density matrix formulation for quantum renormalization groups,” Phys. Rev. Lett. 69, 2863–2866 (1992).
 S. Östlund and S. Rommer, “Thermodynamic limit of density matrix renormalization,” Phys. Rev. Lett. 75, 3537–3540 (1995).
 F. Verstraete and J. I. Cirac, “Renormalization algorithms for quantummany body systems in two and higher dimensions,” arXiv:condmat/040706 (2004).
 Z. Y. Xie, J. Chen, J. F. Yu, X. Kong, B. Normand, and T. Xiang, “Tensor renormalization of quantum manybody systems using projected entangled simplex states,” Phys. Rev. X 4, 011025 (2014).
 E. R. Gagliano and C. A. Balseiro, “Dynamical properties of quantum manybody systems at zero temperature,” Phys. Rev. Lett. 59, 2999–3002 (1987).
 K. A. Hallberg, “Densitymatrix algorithm for the calculation of dynamical properties of lowdimensional systems,” Phys. Rev. B 52, R9827–R9830 (1995).
 Z. G. Soos and S. Ramasesha, “Valence Bond Approach to Exact Nonlinear OpticalProperties of Conjugated Systems,” J. Chem. Phys. 90, 1067–1076 (1989).
 S. K. Pati, S. Ramasesha, and D. Sen, “Lowlying excited states and lowtemperature properties of an alternating spin1˘spin1/2 chain: A densitymatrix renormalizationgroup study,” Phys. Rev. B 55, 8894–8904 (1997).
 T. D. Kühner and S. R. White, “Dynamical correlation functions using the density matrix renormalization group,” Phys. Rev. B 60, 335–343 (1999).
 E. Jeckelmann, “Dynamical densitymatrix renormalizationgroup method,” Phys. Rev. B 66, 045114 (2002).
 P. E. Dargel, A. Honecker, R. Peters, R. M. Noack, and T. Pruschke, “Adaptive lanczosvector method for dynamic properties within the density matrix renormalization group,” Phys. Rev. B 83, 161104 (2011).
 P. E. Dargel, A. Wöllert, A. Honecker, I. P. McCulloch, U. Schollwöck, and T. Pruschke, “Lanczos algorithm with matrix product states for dynamical correlation functions,” Phys. Rev. B 85, 205119 (2012).
 A. Nocera and G. Alvarez, “Spectral functions with the density matrix renormalization group: Krylovspace approach for correction vectors,” Phys. Rev. E 94, 053308 (2016).
 M. A. Cazalilla and J. B. Marston, “Timedependent densitymatrix renormalization group: A systematic method for the study of quantum manybody outofequilibrium systems,” Phys. Rev. Lett. 88, 256403 (2002).
 H. G. Luo, T. Xiang, and X. Q. Wang, “Comment on “timedependent densitymatrix renormalization group: A systematic method for the study of quantum manybody outofequilibrium systems”,” Phys. Rev. Lett. 91, 049701 (2003).
 G. Vidal, “Efficient simulation of onedimensional quantum manybody systems,” Phys. Rev. Lett. 93, 040502 (2004).
 A. J. Daley, C. Kollath, U. Schollwock, and G. Vidal, “Timedependent densitymatrix renormalizationgroup using adaptive effective Hilbert spaces,” J. Stat. Mech.Theory Exp. (2004), 10.1088/17425468/2004/04/P04005.
 S. R. White and A. E. Feiguin, “Realtime evolution using the density matrix renormalization group,” Phys. Rev. Lett. 93, 076401 (2004).
 P. Calabrese and J. Cardy, “Evolution of entanglement entropy in onedimensional systems,” J. Stat. Mech.Theory Exp. 2005, P04010 (2005).
 M. C. Bañuls, M. B. Hastings, F. Verstraete, and J. I. Cirac, “Matrix product states for dynamical simulation of infinite chains,” Phys. Rev. Lett. 102, 240603 (2009).
 Y.K. Huang, P. Chen, Y.J. Kao, and T. Xiang, “Longtime dynamics of quantum chains: Transfermatrix renormalization group and entanglement of the maximal eigenvector,” Phys. Rev. B 89, 201102 (2014).
 T. Barthel, U. Schollwöck, and S. R. White, “Spectral functions in onedimensional quantum systems at finite temperature using the density matrix renormalization group,” Phys. Rev. B 79, 245101 (2009).
 C. Karrasch, J. H. Bardarson, and J. E. Moore, “Finitetemperature dynamical density matrix renormalization group and the drude weight of spin chains,” Phys. Rev. Lett. 108, 227206 (2012).
 A. Holzner, A. Weichselbaum, I. P. McCulloch, U. Schollwöck, and J. von Delft, “Chebyshev matrix product state approach for spectral functions,” Phys. Rev. B 83, 195115 (2011).
 R. Z. Huang, H. J. Liao, Z. Y. Liu, H. D. Xie, Z. Y. Xie, H. H. Zhao, J. Chen, and T. Xiang, “A generalized lanczos method for systematic optimization of tensor network states,” arXiv:1611.09574 (2016).
 A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, “The kernel polynomial method,” Rev. Mod. Phys. 78, 275–306 (2006).
 E. Lieb, T. Schultz, and D. Mattis, “Two soluble models of an antiferromagnetic chain,” Annals of Physics 16, 407 – 466 (1961).
 J.S. Caux and R. Hagemans, “The fourspinon dynamical structure factor of the heisenberg chain,” J. Stat. Mech.Theory Exp. 2006, P12013 (2006).