Spectral functions and time evolution from the Chebyshev recursion

We link linear prediction of Chebyshev and Fourier expansions to analytic continuation. We push the resolution in the Chebyshev-based computation of $T=0$ many-body spectral functions to a much higher precision by deriving a modified Chebyshev series expansion that allows to reduce the expansion order by a factor $\sim \frac{1}{6}$. We show that in a certain limit the Chebyshev technique becomes equivalent to computing spectral functions via time evolution and subsequent Fourier transform. This introduces a novel recursive time-evolution algorithm that instead of the group operator $e^{-iHt}$ only involves the action of the generator $H$. For quantum impurity problems, we introduce an adapted discretization scheme for the bath spectral function. We discuss the relevance of these results for matrix product state (MPS) based DMRG-type algorithms, and their use within the dynamical mean-field theory (DMFT). We present strong evidence that the Chebyshev recursion extracts less spectral information from $H$ than time evolution algorithms when fixing a given amount of created entanglement.

Figure 1

(Top) Typical example of a discontinuous spectral function due to the restriction to a given symmetry sector. In this case, this is the particle contribution of the spectral function of the single-impurity Anderson model (SIAM) with semielliptic bath density of states of half-bandwidth $2v$ and interaction $U/v=4$ [16], taken from Ref. [14]. The spectral function is given by (8) using $|{\psi }_0\rangle =c^{†}|E_0\rangle$, where $|E_0\rangle$ denotes the half-filled ground state and $E_{\text{ref}}$ is the Fermi energy $E_0$. Here, only the shape of the scalar function is of importance, therefore we postpone the model definition to (30). The same spectral function is obtained for the local density of states of the first site for spinless fermions hopping on a semi-infinite chain with tunneling $v$ and an interaction of $U/v=4$ that acts only on the first site. (Bottom) Comparison of convergence of the Chebyshev moments of $A^{>}\left(\omega \right)$ with its redefinitions ${\tilde{A}}^{>}\left(\omega \right)$ and $A\left(\omega \right)$, giving rise to moments ${\mu }_n^{>},{\tilde{\mu }}_n^{>}$, and ${\mu }_n$, respectively. The full spectral function (18) for this example is $A\left(\omega \right)=A^{>}\left(\omega \right)+A^{>}(-\omega )$ as here, $A^{<}\left(\omega \right)=A^{>}\left(\omega \right)$ due to particle-hole symmetry. All of this is for the setup $b=0$ using a rescaling of $a=100v$ in Eq. (10).

Figure 2

(Top) Integrand for computation of moments for the spectral function shown in Fig. 1 in the two setups $b=0$ (left) and $b=-0.995$ (right). (Bottom) Comparison of convergence of moments computed with the integrands shown in the top panels. In all of that, $a=100v$.

Figure 3

(Top) Input spectral function and reconstructed spectral functions using linear prediction for $N=200$ computed Chebyshev moments. We compare our proposal with the original proposal [13], where for the $b=0$ setup the expansion of the full spectral function was extrapolated. (Left) In the first case, we use the Chebyshev expansion of ${\tilde{A}}^{>}\left(\omega \right)$ in the $b=-0.995$ setup (red dashed lines), and in the second, we use the Chebyshev expansion of the full spectral function $A\left(\omega \right)$ in the $b=0$ setup (solid green line). (Right) Error of these functions $|A_{\text{reconst}}\left(\omega \right)-A_{\text{input}}\left(\omega \right)|$ for both setups. (Bottom) The error ${max}_{\omega \ge 0}|A_{\text{reconst}}\left(\omega \right)-A_{\text{input}}\left(\omega \right)|$ vs number $N$ of computed Chebyshev moments. Here, we also show results for using the Chebyshev expansion of ${\tilde{A}}^{>}\left(\omega \right)$ in the $b=0$ setup (light green crossed line). This is different from using the full spectral function $A\left(\omega \right)$ in the $b=0$ setup. Lower error levels only occur for much higher expansion orders than shown in the panel.

Figure 4

Reconstruction of the spectral function of Fig. 1 represented by the discrete Hamiltonian (29). (Top left) Input function and reconstructed functions using $L=80,N=200,a=100v,b=-0.995$. We compare the linear discretization (26) with the cosine (28) discretization. (Top right) Difference of input and reconstructed functions of the top left panel. (Bottom left) The error ${max}_{\omega \ge 0}|A_{\text{reconst}}\left(\omega \right)-A_{\text{input}}\left(\omega \right)|$ vs number $N$ of computed Chebyshev moments using a cosine discretization. (Bottom right) Error for linear discretization.

Figure 5

Comparison of MPS computed spectral functions in the two setups studied in the previous sections, with data by Raas et al. [16]. Solid lines refer to the method that uses $\tilde{A}\left(\omega \right)$, dashed lines refer to the method that uses the full spectral function $A\left(\omega \right)$. (Top left) For $L=80$ and two exemplary expansion orders. (Top right) Errors of comparison in top left panel. (Bottom left) Plot of the maximum error vs expansion order $N$. (Bottom right) Plot of the maximum error vs computer time.

Figure 6

(Left) Time evolution of a particle created on the first site of a chain of length $L=100$ with hopping $v=1$, that is, $|{\psi }_0\rangle =c_0^{†}|\text{vac}\rangle$. We compare the time evolution of the Green's function $i{G}^{>}\left(t\right)=e^{i{E}_{0}t}\langle {\psi }_0|\psi \left(t\right)\rangle$ (shown as blue crosses) with the Chebyshev moments ${\mu }_n=\langle {\psi }_0|{\psi }_n\rangle$ (shown as green dots) obtained when computing the recursion $|{\psi }_n\rangle =2H/a|{\psi }_{n-1}\rangle -|{\psi }_n-2\rangle$. Hopping amplitudes $v_l$ are obtained from the discretization of the spectral function of Fig. 1. A qualitatively equivalent behavior is obtained for a chain with homogeneous hopping $v_l=v$. (Right) Long-time evolution of the Green's function (blue crosses), and difference of the Green's function and the Chebyshev moments (red dots). The horizontal dashed line marks the prefactor of the error estimate (37), which is computed as ${\sigma }^3/a^2=8\times {10}^{-4}$.

Figure 7

Time evolution of the single-particle excitation $c_{0\uparrow }^{†}|E_0\rangle$ in the half-filled single-impurity Anderson model (30) with semielliptic density of states of half-bandwidth $2v$ and interaction $U/v=4$ for $L=40$. Computations using an MPS Krylov algorithm with error tolerance ${\varepsilon }_{\text{kry}}=5\times {10}^{-4}$ and the Chebyshev recursion (40) for different error tolerances ${\varepsilon }_{\text{che}}=5\times {10}^{-4}$ and ${\varepsilon }_{\text{che}}=5\times {10}^{-5}$. (Top left) Time evolution of Greens's function. Both algorithms produce the same result upon using the smaller error tolerance for the Chebyshev algorithm. The legend is found in the top right panel. (Top right) Maximal bond dimension, located at the central bond. The Krylov time evolution leads to a much smaller maximal bond dimension, as its computation produces a faithful result already with the relatively high error tolerance of ${\varepsilon }_{\text{kry}}=5\times {10}^{-4}$, for which the Chebyshev algorithm shows strong errors in the Green's function. (Bottom left) Difference of overlap of Chebyshev and Krylov evolved states, comparing the ${\varepsilon }_{\text{che}}=5\times {10}^{-5}$ with the ${\varepsilon }_{\text{kry}}=5\times {10}^{-4}$ computation. The difference of overlap is bounded by the analytical prediction of (37), except for few exceptions that lie above it. These exceptions are of purely numerical origin as they are not visible in any other quantity. For the highest times shown, truncation errors have accumulated so much that the analytical prediction starts to fail. (Bottom right) Bond entanglement entropy at center bond.

Figure 8

Time evolution of a one-dimensional fermionic Hubbard model on $L=90$ sites, with an interaction of $U/v=4$ and nearest-neighbor hopping $v$ starting from a product state (double occupation in the center of the system). (Left) Chebyshev computation using ${\varepsilon }_{\text{che}}=0.0001$. (Right) Krylov computation using ${\varepsilon }_{\text{kry}}={\varepsilon }_{\text{che}}$. (Bottom) Detailed comparison for the occupation of specific sites. Chebyshev results are shown as dashed lines, Krylov results are shown as solid lines. Deviations are smaller than 1%.

Figure 9

(Top) Test function consisting of two Lorentzians (B1). (Center) Corresponding Chebyshev moments in the three different setups ${\mu }^{>},{\tilde{\mu }}^{>}$, and $\mu$, analogously to the bottom panel of Fig. 2. (Bottom) Error of reconstructed spectral function, analogously to the bottom panel of Fig. 3. All of this is for $a=100$.

Figure 10

(Top) Test function consisting of two Gaussians (B1). (Center) Corresponding Chebyshev moments in the three different setups ${\mu }^{>},{\tilde{\mu }}^{>}$, and $\mu$, analogously to the bottom panel of Fig. 2. (Bottom) Error of reconstructed spectral function, analogously to the bottom panel of Fig. 3. All of this is for $a=100$.

Figure 11

(Top) Particle spectral function of half-filled noninteracting SIAM. (Center) Corresponding Chebyshev moments in the three different setups ${\mu }^{>},{\tilde{\mu }}^{>}$, and $\mu$, analogously to the bottom panel of Fig. 2. (Bottom) Error of reconstructed spectral function, analogously to the bottom panel of Fig. 3.