{'fraction_non_alphanumeric': 0.029633274041554705, 'fraction_numerical': 0.009726374749271467, 'mean_word_length': 4.495840266222961, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reorthonormalization of Chebyshev matrix product states for dynamical correlation functions

The Chebyshev expansion offers a numerically efficient and easy-implement 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 one-dimensional 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 single-particle 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) 1 or other methods related to tensor network states [2][3][4] . In 1987, Gagliano and Baliseiro 5 (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 6 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 7 to improve the continuedfraction method. This correction vector method was adopted in the DMRG calculation first by Ramasesha et al. 8 (1997) and later by Kühner and White 9 (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 10 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 11 (2011) and an MPS-based Lanczos method 12 (2012) that invokes the idea of orthogonalization to improve the accuracy of dynamic correlation functions. More recently, Nocera and Alvarez 13 (2016) proposed a Krylov-space 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 time-dependent 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 finite-size effects such as reflections from open ends. In 2002, Cazalilla and Marston evaluated the time-dependent correlation functions for a one-dimensional (1D) system under an applied bias using the DMRG 14 . In 2003, Luo et al. 15 proposed a pack-keeping approach to optimize the basis states that are retained with the time evolution. In 2004, Vidal introduced the so-call time-evolving blockdecimation (TEBD) method to efficiently simulate timeevolution of 1D systems with short-range interactions using MPS 16 . Based on this approach, an adaptive timedependent DMRG method was explored by Daley et al. 17 and by White and Feiguin 18 .
Long-time dynamics, however, remains challenging due to the linear growth of the entanglement entropy as the state evolves in time 19 , which implies that an exponentially growing bond dimension is required in order to keep the error fixed. In 2009, Bañuls et al. 20 proposed a folding scheme to present more efficiently the entanglement structure of the tensor-network state in the time direction, providing an accurate method for evaluating long-time dynamics of the ground state in the thermodynamic limit. A quantum transfer-matrix DMRG extension of this method to finite temperatures was proposed by Huang et al. 21 . 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 22,23 .
Recently, Holzner et al. 24 (2011) proposed a Chebyshev matrix-product-states (CheMPS) method to calculate dynamical spectral functions for operators A and B at zero temperature where H is the Hamiltonian, E 0 is the ground state energy and |0 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 correction-vector 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 high-energy spectrum accurately, similar as in the continued-fraction 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 Lanczos-MPS calculation 12 ). 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 many-body basis states 25 . Within the truncated Hilbert space spanned by this set of basis states, we re-diagonalize 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 finite-size results to the thermodynamic limit using the eigen-energies 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 A, the ReCheMPS is introduced in Sec. II B. We propose a finite size scaling approach in Sec II C, based on the eigen-spectra 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 A for the XY model and in Sec. III B for the Heisenberg model, and compared with those obtained by the CheMPS. The finite size scaling result is showen in Sec. III C and compared with the exact results obtained in the thermodynamic limit. Finally, Sec. IV gives a summary.

A. The CheMPS method
To carry out a CheMPS calculation 24 , one needs to first calculate the wave function of the ground state |0 and the corresponding energy E 0 using the DMRG or variational MPS methods. In order to determine the spectral function G AB (ω) defined in Eq. (1) in a frequency interval [0, W ], we rescale and shift the frequency and the Hamiltonian such that the nonzero parts of spectral function are completely concentrated on the frequency interval [−W , W ], where W is usually taken to be slightly smaller than 1 for numerical stability 24 . Then the spectral function becomes where a = W/(2W ), and The δ-function in Eq. (2) can be expanded by an infinite Chebyshev polynomial series 26 . In practical calculation, this expansion is approximated by the first N terms of the Chebyshev polynomials 24 , i.e.
where b = π √ 1 − ω 2 and g n 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 26 , here we use the Jackson damping defined by In Eq. (4), T n (x) are the Chebyshev polynomials defined by the recurrence relation with T 0 (x) = 1 and T 1 (x) = x. T n (x) can be also represented as Substituting Eq. (4) into Eq. (2), we obtain where µ n are the Chebyshev moments defined by To obtain the Chebyshev moments µ n , we first calculate the Chebyshev vectors From the recurrence relation Eq. (6), we find that |ψ 0 = B|0 , |ψ 1 = H |ψ 0 , and These Chebyshev vectors |ψ n span a Krylov basis space.
In the CheMPS calculation, the ground state |0 and all the Chebyshev vectors |ψ n 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 Similarly, other Chebyshev vectors are determined by minimizing the cost function From these Chebyshev vectors, the Chebyshev moments are determined by the formula Given a Chebyshev expansion number N , it is known that the CheMPS has a spectral resolution of order O(W/N ) in the interval ω ∈ [0, W ] 24 . When the bandwidth is broadened or the spectrum changes drastically, one has to increase N to increase the resolution.

B. 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 N . 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: 1. Reorthonormalize the Chebyshev vectors |ψ α using the Gram-Schmidt orthogonalization process so that Here |Ψ α with |Ψ 0 = B |0 are the new orthonormalized basis states. c α are the normalization constants.
2. Represent the Hamiltonian and other operators in the basis space spanned by 3. Redo the Chebyshev expansion and generate N r new Chebyshev vectors |ψ α using H in this new basis space.
In this new scheme, the Chebyshev expansions are applied twice. In the first expansion, N 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 H , which is an N ×N matrix defined in the space spanned by {|Ψ α , (α = 0 · · · N − 1)}, and N r Chebyshev vectors that satisfy the standard recurrence relations of Chebyshev polynominals are generated using H . Here N r is the number of Chebyshev vectors used for approximating the δ-function using Eq. (4) for the Hamiltonian H . It can be set either equal to N or as an independent control parameter. The Chebyshev expansion Eq. (4) becomes exact in the limit N r → ∞.
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 H has a finite dimension. We can diagonalize it to obtain its spectrum where λ n and |n are the eigen-pairs of H , which are arranged in an ascending order, i.e. λ n < λ m if n < m.
In this truncated Hilbert space, the spectral function is given by In the second Chebyshev expansion, the above deltafunctions are blurred by N r 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 N r equal or not equal to N , depending on the resolution required. Eq. (20) is recovered in the limit N r → ∞.
Furthermore, as the spectrum of H is known, we can always set the bandwidth of H , represented by W * , exactly equal to the difference between the largest and smallest eigenvalues of H , i.e. W * = λ N − λ 1 , and find all the Chebyshev vectors without taking any approximation in the second Chebyshev expansion. Clearly, W * < W , thus if the same number of Chebyshev vectors is used as in CheMPS, i.e. setting N r = N , the energy resolution, which is qualitatively determined by the ratio W * /N r , should be higher in ReCheMPS.

C. Finite-size scaling
The CheMPS offers a convenient way to control the accuracy and resolution of the spectral function by simply adjusting the expansion order N . This is an advantage that can be used to analyze the finite-size effects 24 . For example, one can mimic the infinite lattice limit (L → ∞) by choosing a small enough N or more precisely a large W/N to smear out the finite-size subpeaks. However, in practical calculation, it is difficult to adjust N to reduce the finite size effect for the systems with different lattice lengths.
An alternative approach is to take a large N 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 24 . 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 N to do the Chebyshev expansion when the lattice size L or the bond dimension D becomes large.
Here we propose a novel approach to perform the finitesize scaling. In Sec. II B, we show that by diagonalizing the Hamiltonian H in the truncated basis space, we can approximately represent the spectral function as where λ n and |n are the eigen-pairs of H , defined in Eq. (19), and is the spectral weight corresponding to λ n .
Let us denote a n as the midpoint between λ n−1 and λ n , i.e. a n = λ n−1 + λ n 2 , The average spectral weight within the interval between a n and a n+l (l > 0) is given by where ∆ n,l = a n+l − a n is the energy difference between a n+l and a n , and ω n,l = a n + a n+l 2 is the average energy of these two points. Here l should be taken such that the interval between a n and a n+l 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 S = 1/2 XY and the S = 1/2 antiferromagnetic Heisenberg models in one dimension. These models are respectively defined by the Hamiltonians and where S r = (S x r , S y r , S z r ) is the spin operator and L is the lattice length. For both cases, open boundary conditions are assumed.
We have calculated the spin-spin correlation function defined by where S z k is the spin operator defined in momentum space and k = nπ/(L+1) (n = 1, ..., L) is the quasi-momentum.
In the calculation of this correlation function, the energy spectrum is bounded by the band width W , which equals 2 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 W 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 W = 10 for both models. By doing this, we do not need to employ the energy truncation method used in Ref. [24].
During the second Chebyshev expansion, as already mentioned in Sec. II B, we take the bandwidth exactly equal to the difference between the largest and the smallest eigen-energies, i.e. W * , of the approximate Hamiltonian H . For all the calculation we have done, we find that W * < 10 and varies slightly with N . The results presented below are obtained by taking N r = N if not particularly specified.

A. The XY-model
The 1D XY model can be mapped onto a free spinless fermion model through the Jordan-Wigner transformation 27 . The spin-spin correlation function S(k, ω) contributes mainly by the particle-hole excitations of spinless fermions, whose energy is bounded respectively from below and above by As k → π, the spin-spin correlation function can be exactly represented as where Ω k = 2|cosk|. In the limit L → ∞, it becomes Figure 1(a) shows S(k = π, ω) as a function of frequency ω for the 1D XY model obtained with ReCheMPS on the L = 24 lattice. As expected, the spectral resolution is improved with the increase of N . The peak positions of S(π, ω), as illustrated in Fig. 1(b), converge rapidly to the exact energies, Ω k = 2| cos k|, again with the increase of N . 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 S(π, ω) changes with D for the system of L = 50. By comparison with the result shown in Fig. 2, we find that a larger D 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 D = 128.
With the further increase of the lattice size, the computational cost (including both computation time and memory) increases rapidly with increasing N , and it is difficult to use a large N to do the Chebyshev expansion. This limits the resolution of the CheMPS calculation. However, in the ReCheMPS calculation, we can use a large N r to increase the resolution. Fig. 4 shows, as an example, how the resolution changes with the increase of N r in ReCheMPS for the XY model with L = 100 and N = 80. When N r = N = 80, the resolutions of both ReCheMPS and CheMPS are very poor. However, when we increase N r to 500, the resolution is greatly improved and the ReCheMPS produces all the low energy peaks accurately.

B. The Heisenberg model
For the spin-1/2 Heisenberg model, the asymptotic behavior of the dynamical spin structure factor is known from the Bethe ansatz solution 28 . It diverges at k = π in the limit ω → 0 and L → ∞ 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 D in a large lattice system. Figure 5(a) compares the results obtained using CheMPS and ReCheMPS with D = 64 for the Heisenberg model of L = 100. 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 low-energy 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 D. The accuracy in the high-energy spectral function, as shown in Fig. 5(b),can be improved by increasing D 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 log-log plot of S(π, ω) 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 S(k, ω) for the 1D Heisenberg model in the whole momentum space. An intensity plot of S(k, ω) is shown in Fig. 6. S(k, ω) 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. ω min Heis = π 2 | sin k|, ω max Heis = π| sin

C. Finite-size scaling
Here we use the method introduced in Sec. II C to perform a finite size scaling analysis for the spectra obtained using the effective Hamiltonian H . In practical calculation of the spectral function using Eq. (24), the energy interval ∆ n,l 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 l = 1. But in the following cases, one should enlarge the value of l so that ∆ n,l 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 ∆ n,l does not change significantly the total spectral weight within this enlarged energy interval. Practically, one can start with l = 1, 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 l to contain all those energy levels, provided that ∆ n,l 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 spin-spin 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 L → ∞. 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 high-resolution 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 eigen-energies 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 eigen-energy. This can reduce the finite size effect and provide a reliable and accurate approach to scale the finite-lattice 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.
Our method works not just in one dimension, it can be generalized to higher dimensions, either using MPS or using other kinds of tensor network states, such as projected entangled pair states (PEPS) 3 or projected entangled simplex states(PESS) 4 .