Modified Pade approach for analytic continuation: Application to the zero-gap Kondo lattice model

TL;DR

This paper introduces a modified Pade method for more stable analytic continuation of correlation functions, successfully applied to models including the zero-gap Kondo lattice, revealing insights into Kondo insulators and energy scales.

Contribution

A novel, more robust Pade-based analytic continuation method applied to complex models like the zero-gap Kondo lattice, enhancing stability and accuracy.

Findings

The modified Pade approach is more stable against statistical errors.

Application to the zero-gap Kondo lattice reveals properties of the Kondo insulating state.

Identification of two energy scales related to direct and indirect gaps.

Abstract

A modified Pade approach is presented as analytic continuation for numerical methods that evaluate correlation functions in imaginary time. Instead of the direct analytic continuation of the correlation functions, the Pade method is applied for the self-energy that is then used for deriving the Green's functions at real energies. We find that this modified, self-energy Pade approach is more stable and robust against statistical errors compared to the direct way. The characteristics and success of the modified Pade approach are analyzed by actual calculations for the illustrative cases of the non-interacting Anderson lattice and Hubbard model. A zero-gap Kondo lattice model with linearly vanishing conduction electron density of states at the Fermi level is also studied, where we use the modified Pade approach as analytic continuation. We investigate the properties of the Kondo insulating state including the dependence of the insulating gap on the Kondo coupling and coherence effects. Furthermore, we identify two energy scales from dynamic and thermodynamic quantities, that are associated as a direct and an indirect gap in a band hybridization picture.

Equations

$${\rho }_c(\epsilon +i\delta )=-\frac{1}{\pi }\mathrm{Im}{G}_c(\epsilon +i\delta ),$$

$$G_c(\epsilon +i\delta )=\int d\omega \rho (\omega ){\left(\epsilon +i\delta -\omega +\mu -{\Sigma }_c(\epsilon +i\delta )\right)}^{-1}$$

$${\rho }_c(\epsilon )$$

$$\rho (\omega )=\frac{2}{\pi D^2}\sqrt{D^2-{\omega }^2}$$

$${\mathcal{H}}_{\mathrm{And},\mathrm{lattice}}=\sum _{\mathbf{k},\sigma }{\epsilon }_{\mathbf{k}}{c}_{\mathbf{k},\sigma }^{†}{c}_{\mathbf{k},\sigma }+{\epsilon }_f\sum _{i,\sigma }{f}_{i,\sigma }^{†}{f}_{i,\sigma }+V\sum _{i,\sigma }\left(f_{i,\sigma }^{†}{c}_{i,\sigma }+c_{i,\sigma }^{†}{f}_{i,\sigma }\right),$$

$${\Sigma }_c(i{\epsilon }_n)=\frac{V^2}{i{\epsilon }_n}$$

$${\Sigma }_c(\epsilon )=\frac{V^2}{\epsilon },$$

$${\rho }_c(\epsilon )=-\frac{1}{\pi }\frac{2}{D^2}\mathrm{Im}\left[\epsilon -{\Sigma }_c(\epsilon )-i\sqrt{D^2-(\epsilon -{\Sigma }_c(\epsilon ){)}^2}\right],$$

$${\rho }_c(\epsilon )=\frac{1}{\pi }\frac{2}{D^2}\sqrt{D^2-{\left(\epsilon -\frac{V^2}{\epsilon }\right)}^2},$$

$$a_1=\frac{V^2}{i{\epsilon }_0},a_2=-\frac{i}{{\epsilon }_0},a_i=0\mathrm{for}i>2$$

$${\Sigma }_c(\epsilon )=\frac{a_1}{1+\frac{a_2(\epsilon -i{\epsilon }_0)}{1+\frac{a_3(\epsilon -i{\epsilon }_1)}{1+....}}}=\frac{a_1}{1+a_2(\epsilon -i{\epsilon }_0)}=-i\frac{V^2}{{\epsilon }_0}\frac{1}{[1-\frac{i}{{\epsilon }_0}(\epsilon -i{\epsilon }_0)]}=\frac{V^2}{\epsilon }$$

$$G_c(i{\epsilon }_n)=\frac{2}{D^2}\left[i{\epsilon }_n-\frac{V^2}{i{\epsilon }_n}-i\sqrt{D^2-{\left(i{\epsilon }_n-\frac{V^2}{i{\epsilon }_n}\right)}^2}\right]$$

$${\Sigma }_c^{(\mathrm{error})}(i{\epsilon }_n)$$

$$G_c^{(\mathrm{error})}(i{\epsilon }_n)$$

$${\mathcal{H}}_{\mathrm{Hubbard}}=t\sum _{⟨i,j⟩}{c}_i^{†}{c}_j+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$${\Sigma }_{c,\mathrm{low}-\mathrm{U}}^{(\mathrm{approx})}(\epsilon )=\sigma (\delta ,z_1,y_1,\epsilon )$$

$${\Sigma }_{c,\mathrm{high}-\mathrm{U}}^{(\mathrm{approx})}(\epsilon )=\sigma (0,1,1,\epsilon )$$

$$\sigma (\delta ,z_1,y_1,\epsilon )\equiv \frac{-i\epsilon z_1\Gamma }{\epsilon y_1\Delta -i{\epsilon }^2+i\delta }.$$

$${\Sigma }_{c,\mathrm{low}-\mathrm{U}}^{(\mathrm{approx})}(i{\epsilon }_n)$$

$${\Sigma }_{c,\mathrm{high}-\mathrm{U}}^{(\mathrm{approx})}(i{\epsilon }_n)$$

$$\rho (\omega )=\frac{|\omega |}{2t^{2}D}{\mathrm{e}}^{-{\omega }^2/(2t^{2}D)},$$

$${\Sigma }_{c,\mathrm{low}-\mathrm{T}}^{(\mathrm{approx})}(\epsilon )=\sigma (0,1,1,\epsilon )+\sigma (\delta ,z_1,y_1,\epsilon )$$

$${\Sigma }_{c,\mathrm{high}-\mathrm{T}}^{(\mathrm{approx})}(\epsilon )=\sigma (0,1,1,\epsilon )=\frac{\Gamma }{\epsilon +i\Delta }$$

$${\Delta }_c\sim T_K^{(\mathrm{scaling})}=D(NJ{\rho }_c{)}^{1/N}{\mathrm{e}}^{-1/(NJ{\rho }_c)}$$

$$\mathrm{Re}\sigma (\nu )$$

$$C_N(n)=\frac{A_N(n)}{B_N(n)}=\frac{a_1}{1+\frac{a_2(n-n_1)}{1+\frac{a_3(n-n_2)}{1+...}}},$$

$$\mathrm{for}r=1:g_1(n_i)=G(n_i),$$

$$\mathrm{for}r>1:g_r(n)=\frac{g_{r-1}(n_{r-1})-g_{r-1}(n)}{(n-n_{r-1})g_{r-1}(n)},$$

$$A_{q+1}(n)$$

$$B_{q+1}(n)$$

$${\overline{G}}_f(i{\epsilon }_n)=\int d\omega \rho (\omega ){\left(i{\epsilon }_n-{\epsilon }_f-{\Sigma }_f(i{\epsilon }_n)-\frac{V^2}{i{\epsilon }_n-\omega +\mu }\right)}^{-1}$$

$${\Sigma }_c(i{\epsilon }_n)=\frac{V^2}{i{\epsilon }_n-{\epsilon }_f+{\Sigma }_f(i{\epsilon }_n)},$$

Full text

I Introduction

The physics of strongly correlated electron systems is very rich due to the interplay and competition of different energy scales. In the case of rare-earth systems, the Kondo lattice model (KLM) is the simplest attempt to investigate heavy electrons as a consequence of collective Kondo effect. The easiest way to study the KLM is the application of dynamical mean-field theory (DMFT), although it becomes exact only at infinite dimension. Quantum Monte Carlo (QMC) approaches are essential computation methods in condensed matter physics, and are often used as impurity solver for DMFT calculations as well.

In QMC methods that evaluate imaginary time path integrals quantities like Green’s functions and response functions are obtained in imaginary time. In order to interpret physical quantities measured experimentally at real energies, the Matsubara frequency ($i\varepsilon_{n}$) data obtained by Fourier transform from the imaginary-time results have to be analytically continued to real energies, $\varepsilon$. The subject of analytic continuation is an important issue in numerical physics, however, it is very difficult because of its ill-posed nature. When the function ${\cal F}(i\varepsilon_{n})$ under question is a rational function of $i\varepsilon_{n}$, the replacement $i\varepsilon_{n}\rightarrow\varepsilon+i\delta$ can be done. However, it is not correct when ${\cal F}(i\varepsilon_{n})$ is not a rational function. Furthermore, from the ill-posed nature it follows that the output depends sensitively on the input, therefore, statistical errors that are always present in numerical data make the analytic continuation scheme instable. The simplest and widely used approach to perform analytic continuation is the Padé approximationPade where the real-energy function in question is approximated by a continued fraction using the Matsubara frequency data points. Other methods are also often used for analytic continuation in numerical physics such as the Maximum Entropy methodMaxEnt .

Efficient QMC impurity solvers such as the continuous-time quantum Monte Carlo (CT-QMC)CTQMC method provide very accurate numerical results in imaginary time. Therefore, the simple Padé method for analytic continuation can be successfully applied for most impurity problems. This widely accepted method is very popular because of its simplicity or easy programmability, and it is often used also for less efficient, conventional QMC methods. However, it has disadvantages as well including the overestimation of statistical errors or cut-off effects due to the limited number of Matsubara data points. Another problem would be that when we solve lattice problems by DMFT with considerable number of self-consistency steps, less accurate results are derived compared to impurity cases with a given amount of computational time even when we use efficient QMC methods. Thus, the application of Padé method as analytic continuation is often questionable.

In this paper we present a modified Padé approach which is based on the analytic continuation of the self-energy instead of the Green’s function. Namely, we apply the Padé analytic continuation method for the conduction electron self-energy known at imaginary-frequency points, and then calculate the Green’s functions and $t$-matrix at real energies from the analytically continued, real-energy self-energy. We show that this modified, self-energy Padé approach works better than the standard way based on the direct analytic continuation of the Green’s function especially for systems that have bare conduction electron density of states (DOS) with singular behavior at band edges or around the Fermi level. We found that in Ref. Wang, while the effect of antiferromagnetic ordering was studied on the properties of a Mott insulator, the analytic continuation of the self-energy was proposed. However, the authors neither used the Padé method for analytic continuation nor discussed the properties and efficiency of the self-energy analytic continuation. We discuss the reasons for the success of the modified Padé approach including its robustness against statistical errors or stability by actual calculations for the non-interacting Anderson lattice model, and also for the Hubbard model. Finally we apply the modified, self-energy Padé method for the zero-gap Kondo lattice with linearly vanishing bare conduction electron DOS at the Fermi level as well, and discuss the properties of its Kondo insulating state.

The so-called pseudogap or zero-gap Kondo impurity models with vanishing conduction electron DOS at the Fermi level as $\rho(\omega)\sim|\omega|^{r}$ have been extensively studiedWithoff ; Chen ; Fritz_2004 ; Bulla1 ; Buxton ; Vojta ; Vojta_Fritz ; Bulla_Glossop ; Ingensent in the last two decades. Renormalization group calculations have given comprehensive analytical account for the fixed-point structure, phase transitions and critical propertiesWithoff ; Chen ; Fritz_2004 ; Buxton ; Vojta_Fritz , while numerical renormalization group (NRG) studies have provided numerical calculations of dynamic and thermodynamic quantitiesIngensent ; Bulla_Glossop ; Bulla1 ; Vojta . General understanding for the pseudogap Kondo impurity model involves a quantum phase transition at a critical Kondo coupling $J_{c}$ above which the impurity is Kondo screened while it is in a local moment state for $J<J_{c}$. For $r\geq 1/2$ the transition disappears for particle-hole symmetric case and the impurity is always unscreened independent of the value of Kondo coupling $J$. However, with particle-hole asymmetry the phase transition is present for all $r>0$. In contrast to the impurity version of the zero-gap Kondo model, less is known about its lattice version. The Dirac Kondo lattice model as the extension of the impurity pseudogap Kondo model with exponent $r=1$ has been studied by mean-field methodsPrincipi ; Zhong and QMC simulationsLiu ; Fu . Refs. Zhong, ; Liu, ; Fu, have investigated the phase diagram of the Dirac Kondo lattice with possible electronic orders by extended mean-field studyZhong and QMC methodLiu ; Fu , however, less attention has been paid to the Kondo insulating state. In our study on the zero-gap Kondo lattice with $r=1$ we concentrate on its insulating phase. This $r=1$ case might have relevance for the case of Kondo physics in graphene or Kondo correlation in $d$-wave superconductorsFritz . Based on the accurately derived real-energy Green’s function and $t$-matrix obtained by the modified Padé approach, we discuss dynamical and thermodynamical properties of the Kondo insulating phase including the optical conductivity, static local susceptibility, dependence of the insulating gap on the Kondo coupling $J$, and coherence effects. We identify two energy scales as well from the dynamic and thermodynamic properties, that are associated as a direct and an indirect gap in a band hybridization picture.

II The modified analytic continuation method

In QMC methods that evaluate imaginary-time path integrals the conduction electron Green’s function is obtained in the imaginary time domain. Interpretation of physical quantities measured experimentally requires the knowledge of Green’s function at real energies, which is mathematically carried out by analytic continuation of the Matsubara data down to the real-energy axis by the procedure $i\varepsilon_{n}\rightarrow\varepsilon+i\delta$.

Numerically the conduction electron DOS, $\rho_{c}$, is ordinarily derived at real energies, $\varepsilon$, through the direct analytical continuation of the Matsubara Green’s function $G_{c}(i\varepsilon_{n})$ as

[TABLE]

which will be addressed as Green’s function approach, or direct analytic continuation method, in the following.

Instead of this direct way, there is a modified approach to obtain the spectral density $\rho_{c}(\varepsilon+i\delta)$ by using the analytically continued conduction electron self-energy, $\Sigma_{c}(\varepsilon+i\delta)$, derived from the Matsubara data $\Sigma_{c}(i\varepsilon_{n})$. Namely, the Green’s function $G_{c}(\varepsilon+i\delta)$ is expressed by the analytically continued self-energy at real energies as

[TABLE]

with $\rho(\omega)$ being the bare conduction electron DOS, which gives $\rho_{c}$ as

[TABLE]

for ${\rm Im}\,\Sigma_{c}(\varepsilon)\gg\delta$, and we took $\mu=0$ for simplicity. Simplifying further the notation, we take the convention in this paper that energy including an infinitesimal imaginary part, $\varepsilon+i\delta$, is simply written as $\varepsilon$. We will address the calculation of $G_{c}(\varepsilon)$ or $\rho_{c}(\varepsilon)$ from $\Sigma_{c}(\varepsilon)$ as self-energy approach, or modified analytic continuation method, in the following.

So far, both the Green’s function and self-energy approaches were introduced generally. From now on we apply these two approaches with the Padé analytic continuation method, whose details is given in Appendix A. As a test for the modified Padé approach we take the semicircular bare conduction electron DOS in this section given as

[TABLE]

with has sharp change at band edges, and examine the properties of the self-energy Padé approach through the example of the non-interacting Anderson lattice and Hubbard model.

II.1 Exemplary case of the non-interacting Anderson lattice model

II.1.1 Padé approximation

First we take the non-integrating Anderson lattice model (ALM) and discuss the main features of the Green’s function and self-energy Padé approaches. Since this model can be solved exactly, we can compare the numerically obtained resultsnote-ctqmc with statistical errors to the exact ones. Namely, first we derive the local ($G_{f}$) and conduction electron ($G_{c}$) Green’s functions together with the self-energy $\Sigma_{c}$ numerically in the imaginary-time domain within DMFT. Then, we obtain the local and conduction electron spectral densities by the Green’s function and self-energy Padé analytic continuation approaches and compare them to the exact solutions.

The non-interacting ALM is given by the Hamiltonian

[TABLE]

where $f_{i,\sigma}$ ($f^{{\dagger}}_{i,\sigma}$) and $c_{{\mathbf{k}},\sigma}$ ($c^{{\dagger}}_{{\mathbf{k}},\sigma}$) are annihilation (creation) operators for the local and conduction electrons with the cite $i$ and spin $\sigma$ indices. The dispersion $\varepsilon_{{\mathbf{k}}}$ of the conduction electrons in Eq. (5) corresponds to the bare conduction electron DOS as $\rho(\omega)=1/N_{0}\sum_{\mathbf{k}}\delta(\omega-\varepsilon_{{\mathbf{k}}})$ with $N_{0}$ being the number of lattice sites, where we take the semicircular conduction electron DOS given in Eq. (4) with $D=1$ in the numerical calculations. Analytic results for the non-interacting ALM that we will use in the following are summarized in Appendix B.

Figure 1 shows the numerically calculated conduction electron and $f$-electron DOS at real energies derived by the Green’s function Padé approach (left panel) and self-energy Padé approach (right panel) together with the exact solution. It is clear that the real-energy results obtained by the self-energy Padé approach show excellent agreement with the exact solution in spite of the presence of statistical errors, while in the case of the Green’s function Padé approach the obtained densities of states are less accurate.

Now we discuss the characteristics of the two Padé approaches. Considering the self-energy Padé approach first, we start from the exact Matsubara self-energy given as

[TABLE]

from Eq. (32) by taking $\varepsilon_{f}=0$, which is a rational function. Therefore, it is immediately expected that the procedure $i\varepsilon_{n}\rightarrow\varepsilon+i\delta$ is correct and leads to the real-energy result

[TABLE]

which gives

[TABLE]

obtained after the $\omega$-integration in Eq. (3), as

[TABLE]

i.e. the exact result given in Eq. (36) is recovered. Actually, following the Padé method described in Appendix A, the Padé coefficients $a_{i}$ for the self-energy are calculated as

[TABLE]

by taking Eq. (6) for $\Sigma_{c}(i\varepsilon_{n})$, which gives the real-energy Padé approximant as

[TABLE]

as expected from Eq. (7). The rational function behavior of the self-energy guarantees the success of the modified Padé approach as it was demonstrated in the right part of Fig. 1.

Considering next the Green’s function Padé approach, we start from the imaginary-frequency Green’s function given as

[TABLE]

from Eq. (33) by taking Eqs. (4) and (6), which is not a rational function. Thus, we cannot simply use the replacement $i\varepsilon_{n}\rightarrow\varepsilon+i\delta$, which implies that the exact result for $\rho_{c}(\varepsilon)$ given in Eq. (36) will not be fully reproduced. Indeed, when we consider the coefficients $a_{i}$ analytically for the Padé approximant of the Green’s function, we find that all coefficients $a_{i}$ are non-zero in contrast to the case of the self-energy. In left part of Fig. 2 we show the real-energy Padé approximants for the Matsubara Green’s function given in Eq. (12) with increasing number of Matsubara frequency points $N$ used in the calculation. As we increase the number $N$, which is the only adjustable parameter in the Padé method, the exact solution is approximated more and more. However, we find that for $N$ larger than about $1000$ the result does not change apparently, and this converged Padé result shown separately in the right part of Fig. 2 differs from the exact result, especially at the sharp band and gap edges.

II.1.2 Effect of statistical errors

So far we did not considered the presence of statistical noises in the analytic calculations, however, numerical calculations always contain statistical errors which are typically overestimated in the real-energy Padé approximant giving rise to inaccuracy in the Padé method. In order to examine this effect, we model the presence of statistical noises in the conduction electron self-energy and Green’s function by introducing random deviations from the exact results $\Sigma_{c}(i\varepsilon_{n})$ and $G_{c}(i\varepsilon_{n})$ for the ALM given in Eqs. (6) and (12), respectively, at each Matsubara frequency point as

[TABLE]

Here, $\delta_{r}$ is a random number in the range $\{-1,1\}$ and $q$ is a control parameter which determines the "largeness" of the noise (error) in this toy model. In Eqs. (13) and (14) we add error only to the imaginary part of $\Sigma_{c}$ and $G_{c}$ because they are pure imaginary due to the presence of particle-hole symmetry.

We take the functions $\Sigma_{c}^{(\rm error)}(i\varepsilon_{n})$ and $G_{c}^{(\rm error)}(i\varepsilon_{n})$ with noises and derive the analytically continued self-energy and Green’s function by the Padé method. Figure 3 shows the spectral densities calculated from $\Sigma_{c}^{(\rm error)}(\varepsilon)$ and $G_{c}^{(\rm error)}(\varepsilon)$ for different values of the error parameter $q$. We find that increasing the largeness of the random noises by increasing the value of parameter $q$ acts differently on the Green’s function and self-energy Padé approaches. Namely, the self-energy Padé approach is less sensitive to the error compared to the Green’s function approach, and remains stable even for "large" errors. On the contrary, even small error causes large instability of the Padé approximant in the Green’s function method as it can be seen in the upper panel of Fig. 3. This toy model with $q\sim 10^{-6}$ depicts a situation similar to the numerical solution of the ALM shown in Fig. 1.

II.2 Application of the modified Padé method to the Hubbard model

In this Section we take another example which is the half-filled Hubbard model given by the Hamiltonian

[TABLE]

and obtain the spectral density $\rho_{c}(\varepsilon)$ by both Padé approaches numericallynote-ctqmc within DMFT. The bandwidth is $D=2t$ in Eq. (15), and we take the semicircular bare conduction electron DOS given in Eq. (4).

Figure 4 shows the spectral densities as a function of real energy for Coulomb interaction values $U=2$ and $U=2.8$ obtained by both Padé methods, while result for $U=3.5$ is shown in the inset of Fig. 5. We intentionally chose relatively small number of Monte Carlo steps in the simulations in the left and center panels of Fig. 4 in order to test the accuracy of the two Padé approaches. In the right panel of Fig. 4, on the other hand, we show result with larger number of Monte Carlo steps compared to the center panel as comparison. For $U=2$ the ground state is a correlated metal which is indicated by a quasiparticle peak at the Fermi level. As we increase the value of $U$, the quasiparticle peak narrows, as for the case of $U=2.8$ in Fig. 4, and finally disappears in the insulating phasenote-3 for which an example is shown in Fig. 5 with $U=3.5$. The position of the metal-insulator transition is estimated as $U_{\rm cr}/D=2.6$ by iterated perturbation theorydmft , while it is found as $U_{\rm cr}/D=2.94$ from DMFT calculation solved by NRGBulla_Hubbard_NRG for the Bethe lattice (semicircular DOS). We find $U_{\rm cr}/D=2.9\pm 0.05$ in our numerical calculations from the sharp change in the wave-function renormalization factor as a function of $U$, which coincides with the disappearance of the quasiparticle peak in the spectra.

The discrepancy between the results for the spin-up and spin-down components shown by dashed and continuous lines in Fig. 4 mainly reflects the presence of statistical errors. We find again that the self-energy Padé approach is less sensitive for the statistical errors while they cause large instabilities in the case of the Green’s function Padé approach. Moreover, the Green’s function approach gives even unphysical negative value for the spectral density as it can be seen in the inset of the center panel of Fig. 4. When the number of the Monte Carlo steps is increased, the difference between the results obtained by the two Padé approaches reduces.

To investigate the success of the self-enegy Padé approach for the Hubbard model, we examine the properties of the conduction electron self-energy. Fig. 5 shows the analytically continued self-energy both in the correlated metallic ($U=2$) and insulating ($U=3.5$) phases. Since the transition between these two phases is first-order, the self-energy has different shape in the two phases. We find that the self-energy shows smooth energy dependence in both phases and can be described well phenomenolocically. Namely, the phenomenological form

[TABLE]

describes the self-energy in the correlated metallic phase with two-peak structure in its imaginary part, while formula

[TABLE]

describes it in the insulating phase with a one-peak structure in the imaginary part, where the function $\sigma$ is defined as

[TABLE]

The fits of the numerical data with these formulas are also shown in Fig. 5.

Since the $\varepsilon$-dependence of the self-energy can be excellently approximated by rational function in both phases, i.e. with formulas (16), (17), it is expected that the Padé approximant for the self-energy obtained from the Matsubara numerical data will be accurate. Actually, taking the Matsubara phenomenological forms

[TABLE]

from Eqs. (16), (17) by the inverse of the replacement $i\varepsilon_{n}\rightarrow\varepsilon+i\delta$, the mathematically calculated real-energy Padé approximants quickly converge and fully recover the real-energy phenomenological forms (16) and (17) as it is shown in Appendix C.

III Zero-gap Kondo lattice model

Last, we consider a zero-gap Kondo lattice model with bare conduction electron DOS that linearly vanishes at the Fermi level and solve it numericallynote-ctqmc within DMFT. This model corresponds to the lattice version of the pseudogap Kondo modelWithoff ; Chen ; Fritz_2004 ; Bulla1 ; Buxton ; Vojta ; Vojta_Fritz ; Bulla_Glossop ; Ingensent with exponent $r=1$ which might have relevance for Kondo physics in grapheneFritz . Since DMFT is an infinite dimensional theory, we take the infinite dimensional extension of the two-dimensional honeycomb and the three-dimensional diamond lattice as DOS, which is given as

[TABLE]

where we take the bandwidth $D=1$ as the unit of energy. This DOS has an extra parameter $t$ compared to the well-studied case $\rho(\omega)\sim|\omega|^{r}$, which controls the slope of the linear dispersion, i.e. the density of conduction electrons being around the Fermi level.

First, we compare spectral densities calculated at real energies by the Green’s function and self-energy Padé approaches shown in Fig. 6 at selected values of the Kondo coupling $J$, $t$, and temperature in the Kondo insulating phase. The insulating gap is smaller for $t=1$ than for $t=0.5$ since the slope of the linear dispersion, which is given as $1/(2t^{2}D)$ from Eq. (21), is smaller for the $t=1$ case. As a comparison, we show also example for the ordinary Kondo model with semicircular, i.e finite conduction electron DOS at the Fermi level. We find in all cases again that the modified Padé approach is more stable than the Green’s function Padé approach, which is indicated by the almost complete overlap of the spin-up and spin-down components. Although both methods give similar result near the Fermi level, the Green’s function Padé approach becomes uncertain at larger energies.

Figure 7 shows the analytically continued self-energy as a function of real energy both in the high- and low-temperature phases. The imaginary part of the self-energy shows a simple one-peak structure around the Fermi level at high temperatures, while two additional peaks evolve at higher energies with DOS origin as the temperature is decreased. Furthermore, the central peak narrows with decreasing temperature which indicates the continuous evolution of a small energy scale responsible for the Kondo effect. We find again that the smooth energy dependence of $\Sigma$ can be described well by the phenomenological form

[TABLE]

in the low-temperature phase, and as

[TABLE]

in the high-temperature phase, where the function $\sigma$ is defined in Eq. (18). These fits are also included in Fig. 7. We argue that the rational function nature of expressions (22) and (23) is responsible for the success of the self-energy Padé method similarly to the case of the Hubbard model.

We note that the conduction electron self-energy has different origin in the Kondo lattice and in the Hubbard model, although the DMFT equations are the same in these two cases. In the case of the Kondo model the self-energy arises from the impurity scattering while in the Hubbard model it emerges because of the Coulomb interaction.

Figure 8 shows the temperature dependence of the conduction electron DOS and $t$-matrix at real energies in the Kondo insulating phase, i.e. where an energy gap develops as the temperature is decreased. The $J$-dependence of this insulating gap, $\Delta_{c}$, is shown in the left part of Fig. 9 for different values of the DOS parameter $t$. The gap $\Delta_{c}$ is measured at the half-peak value of $\rho_{c}(\varepsilon)$ since it becomes almost temperature independent at low temperatures. We find a linear dependence of $\Delta_{c}$ on the Kondo coupling $J$ for all values of parameter $t$, which is in contrast to the normal Kondo case with finite DOS at the Fermi level. Namely, we show also the gap for the Kondo lattice with semicircular bare conduction electron DOS in Fig. 9 as comparison, which is in agreement with the scaling formula for the Kondo temperature as

[TABLE]

obtained for the $N$-fold degenerate Coqblin-Schrieffer modelhewson with $D$ being the band half-width and $\rho_{0}$ being constant conduction electron DOS. In accord with the scaling formula (24) Kondo insulating state appears for any value of the Kondo coupling $J$ at low enough temperature due to the finite DOS at the Fermi level. In contrast, there is a critical value $J_{c}$ of the Kondo coupling for the formation of the low-temperature insulating phase for the zero-gap KLM as it is apparent from Fig. 9. The critical values $J_{c}$ are shown in the inset of the left part of Fig. 9 for different $t$ values, which are derived from the linear fit to the gap values $\Delta_{c}$. In the limit of $t\rightarrow 0$ the critical coupling $J_{c}$ tends to zero, i.e. the case of finite DOS at the Fermi level is recovered since the slope of the dispersion is infinite in the $t=0$ limit.

The temperature dependence of the static local susceptibility is shown in the left part of Fig. 10 for different values of the Kondo coupling $J$. For $J<J_{c}$ the susceptibility increases with decreasing temperature approximately as $\sim 1/T$, while it saturates at low temperatures for $J\geq J_{c}$. The critical Kondo coupling is estimated as $J_{c}\approx 0.23$ for $t=0.5$ in our numerical calculations (see left part of Fig. 9). We define the Kondo temperature from the susceptibility as $T_{\rm K}^{-1}=\chi(T\rightarrow 0)/C_{N}$, where $C_{N}$ is the Curie constant. Beside $T_{\rm K}$, another energy scale $T_{\rm c}$ can be identified which is the temperature where the susceptibility has maximum in its $T$-dependence. Below this temperature, coherence between the impurities appear in a Kondo lattice by coherent scattering off the Kondo singlets, which is shown by the behavior $\chi\sim\,{\rm constant}$ for $T<T_{\rm c}$. The $J$-dependence of the coherence temperature $T_{\rm c}$ is shown in the right part of Fig. 9.

We expect that coherence effects are reflected in further thermodynamic and dynamic quantities as well. Actually we find that the temperature values where the insulating gap appears in the spectral density at a given value of $J$ coincide well with $T_{\rm c}$ as it is shown in Fig. 9. Furthermore, if we obtain the inflextion pointsnote-2 $\{J_{\rm infl},T_{\rm infl}\}$ of the wave-function renormalization factor $z=(1-\partial{\rm Re}\Sigma_{c}(\varepsilon)/\partial\varepsilon|_{\varepsilon=0})^{-1}$ from its $J$- and $T$-dependence, and plot $T_{\rm infl}$ as a function of $J$, we find that this curve also coincides with the previous ones as it can be seen in Fig. 9. All these quantities are finite only above a critical value of the Kondo coupling which matches well with the critical value $J_{c}\approx 0.23$ obtained previously from the vanishing of the insulating gap $\Delta_{c}$. On the other hand, there is no such critical value for the case of the normal Kondo lattice model as it is shown in the inset of the right part of Fig. 9.

We note that in contrast to the lattice case, there is no Kondo effect in the corresponding zero-gap Kondo impurity problem even for infinitely strong Kondo coupling, which is illustrated by the temperature dependence of the static local susceptibility shown in the right part of Fig. 10. Namely, the susceptibility follows $1/T$ dependence with a $J$-dependent constant as the temperature is lowered for even large values of $J$ (e.g. $J=0.85$ in Fig. 10), which reflects that the impurity is in local moment state without screening. Kondo screening of the local moment in the impurity case is possible in particle-hole asymmetric case which can be achieved by the introduction of non-zero chemical potential, for example, or including potential scattering.

Right part of Fig. 11 shows the temperature evolution of the optical conductivity as a function of frequency, where the optical conductivity is calculated as pruschke

[TABLE]

with $f(\omega)=1/(1+{\rm exp}(\beta\omega))$ being the Fermi function. We can observe that a peak emerges at a finite frequency as the temperature is decreased which indicates the presence of a direct gap in the energy spectrum. The positions of these peaks, that we identify as direct gap, are shown in the inset as a function of Kondo coupling $J$ for different values of parameter $t$. Beside this direct gap, we identify the insulating gap $\Delta_{\rm c}$ from dynamics as indirect gap. Right part of Fig. 11 shows this indirect gap together with $\Delta_{t}$ obtained from the $t$-matrix, and the Kondo temperature $T_{\rm K}$ as well.

In order to interpret the direct and indirect gaps as energy scales, we recall the non-interacting Anderson lattice model. From the exact energy spectrum given in Eq. (39) for the ALM, these two energy gaps can be defined. The indirect gap is written as $\Delta_{\rm indirect}=E_{+}(k\rightarrow\infty)-E_{-}(k\rightarrow-\infty)\approx 2V^{2}/D$ that emerges due to the hybridization in dynamical quantities such as conduction electron or local DOS. The direct gap is expressed as $\Delta_{\rm direct}=E_{+}(k=0)-E_{-}(k=0)\approx 2V$. We recall Eq. (6) which says $\Sigma_{c}(i\varepsilon_{n})=V^{2}/(i\varepsilon_{n})$ for the self-energy in the case of the non-interacting ALM. Thus, we obtain the associated term "$V^{2}$" for the KLM by extracting numerically the coefficient of $1/\varepsilon_{n}$ in the self-energy $\Sigma_{c}(i\varepsilon_{n})$ in the Kondo lattice model, and with the use of this associated term "$V^{2}$" we estimate the gaps as $\Delta_{\rm indirect}\sim 2V^{2}/D(=2V^{2})$ and $\Delta_{\rm direct}\sim 2V$ which are also shown in Fig. 11. We find an excellent agreement between $2V$ and the direct gap of KLM obtained from the optical conductivity data as it is shown in the inset of the left part of Fig. 11. For the indirect gap we expect that $\Delta_{c}/2\sim V^{2}/D$, but instead we find from the numerical analysis that $\Delta_{c}/2\sim 0.5(V^{2}/D)$ for parameter value $t=0.5$ as it is shown in the right part of Fig. 11. The reason for the discrepancy might be that our identification between the indirect gaps of the non-interacting ALM and KLM is just approximate and does not contain interaction effects which are expected to influence the indirect gap stronger than the direct gap. Nevertheless, the hybridized bands picture with direct and indirect gaps seems to be held for the zero-gap Kondo lattice model as well.

IV Summary

In this paper we have presented and discussed a modified Padé approach for numerical methods in which quantities like correlation functions are obtained in imaginary time. The modified Padé approach is based on the analytic continuation of the self-energy which is then used for deriving the Green’s function and $t$-matrix at real energies. We showed that this modified, self-energy Padé approach is more stable and robust against statistical errors than the conventional way with the direct analytic continuation of the Green’s function, where the latter property is demonstrated within a toy model. We found that the success of the modified Padé approach lies in the fact that the self-energy has smooth real-energy dependence and this dependence can be approximated excellently by rational function, therefore, the Padé approximation works well for the self-energy. We considered also the zero-gap Kondo lattice with linearly vanishing conduction electron DOS at the Fermi level where we applied the modified Padé approach, and investigated the properties of its Kondo insulating phase. We found linear dependence of the insulating gap on the Kondo coupling $J$ with a finite critical value $J_{c}$ above which Kondo insulating state appears. Besides the Kondo temperature $T_{\rm K}$, we identified the coherence temperature $T_{\rm c}$ as well from thermodynamic and dynamic quantities below which coherence between the impurities appear. Finally, we have identified two additional energy scales from dynamic and thermodynamic properties, that are associated as a direct and an indirect gap in a band hybridization picture.

Acknowledgment

The author would like to acknowledge helpful advices and discussion with Professor Yoshio Kuramoto. The author acknowledges the Bolyai Program of the Hungarian Academy of Sciences and the Hungarian Scientific Research Funds No. K124176.

Appendix A Padé approximation

The principal problem in analytic continuation by using Padé approximation is to find a rational function $C_{N}(n)$ which approximates the function $G(n)$ at arbitrary real energies $n$, $G(n)\sim C_{N}(n)$, on condition that we know the function $G$ at $N$ Matsubara frequency points $n_{1},n_{2},...,n_{N}$. The rational function $C_{N}(n)$ is expressed as

[TABLE]

where the coefficients $a_{i}$ are determined from the condition $C_{N}(n_{i})=G(n_{i})$ for $i=1,...,N$.

The simplest way for obtaining the coefficients $a_{i}$ is the use functions $g$ as $a_{i}=g_{i}(n_{i})$, where $g_{r}(n)$ can be obtained by the following recursion formula:

[TABLE]

which gives the functions $A_{N}(n)$ and $B_{N}(n)$ in Eq. (26) as

[TABLE]

with the starting conditions $A_{0}=0$, $A_{1}=a_{1}$, $B_{0}=B_{1}=1$.

Appendix B Non-interacting Anderson lattice model

In the DMFT treatment of the Anderson lattice model given in Eq. (5) the site-diagonal Green’s function is calculated as

[TABLE]

in the Matsubara frequency domain with $\varepsilon_{n}=(2n+1)\pi T$, and $\rho(\omega)$ is the bare conduction electron DOS of the conduction electrons. We quote also the results for the conduction electron self-energy and Green’s function, which read as

[TABLE]

and

[TABLE]

In the non-interacting case, the local self-energy vanishes, i.e. $\Sigma_{f}=0$, and thus from Eqs. (31), (32) and (33) the local and conduction electron Green’s functions can be evaluated exactly in the Matsubara frequency domain as

[TABLE]

The conduction and local electron DOS are obtained at real energies as

[TABLE]

after performing the ${\mathbf{k}}$-summations, which can be expressed by the self-energy as well. Namely, $\rho_{c}(\varepsilon)$ was already expressed in Eq. (3), and the local DOS is expressed as

[TABLE]

where we used Eqs. (31) and (32).

The energy of the hybridized bands between the local and the conduction electron states in the non-interacting Anderson lattice model is expressed as

[TABLE]

Appendix C Padé approximant for the self-energy in case of the Hubbard model

We calculate the Padé approximant for the phenomenological self-energy both for the metallic and insulating phases given in Eqs. (19) and (20). From Eq. (20) we calculate the Padé coefficients as

[TABLE]

based on Appendix A, which gives the real-energy Padé approximant as

[TABLE]

i.e. the phenomenological form (17) is reproduced.

In the correlated metallic phase the Padé coefficients $a_{i}$ are calculated for Eq. (19) as

[TABLE]

With these coefficients the analytically continuated $\sigma$ function is obtained as

[TABLE]

Thus the real-energy phenomenological form (16) is also reproduced.