{'fraction_non_alphanumeric': 0.0334865397242285, 'fraction_numerical': 0.010735390676296783, 'mean_word_length': 4.496391194514616, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spectral Functions with the Density Matrix Renormalization Group: Krylov-space Approach for Correction Vectors

Frequency-dependent correlations, such as the spectral function and the dynamical structure factor, help understand condensed matter experiments. Within the density matrix renormalization group (DMRG) framework, an accurate method for calculating spectral functions directly in frequency is the correction-vector method. The correction-vector can be computed by solving a linear equation or by minimizing a functional. This paper proposes an alternative to calculate the correction vector: to use the Krylov-space approach. This paper then studies the accuracy and performance of the Krylov-space approach, when applied to the Heisenberg, the t-J, and the Hubbard models. The cases studied indicate that Krylov-space approach can be more accurate and efficient than conjugate gradient, and that the error of the former integrates best when a Krylov-space decomposition is also used for ground state DMRG.


I. INTRODUCTION
In the last two decades, several approaches have been proposed for the calculation of frequency dependent correlations or spectral functions of low dimensional strongly correlated systems within the density matrix renormalization group (DMRG) method. 1,2 In 1995, Hallberg proposed a Lanczos-vector method based on a continued fraction expansion. 3 Later, Kühner and White 4 showed that this method is suitable only for spectra consisting of few discrete peaks. The same paper introduced the correction vector method to calculate spectral functions. Jeckelmann later developed a variational improvement to the correction-vector method-the dynamical DMRG. 5 It has since been used successfully in several studies 6 to calculate the dynamical properties, directly in frequency. The main disadvantage of the correction-vector method is two-fold: (i) dynamical properties need to be computed in small intervals of frequency; (ii) the artificial broadening, η, which always needs to be introduced at some point in the calculation, is set from the start, and to change it a new run needs to be carried out. The η broadening can be viewed as convolution of the exact spectral function with a Lorentzian of the same width. Because one is interested in the exact spectrum, one needs to perform a deconvolution of the spectrum-an ill-defined operation. Many works have presented different and successful deconvolution procedures for correction-vector DMRG spectra of one dimensional 7,8 and quantum impurity problems. [9][10][11][12] Recently, ref. 13 has proposed a blind deconvolution algorithm, having the advantage of reproducing quite well sharp singularities such as power-law band edges and excitonic peaks, but the drawback of introducing artificial shoulder-like spectral features.
Within the DMRG framework, another method for calculating the dynamical properties consists of first computing the correlation functions using time-dependent DMRG, [14][15][16] and then space-time Fourier transform-ing them into real frequency. Recently, ref. 17 presented a highly efficient and accurate adaptive method using Chebyshev polynomials in combination with matrix product states (MPS). In order to compute timedependent correlation functions over a long time interval, ref. 18 has shown that one can apply an efficient prediction method. In addition, ref. 19 has shown that linear prediction can be used as a method to extrapolate Chebyshev or Fourier expansions of spectral functions. Within the MPS formulation of DMRG, other powerful methods have been developed. 17,20,21 Ref. 20 and 21 have proposed an improvement of the continued fraction expansion method of Hallberg in the MPS language.
In this paper, we focus on the correction-vector DMRG method, and propose an alternative method to calculate the correction vector. In its traditional formulation, the imaginary part of the correction-vector, being directly proportional to the spectral function, is the solution of an inhomogeneous set of linear equations. This set of equations can be solved explicitly with the conjugategradient method. The conjugate-gradient method's convergence error is independent of the ground state DMRG error. Dynamical DMRG 5 recasts the calculation of the correction-vector as an elegant variational principle, defining a functional that one then minimizes. The main advantage of dynamical DMRG is the possibility of finding a solution with an error smaller than that of the traditional correction-vector method.
We here propose an alternative approach: to calculate the correction vector using a Krylov-space 22 decomposition instead of solving a linear equation. Krylov space decomposition approaches have found application in several fields of science and engineering, where the computation of functions of large sparse matrices is often the main problem to be solved. In physics, they have been used in several contexts 23 as a technique to solve partial or ordinary differential equations, such as numerical general relativity and electrodynamics. 24 The correction vector is not then calculated with a separate algorithm, be it conjugate gradient as originally proposed, or minimization of an auxiliary functional as in dynamical DMRG, but with the following algorithm instead. At each DMRG step (i) tri-diagonalize and successively diagonalize the Hamiltonian of the problem in the current basis; (ii) calculate the correction-vector in the diagonal basis and rotate it back to the current basis. The error profile is then the same as standard DMRG if Krylov-space decomposition is used for the ground state computation, because Krylov-space decomposition is needed at each step of the ground-state DMRG algorithm. Our proposal therefore integrates better with ground-state DMRG. Translated into MPS language, our method works on the same space spanned by the effective Hamiltonian where all but one site of the tensor network (MPO-MPS product) has been contracted, where MPO stands for matrix product operator. Moreover, to avoid Lanczos ghost states, 25,26 we have been careful not to overconverge the Lanczos iterative process by not setting the Lanczos error to be extremely small. This paper is organized as follows. Sec. II presents the Krylov-space approach for the calculation of the correction-vector. Sec. III applies our equations to the Heisenberg spin chain, to the t-J model, and to the Hubbard model. We have calculated the S(k, ω) for the Heisenberg and Hubbard models, while the A(k, ω) for a t-J chain. Section III D shows that Krylov approach can efficiently calculate the dynamical spin structure factor of the t-J model on a ladder geometry. Section IV compares the frequency resolution and the computational performance of the conjugate gradient method with the Krylov-space approach for the simulated models. The last section presents a summary and conclusions.

II. CONJUGATE GRADIENT AND KRYLOV METHODS
This section briefly recalls the basics of the correctionvector method. One is interested in the calculation of the Green's function where |ψ 0 is the ground state of some HamiltonianĤ with ground-state energy E 0 ,Â andB are operators associated with the dynamical correlation function to be calculated, and z ≡ ω+iη, where ω is the real frequency and η is a positive constant that provides a finite broadening of the Green's function peaks. This procedure is made rigorous by realizing that what we call Green's functions are actually distributions. The correction-vector associated with Eq. (1) is defined by where the vector |A ≡Â|ψ 0 is assumed to be real. Within the traditional formulation of DMRG, 4 the correction-vector method uses a multi-target approach: at each step of the DMRG algorithm, one targets the ground state of the system |ψ 0 , the vector |A and the |x(z) in the reduced density matrix, for each frequency value ω and broadening η. After one obtains the correction vector, one can calculate the Green's function G(z) using A direct approach to calculate the correction-vector is to solve the following set of linear equations in the local DMRG basis: Eq. (4) could be solved with the generalized minimal residual method (GMRES), which uses the Arnoldi algorithm to find a generalized decomposition of the matrix M ≡ z + E 0 −Ĥ onto a Krylov subspace of much smaller dimension than the local Hilbert space. 27 But such approach would have the drawback that the matrix M ≡ z + E 0 −Ĥ is not Hermitian due to the presence of the factor η > 0 in z = ω + iη. The convergence error of the algorithm is given by the condition number of the matrix M , κ(M ). 27 The larger the condition number, the greater is the number of iterations needed for solving the set of equations, and the smaller is the improvement of the solution at each iteration step.
To avoid dealing with non Hermitian matrices, the imaginary part of the correction-vector, |x(z) Im , is calculated by solving the following set of linear equations instead of Eq. (4). The real part of the correction-vector |x(z) Re is then calculated from its imaginary part using In this case, the matrix M ≡ [(E 0 + ω −Ĥ) 2 + η 2 ] is real, symmetric, and positive definite, therefore at each DMRG step, Eq. (5) can be solved with the conjugate gradient method. 27 This method can be thought of as a particular case of GMRES for real symmetric and positive definite matrices. 27 The main problem with this second approach is that the condition number of the matrix M , κ(M ) in Eq. (5) is larger: roughly the square of the condition number of Eq. (4). Recall that, for a system of inhomogeneous set of equations M x = b for the unknown vector x, the error in the conjugate-gradient algorithm at step k is given by the norm of the residual r k = b − M x k , where x k is the approximate solution at that step 27 . In the next sections, we assume that a measure of the error for the conjugategradient algorithm is given by the number of iterations required for reaching a certain tolerance on the solution.
In this paper, we propose an alternative for calculating the correction-vector in Eq. (2). At each step of the DMRG, we perform a Lanczos tri-diagonalization of the Hamiltonian and a successive diagonalization in the current DMRG basis; the vector |x(z) is then calculated directly as where D is the diagonal form of the Hamiltonian operator H in the current basis. We have assumed that is a faithful representation of the Hamiltonian when applied to the starting vector. V represents the matrix of the Lanczos vectors spanning the Krylov space, and T the representation of the Hamiltonian in tridiagonal form. The tri-diagonalization is then followed by a small full diagonalization of T , yielding the matrix of eigenvectors S. The first equality in Eq. (8) holds with the Lanczos error, unlike in the conjugate gradient method, where the error is separate from the DMRG. Therefore, for the Krylov method the error coincides with the Lanczos error which occurs in the tri-diagonalization of the Hamil-tonianĤ. The accuracy of the approximation Eq. (7) can be estimated from the high frequency expansion of the Green's function Eq. (1) in the limit of η → 0 28 Using approximation Eq. (8) in the above equationusually satisfied after n Lanczos iterations for a certain required accuracy-the approximated Green's function reproduces the first n coefficients of the expansion within Lanczos error because the vector belongs to the Krylov space for all r ≤ n − 1. We have here neglected the error in the determination of the ground state |ψ 0 . The expansion Eq. (10) therefore produces a good approximation for the first n moments of the spectral func- 28 As will be confirmed numerically in the next sections, the accuracy and efficiency of this approach is then best at low frequencies, and worsens for high excitation energies. In our later analysis, we shall find that the number of iterations required to reach a certain accuracy in our Krylov approach is of the order of κ(M ), showing evidence that that the conjugate gradient method is less efficient. When the number of Lanczos iterations is sufficiently small, then the dimension of the matrix T is of the order of a few hundreds; for all the models and geometries investigated in this paper, we find that about 150 iterations n are at most needed in our simulations. The computational cost of the diagonalization S is then negligible. On the contrary, the number of conjugate gradient iterations necessary to obtain the same accuracy is much larger.
Our method is similar to DMRG continued-fractionexpansion or CFE in that it uses a tridiagonal decomposition of the Hamiltonian to compute the correction vector. The approach is therefore equivalent to fully diagonalizing the tridiagonal Hamiltonian and obtaining the matrix S on the one hand, or doing a CFE on the other. Our method is dissimilar to DMRG CFE in that it is a correction-vector method, and thus Lanczos vectors are not targeted. In fact, DMRG CFE refers to the method where all or some of the Lanczos vectors are targeted in the reduced density matrix at each DMRG step. In the original paper 3 where S i denotes the spin operator at site i. We choose J = 1 as unit of energy for this model. As in the original paper of Kühner and White, 4 the longitudinal dynamical spin structure factor is calculated with correction-vector DMRG as where the operator S z c is applied to the ground state at the center of the lattice for calculating the correction vector |x(ω + iη) . Then, the correlator is computed for all the sites j of the lattice, and for a fixed value of frequency ω and broadening η. The above quantity is finally transformed to momentum space as where the quasi-momenta k = πn L+1 with n = 1, .., L are appropriate for open boundary conditions. The DMRG implementation used throughout this paper is discussed in the supplemental material, which can be found at https://drive.google.com/open?id= 0B4WrP8cGc5JHX3h6M25zUVhiQmc.
It is known from the Bethe ansatz solution 29-32 of Hamiltonian Eq. (11) that the upper and lower boundary of the spin excitation manifold of the infinite system, called the des Cloiseaux-Pearson (dCP) dispersion relations, are given by where q represents the momentum and not the quasimomentum. Moreover, it is known that the S(k, ω) diverges as as ω approaches the lower boundary ω l from above. This divergence has its profound origin in the Luttinger liquid nature of the ground state, and describes the instability of the model toward antiferromagnetic ordering. Numerically, because one has always finite size systems, one usually cuts off the divergences at ω − ω l 1/L, so that one has peaks of finite height Fig. 1 shows the dynamical spin structure factor of an antiferromagetic (J = 1) Heisenberg spin chain of size L = 64. The Krylov-space correction-vector method described in the previous section has been applied. We use m min = 64 and a maximum of m = 1000 DMRG states by keeping the truncation error no bigger than 10 −8 . Moreover, we set the maximum number of Lanczos iterations necessary to calculate the correction-vector to 1000, keeping the tri-diagonalization error no bigger than 10 −7 . For all the models and geometries investigated in this paper, we have found that the number of Lanczos iterations needed is at most about 150. Even at a finite system size, the dCP relations, which are exact in thermodynamic limit, are very well reproduced by the DMRG data. As explained thoroughly in ref. 5, finite size scaling of dynamical correlation functions with correction-vector DMRG should be done carefully: the artificial broadening η of the spectral function peaks should be rescaled as function of system size, such that lim L→∞ η(L) = 0. We find that the scaling of the peaks maxima in Eq. (17) is well reproduced by imposing η = c/L, where c is a constant of the order of the width of the full spectrum. Panels (a) and (b) of fig. 2 show the cuts of the S(k, ω) spectrum at k = π and k = π/2 as a function of ω for different system sizes. When using open boundary conditions, we use the convention that k = π refers to k = πL/(L + 1), while k = π/2 refers to k = πL/(2(L + 1)). We have is approximately constant as a function of the system size and approximately close to the expected thermodynamic limit value ω l (π/2) π/2. Panels (c) and (d) show the maxima of the peaks obtained in panel (a) and (b) as a function of the system size. With a linear fit, we have verified that the peaks maxima have the correct scaling as described by the relations in Eq. (17). Fig. 3, panels (a) and (c) show the dynamical structure factor for a Heisenberg chain with L = 48 sites, calculated with the Krylov-space and the conjugate gradient method, respectively. In this figure, η = 0.1 and m = 800 DMRG states were kept in our simulations. The frequency resolution of the Krylov-space approach is much better than that provided by the conjugate gradient method. While in the Krylov-space approach the spectral weight outside the region defined by the dCP relations is practically zero, the conjugate gradient method gives a spectral weight spread everywhere in the frequency interval investigated, with much less defined peak features. We have tried to use the same parameters for both approaches, but could not do so precisely, as we now explain.
We have set the maximum number of conjugate gradient iterations to 1000 trying to keep the error no bigger than 10 −7 . Unfortunately, we have found that, for most of the frequency interval investigated, a number of iterations much larger than 1000 is necessary to reach convergence. The conjugate-gradient error can thus be as high as 10 −1 . This explains why the spectrum shown in panel (c) has much less frequency resolution than the spectrum calculated with the Krylov-space approach. Panel (e) shows a more detailed comparison between the spectra calculated with the two approaches. Cuts at k = π and k = π/2 are examined, indicating that conjugategradient is able to capture only qualitatively the main features of the spectra. For the cut at k = π the position of the peak at low frequency ω 0 is correctly reproduced, but the spectral weight is smaller than that provided by the Krylov-space approach by a factor of 3. A different behavior is observed for the peak at k = π/2. Here, the position of the spectral peak is shifted to lower frequency, while the spectral weights are also redistributed to lower frequencies; see dashed (red) line in panel (e).

B. Dynamical spin structure factor of Hubbard chains
In the same figure 3, panels (b) and (d) show the dynamical spin structure factor calculated with the two methods outlined above for a Hubbard chain of L = 48 sites at half-filling and at U/t = 4. The Hubbard model Hamiltonian is given bŷ where we have used standard notation. In these panels, t = 1 is assumed as unit of energy. Notice the similar-  (a-b), and the conjugate gradient method (c-d) for a Heisenberg (panels (a-c)) and a Hubbard (panels (b-d)) chain of L = 48 sites, using η = 0.1 and m = 800 DMRG states (mmin = 64). Panel (e-f) shows cuts at k = π and k = π/2 of the spectra in (a-c) and (b-d) as a function of frequency. In the conjugate gradient method, a maximum of 1000 iteration steps has been imposed, regardless of the error.
ity between the Hubbard and the Heisenberg chain spectrum in panels (a) and (b) calculated with Krylov-space approach. The Bethe ansatz solution 33 at half-filling and in the limit of large Coulomb repulsion U t, confirms that the Hubbard model on a chain is gapless in the spin sector. In this case, the Hubbard model behaves as the Heisenberg model. Similar to the Heisenberg case, the quality of the spectrum obtained with the conjugate gradient is much worse than that obtained with Krylov. Here again, the spectral weight is distributed everywhere outside the region delimited by the dCP relations. A detailed comparison of cuts of the spectrum at k = π and k = π/2 in panel (f) show again that the position of the spectral peaks is only qualitatively captured by the conjugate gradient.

C. Spectral function of a t-J chain
where we have used standard notation. For this model, we choose the unit of energy t = 1, and study the case where J = 0.5. The spectral function A(k, ω) = ξ=± A ξ (k, ω) consists of two branches: the photoemis- sion ξ = −, and the anti-photoemission ξ = +.
x (construction operator) at site x. The two branches of the spectrum are calculated with correctionvector DMRG and Fourier transformed to momentum space, as described in the previous subsection for the dynamical spin structure factor S(k, ω).
Panels (a) and (b) of Fig. 4 show the spectral function of the t-J chain calculated with the Krylov-space and conjugate gradient approaches with η = 0.1 and keeping a small number m min = m = 100 of DMRG states. The spectral function of the t-J model contains the phenomenon of spin-charge separation. Below the Fermi level (indicated by a solid black line), the spectral weight is concentrated on the spinon and holon bands. The spinon band forms an arch of amplitude J 0.5 connecting the two Fermi points (−k F , ω − µ = 0) and (k F , ω − µ = 0). Because the filling is N/L = 2/3, then k F = π/3. The holon bands depart from the spinon band at the Fermi points in two approximately straight lines with slope 2t/k F . As expected for a Luttinger liquid, the shadow bands extend in frequency well beyond |k F |.
Panel (c) and (d) of Fig. 4 show the spectral function with a larger number of DMRG states, m = 400. Regardless of the method used, the frequency resolution of the spectrum is improved. The frequency resolution of the holon and spinon branches of the spectrum is worse in the conjugate gradient case. Nevertheless, the qualitative features of the spectrum are captured by the conjugate gradient method too.

D. Dynamical spin structure factor of a t-J ladder
We now show that our Krylov based approach can be used to calculate efficiently the magnetic excitation spectrum of systems with more complex geometries than chains. As a case study, we analyze a t-J model on a ladder geometry, with Hamiltonian Recently, the ground state properties and the spectral function A(k, ω) in the limit of one hole doping have been studied with DMRG on large system sizes. 34,35 Yet it is well known that finite dopings are more difficult to treat with DMRG, therefore our calculation required a comparable effort. In the past, this model was thoroughly studied in the context of cuprates. [36][37][38] In the undoped limit, it has been well established that the t-J model has a spin gap due to the particular ladder geometry, and, as in the chains' case, the physics can be described in terms of the Heisenberg ladder. Upon doping, superconductivity mainly occurs in the d-channel, as described in refs. 36,39-42. As suggested in ref. 43, neutron scattering data could provide important evidence for a pairing mechanism based on the exchange interaction J. The physics of t-J (and Hubbard) two leg ladders has been studied with many techniques ranging from bosonization to DMRG to exact diagonalization. The mostly studied case has been the isotropic regime, where t x = t y = t. The possibility of pairing has been also investigated in the anisotropic limit, where J y /J x = 1 and t y /t x = 1. At half-filling a Heisenberg description applies, and the spin gap has been shown to persist until very small values of the J y parameter. 36 Away from half-filling, the authors of ref. 44 studied the spin gap and the superconducting binding energy of the hole pairs showing that they can be maximized by tuning the anisotropic ratios to t y /t x 1.25, and J y /J x 1.56. Fig. 5 shows our results obtained with DMRG for the spin structure factor. In the case of the ladder, the dynamical structure factor in momentum space has two components because the momentum in the y direction has only two possible values: k y = 0 and k y = π.
Let t x = 1 be our unit of energy, and let us study the spectrum for the set of parameters t y = J y = J x = 1. The characteristics of the spectrum are completely different from the case of decoupled chains (not shown). Indeed, a spectral gap appears clearly, both in the k y = 0 and in the k y = π components of the spectrum. In particular, a very dispersive gap arises in the k y = 0 component with amplitude at k = 0 ω gap (0) 0.75. In the same component, the spin excitations form a triangular structure where all the spectral weight is concentrated. Two separate dispersive arcs can be noticed. The S(k, π, ω) component shows a very defined and peculiar shape. A spin gap can be observed, with ω gap 0.5, a value smaller than the gap in the other component of the spectrum. Above that gap, the dispersion curves of the spin excitations describe the low energy boundary of the non-interacting tight-binding chain results. As in the case of a t-J chain, the spectral weight is redistributed to low energy, but the weight decreases rapidly away from the low boundaries.

IV. COMPUTATIONAL PERFORMANCE
This subsection compares the computational performance of the Krylov-space and conjugate gradient methods for the dynamical spectra studied in the previous sections. Panel (a) of fig. 6 shows the CPU time needed at each frequency ω for the two methods, when calculating the dynamical spin structure factor of the Heisenberg model, investigated in section III A. With the Krylovspace approach, the CPU time is smaller in the low frequency regime than in the high frequency one. At larger frequencies, larger CPU times are needed. As outlined in the previous section, we have set a maximum number of Lanczos iterations equal to 1000, and kept the error no larger than 10 −7 . This accuracy is reached in about 50 iterations in the low frequency regime 0 < ω < 0.5, while the number of iterations needed increases up to 150 in the large frequency regime. The conjugate gradient method has a similar CPU time profile; see the circle (red) symbols. As already mentioned in the previous section, we have set the maximum number of conjugate gradient iterations to 1000, and have considered the method converged if the error is smaller than 10 −7 . However, for most of the frequency interval investigated, the error can be as high as 10 −1 with 1000 iterations or less. We have found that increasing the maximum number of iterations needed by the algorithm to reach an error less than 10 −7 increases the CPU time even more that what we have shown in this paper. As seen in panel (a) of fig. 7, the Krylov method is faster than the conjugate gradient by more the one order of magnitude in the very low frequency regime. But it is not just faster, the conjugate gradient has not even converged despite having used a large amount of CPU time and iterations. Panel (b) of fig. 7 compares the CPU time performance of the two methods for the spectrum of the Hubbard model. Here the performance of the conjugate gradient is better, but still at least a factor of 3 smaller than that of the Krylovspace approach. fig. 6 shows the CPU time as function of the frequency for the spectral function A(k, ω) of the t-J model studied in sec. III C. The Krylov-space approach is again substantially faster than the conjugate gradient. Even if slightly visible, the Krylov approach has its best performance at low energy, that is, for ω−µ 0, where it is 20% faster than at higher frequencies. On the contrary, because of the convergence issues, we have not been able to obtain a well defined CPU dependency on frequency for the conjugate gradient method. Panel (c) of fig. 7 shows the comparison between the two methods in the low frequency regime. As with the other models investigated, the Krylov method is faster than the conjugate gradient by a factor of 3. In panels (a) and (b), the grey and red bars indicate the CPU times for ω = 0 and ω = 2, respectively. In panel (c), the grey and red bars indicate the CPU times for ω − µ = 0 and ω − µ = 1, respectively.

V. SUMMARY AND CONCLUSIONS
This paper proposes an alternative method for computing correction-vectors based on a Krylov-space decomposition. We have tested the quality of our approach by studying the dynamical spin structure factor of a Heisenberg and Hubbard chain, and the spectral function of a t-J chain. We have also shown that the method is general, and applicable without restriction or further approximations to both more complex models and geometries. Beyond chains, we have shown in section III D that our method can be successfully applied to the calculation of the dynamical spin structure factor of a t-J model on a ladder. The supplemental material, which can be found at https://drive.google.com/open?id= 0B4WrP8cGc5JHX3h6M25zUVhiQmc, provides a pointer to the full open source code, input decks and additional computational details.
In all cases investigated, we have found that the Krylov-space approach not only provides frequency spectra with much higher frequency resolution, but that it also requires much less CPU time than the the conjugate gradient method. When DMRG ground state itself uses Krylov space, then the calculation of the correction vector integrates better with DMRG. In those implementations that use the Davidson method 45 for ground state, the Davidson method could perhaps be used to calculate the correction vector, as we have used Krylov space here. On the other hand, Chebyshev or Fourierbased methods for the computation of spectral functions with DMRG should be computationally less expensive than the correction-vector method. But correction vector should be better suited for problems where a constant resolution is needed, and where Chebyshev or Fourier-based methods might have limited resolution in parts of the spectrum. For example, in neutron scattering and photoemission spectra-properties much needed to understand quantum magnets, superconductors, and transition metal oxides-resolving fine features of the low-frequency spectrum is of great importance. The correction-vector method can be very precise in the estimation of these small energy gaps, with the only limitation given by the finite broadening η.
Using DMRG to obtain dynamical functions is meaningful because the only other unbiased method, quantum Monte Carlo, does not work directly in real frequency. In the future, we plan to extend our present work on dynamical properties to finite temperature. 46,47