Finite-temperature linear conductance from the Matsubara Green function
without analytic continuation to the real axis
We illustrate how to calculate the finite-temperature linear-response conductance of quantum impurity models from the Matsubara Green function. A continued fraction expansion of the Fermi distribution is employed which was recently introduced by Ozaki [Phys. Rev. B 75, 035123 (2007)] and converges much faster than the usual Matsubara representation. We give a simplified derivation of Ozaki’s idea using concepts from condensed matter theory and present results for the rate of convergence. In case that the Green function of some model of interest is only known numerically, interpolating between Matsubara frequencies is much more stable than carrying out an analytic continuation to the real axis. We demonstrate this explicitly by considering an infinite tight-binding chain with a single site impurity as an exactly-solvable test system, showing that it is advantageous to calculate transport properties directly on the imaginary axis. The formalism is applied to the single impurity Anderson model, and the linear conductance at finite temperatures is calculated reliably at small to intermediate Coulomb interactions by virtue of the Matsubara functional renormalization group. Thus, this quantum many-body method combined with the continued fraction expansion of the Fermi function constitutes a promising tool to address more complex quantum dot geometries at finite temperatures.
pacs:71.27.+a, 72.10.-d, 73.21.La
Computing the finite-temperature linear conductance for widely used model systems like an Anderson impurity connected to non-interacting leads involves integrals over the impurity single-particle spectral function multiplied by the derivative of the Fermi distribution.pinkbook () At , however, most quantum many-particle methods are set up in imaginary time and yield the relevant one-particle propagator at the Matsubara frequencies on the imaginary axis. Thus, an analytic continuation is required in order to obtain the corresponding spectral functions on the real axis and to compute transport properties at finite temperatures. If is only known numerically – which is the generic scenario – the analytic continuation can in principle be achieved, e.g., by Padé approximation of the Matsubara data,pade1 (); pade2 () but in general such a procedure is ill-controlled (see, e.g., Ref. frequenzen, ). Alternatively, one can use the analytic properties of the Fermi function and resort to the calculus of residues to perform the integration directly on the imaginary axis. Even though this might involve the propagator (or, aiming at the conductance, its derivative) not only at the physical Matsubara frequencies but also at different arguments, computing such data by interpolation is expected to be significantly more stable than carrying out an analytic continuation.
In order to rotate the path of integration from the real to the imaginary axis, one can employ the Matsubara expansion for the Fermi distribution , but this is not useful for practical applications in fermionic many-body problems because of its slow convergence.FW () A much more efficient representation was recently proposed by Ozaki.Ozaki () It is based on a continued fraction expansion and ultimately yields an infinite sum over simple poles which are all located on the imaginary axis. In its infinite form, this series is equivalent to the usual Matsubara expansion. Truncation, however, modifies the position of the poles (but they remain on the imaginary axis) as well as the corresponding residues, and it turns out that the resulting finite series yields an excellent approximation to for not too large values of , even if only a few terms are accounted for (all this will be quantified below).
The motivation for the present work is three-fold. First, we give a simplified derivation of Ozaki’s idea which relies on concepts taken from condensed matter theory and present results for the rate of convergence of the truncated series (Section II). Secondly, we employ an exactly-solvable model – a non-interacting tight-binding chain with a single site impurity – and demonstrate that interpolation of the Matsubara Green function is much more stable than continuation to real frequencies, illustrating that it is advantageous to perform calculations on the imaginary axis (Section III). Thirdly, the formalism is applied to the single impurity Anderson model, and the Coulomb interaction is tackled by an approximation scheme which (in thermal equilibrium) is most easily set up in Matsubara space – the functional renormalization group (FRG).salmhofer (); dotsystems (); frequenzen (); severinsiam () By comparison with very accurate numerical renormalization group (NRG) data, we demonstrate that the FRG allows for reliably extracting the finite-temperature linear conductance at small to intermediate values of (Section IV). This was previously not possible due to the instability of the analytic continuation.frequenzen () Whereas within the NRG framework the numerical effort grows exponentially when considering general models with several correlated degrees of freedom and one is generically restricted to situations of high symmetry, the Matsubara FRG can easily be generalized to multi-dot or multi-level quantum dot geometries without specific symmetries with only a mild increase in computation time.frequenzen () Combined with the continued fraction expansion of the Fermi function, the FRG thus holds the promise for future applications addressing transport properties of such systems at .
Ii Continued-fraction representation of the Fermi function
In this Section we give a simplified derivation of Ozaki’s idea.Ozaki () The first step is to express the Fermi distribution in the form
where later on , with and denoting the inverse temperature and the chemical potential, respectively. In order to obtain a continued fraction expansion for , Ozaki introduced the (generalized hypergeometric) auxiliary functions
It is straightforward to show that
The first identity, e.g., is easily confirmed by establishing by induction. Thus, is determined by a ratio of the form , rendering it reasonable to seek for a recurrence relation. To this end, one calculates the difference of the power series of Eq. (2) for and :
Dividing by , taking the inverse and replacing by gives
where we have introduced the ratio
For a practical application it is not necessary to write down the continued fraction explicitly. Rather, one can directly exploit the fact that the very same recursion relation of Eq. (6) appears if one aims at computing the inverse of an infinite tridiagonal matrix associated with a semi-infinite one-dimensional tight-binding chain with nearest neighbor hopping and starting at site . Using the well-known Feshbach projection method,taylor () one obtains (for )
with the only non-vanishing matrix elements of (which decrease like for large ) given by
Thus, one finally faces the ordinary eigenvalue problem
All come in pairs owing to the fact that the state with components is an eigenvector of with eigenvalue . Then, Eq. (9) can be recast as
with , and symbolizing that the sum extends over the positive part of the spectrum (i.e., all for which ) only. This expansion of the Fermi function should be compared to the well-known () Matsubara sum decompositionFW (); Ozaki ()
which is slowly converging with w.r.t. . If the infinite matrix is not replaced by a finite one, Eqs. (12) and (13) should coincide, implying that and the positive eigenvalues given by . The fact that the latter accumulate at zero is obvious from the large decay of the matrix elements .
Ozaki pointed out the enormous acceleration of convergence one obtains for not too large arguments by replacing the infinite chain by one of finite length which corresponds to a finite termination of the infinite continued fraction (see also Ref. hm, ). In order to keep the property that all eigenvalues of come in pairs , one always truncates at an even number of sites ().
The eigenvalue problem posed by Eqs. (10) and (11) can be solved numerically using standard routines for all values of interest. It turns outOzaki (); hmfn () that the lowest percent of the are very close to the Matsubara values , and the associated residues are approximately one. The other inverse eigenvalues increase much faster than , and their residues become large. As , all poles of the function obtained from restricting the sum in Eq. (12) to terms are located on the imaginary axis. This is an important property for the applications discussed in Section IV which is not shared by another recent partial fraction decomposition of the Fermi function.Saalmann ()
Already for rather small values of , yields an excellent approximation for for not too large values of . This is illustrated in Fig. 1. Even at , the approximation to the Fermi distribution is excellent up to . Likewise, coincides with to the drawing accuracy for arguments . For very large values of , all functions eventually tend to .
In Fig. 2 we show the convergence of to one for very large negative values of . Results are shown for four different arguments, increasing in absolute value proportional to . By rescaling with , the four curves collapse, indicating that has to increase only as in order to obtain the same rate of convergence. In passing, we note that it might be possible to analytically analyze the truncation error along the lines of Ref. hm, .
Iii Stability of interpolation and continuation of numerical data
The linear-response finite-temperature conductance of the single impurity Anderson modelsiam () is determined by an energy integral of the derivative of multiplied by the local single-particle spectral function :wingreen ()
where we have chose units of (with being the elementary charge) and defined the modified Matsubara frequencies . The prefactor is a measure for the level-lead couplings (see below). The second equality of Eq. (14) follows from rotating to the imaginary axis using Eq. (12) – it will be derived explicitly in Section IV. If is only known numerically at the physical Matsubara frequencies , both and need to be determined, e.g., from computing the Padé approximationpade1 (); pade2 () to . Whereas it seems reasonable to expect interpolation (required for extracting ) to be rather stable, an analytic continuation (necessary to calculate ) is in general an ill-controlled procedure (see, e.g., Ref. frequenzen, ). In this Section, we address this issue explicitly by considering an infinite tight-binding chain with a single site impurity as an exactly-solvable test system. The latter constitutes a special realization of the single impurity Anderson model in absence of Coulomb interactions and is associated with a rich local density of states (featuring both a continuum as well as one bound state).
The infinite tight-binding geometry is governed by the Hamiltonian
where is the hopping amplitude between nearest neighbors, and denotes the strength of an impurity located at site . We assume the system to be in thermal equilibrium associated with temperature and zero chemical potential. Using standard projection techniques,taylor () the local impurity Green function can be calculated as
where is the local propagator of an isolated semi-infinite chain:
The imaginary part of the derivative of (which determines the linear conductance; see Eq. (14)) reads
and the retarded Green function can be calculated straightforwardly from Eq. (16) by replacing :
Thus, the local density of states is given by
It features a continuum as well as a single pole.
We can now test the stability of both the interpolation on the imaginary axis and continuation to real frequencies of the Matsubara Green function. To this end, we compute the Padé approximation (as outlined in Refs. pade1, ; pade2, ; frequenzen, ) to and from that both as well as and compare with the exact results of Eqs. (18) and (20), respectively. The Padé approximation is calculated numerically for different meshes of Matsubara frequencies determined by the temperature (as well as additional discretization parameterscommentdiscr ()). For our non-interacting tight-binding chain, both the (exact) local density of states and the derivative do not depend on .
The results for and obtained from the Padé approximation to the propagator of Eq. (16) at different temperatures are shown in Figs. 3 and 4, respectively. It turns out that the former agrees well with the exact curve and is in particular independent of . In contrast, the Padé data for depends strongly on temperature (as well as on other discretization parameterscommentdiscr ()) and is plagued by severe artifacts such as . As expected, analytic continuation of a numerically known Matsubara Green function is much more unstable than interpolating on the imaginary axis. This demonstrates the advantage of evaluating physical quantities (e.g., the linear conductance which is determined either by or ) directly in Matsubara frequency space.
Iv Functional RG for the single impurity Anderson model
In this Section, we apply the formalism introduced above to compute the linear conductance of the single impurity Anderson model (SIAM) at finite temperatures. To this end, we show how to express in terms of the Matsubara Green function using the continued fraction expansion of the Fermi function and demonstrate that in the exactly-solveable non-interacting (test) case, only a few poles () need to be accounted for in order to obtain the conductance in agreement with the one computed on the real axis. For finite (and not too large) Coulomb repulsions, the functional renormalization group – which in equilibrium is most easily implemented in Matsubara frequency space – allows for calculating in good agreement with NRG data. This was previously not possible due to the instability of the analytic continuation.commentfrg ()
iv.1 Model Hamiltonian and linear conductance from the imaginary axis
The Anderson model is governed by the Hamiltoniansiam ()
where and denote fermionic annihilation operators associated with (left and right) baths as well as a single impurity, respectively. The latter features an onsite energy and a Coulomb repulsion between electrons of different spin directions . The baths are assumed to be in the wide-band limit leading to an energy-independent hybridization .
Given the exact local single-particle spectral function , the linear-response conductance of the SIAM at temperature can be calculated aswingreen () (in units of )
where we have exploited that the retarded and advanced Green functions are analytic in the upper and lower half of the complex plane, respectively.
and the linear-response conductance as a function of the level position and temperature can be calculated (numerically) from Eqs. (22) and (23), respectively. Comparing both approaches provides another demonstration of the rapid convergence of the imaginary-axis expression based on the continued fraction expansion of the Fermi function. It turns out (see Fig. 5) that if the series of Eq. (23) is truncated with (rather small) values , the resulting conductance nicely agrees with that of Eq. (22) for arbitrary temperatures from to .
In presence of Coulomb interactions, ‘numerically exact’ results for the spectral function of the SIAM can be obtained from the NRG framework.nrg () If one aims at describing more complex geometries, however, it is desirable to devise approximate schemes. A recent approach is provided by the Matsubara functional renormalization groupsalmhofer () which re-formulates a given many-particle problem in terms of an infinite set of coupled flow equations for single-particle irreducible vertex functions with an infrared cutoff as the flow parameter. Truncation of this hierarchy renders the FRG approximate w.r.t. the two-particle interaction and can hence a priori be justified only in the limit of small . Despite this fact, application of the most simple truncation scheme – which yields flow equations for effective system parameters and can thus be regarded as a kind of ‘RG enhanced Hartree-Fock’ approach – to various quantum dot geometries in equilibrium turned out to give accurate results for the zero-temperature linear conductance even at fairly large Coulomb interactions.dotsystems () If one employs a more elaborate truncation (which accounts for the frequency dependence of the two-particle vertex), one can in principle compute the spectral function of the SIAM in agreement with NRG data for small to intermediate values of .frequenzen (); ralf () However, the latter requires an analytic continuation of the Matsubara Green function,commentfrg () which for the SIAM was observed to be particularly ill-controlled if both and .frequenzen () For this reason, it was previously not possible to address the linear conductance as a function of the level position at finite temperatures. The continued fraction expansion for the Fermi function now allows for calculating directly from the imaginary axis. To this end, we extract the FRG approximationfrgappr3 () to the Matsubara Green function of the SIAM using precisely the formalism outlined in Ref. frequenzen, . Thereafter, we compute the derivative from interpolating by virtue of a Padé approximation – which again turns out to be far more stable than continuation to the real axis – and ultimately the linear conductance by Eq. (23). The result is shown (and compared to NRG datanrg ()) in Fig. 6. The FRG correctly describes the widening of the Lorentzian lineshape of due to Coulomb correlations (signaling the development of a Kondo plateau) and for small (intermediate ) agrees nicely (decently) with the NRG reference. Moreover, the temperature-dependence of at particle-hole symmetry is purely quadratic (see the inset to Fig. 6(a)) as expected from Fermi-liquid theory.hewson () For larger values of , the agreement with the numerical RG data becomes worse (which is in line with the results of Ref. frequenzen, ).
In this Paper we illustrated how to compute the finite-temperature linear-response conductance of quantum impurity models from the Matsubara Green function using a rapidly-converging continued fraction expansion of the Fermi function recently introduced by Ozaki. In case that is only known numerically, this formalism allows to circumvent the need to carry out an (ill-controlled) analytic continuation to the real axis. As an application, we studied the single impurity Anderson model at finite temperatures within the framework of the Matsubara functional renormalization group and showed that can be obtained accurately in comparison with numerical RG data for not too large Coulomb interactions. The FRG therefore holds the promise for future treatments of more complex quantum dot geometries which cannot be easily addressed within the NRG framework.
We are grateful to the Deutsche Forschungsgemeinschaft for support via FOR723 and thank Theo Costi for providing his code to carry out NRG calculations. Useful discussions with Jan von Delft and Harmut Monien are acknowledged.
- (1) H. Bruus and K. Flensberg, Many-Body Quantum Theory in condensed Matter Physics (Oxford University Press, Oxford, 2004).
- (2) G.A. Baker Jr., Essentials of Padé Approximants (Academic Press, New York, 1975).
- (3) H.J. Vidberg and J.W. Serene, J. Low. Temp. Phys. 29, 3-4 (1977).
- (4) C. Karrasch, R. Hedden, R. Peters, Th. Pruschke, K. Schönhammer, and V. Meden, J. Phys.: Condensed Matter 20, 345205 (2008).
- (5) A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (McGraw-Hill, New York, 1971).
- (6) T. Ozaki, Phys. Rev. B 75, 035123 (2007).
- (7) M. Salmhofer, Renormalization (Springer, Berlin, 1998.)
- (8) C. Karrasch, T. Enss, and V. Meden, Phys. Rev. B 73, 235337 (2006).
- (9) S. Jakobs, M. Pletyukhov, and H. Schoeller, Phys. Rev. B 81, 195109 (2010).
- (10) J.R. Taylor, Scattering Theory (Wiley, New York, 1972).
- (11) H. Monien, Math. Comp. 79, 857 (2010).
- (12) In the limit the fraction of poles which agree well with their Matsubara counterparts approaches (H. Monien, private communication). In the context of Gaussian quadrature of sums the truncation of the continued fraction for plays a central role. Analytical results for this issue are presented in Ref. hm, .
- (13) A. Croy and U. Saalmann, Phys. Rev. B 80, 073102 (2009).
- (14) P. W. Anderson, Phys. Rev. 124, 41 (1961).
- (15) Y. Meir and N. Wingreen, Phys. Rev. Lett. 68, 2512 (1992).
- (16) For finite temperatures, the Matsubara frequencies are discrete but extend to infinity nevertheless. Thus, one needs to introduce a finite frequency mesh in order to numerically compute the Padé approximation to the Matsubara Green function. There are various possibilities for realizing such mesh (e.g., linearly or ‘logarithmically’ spaced groups of frequencies). While the derivative turns out to be independent of those details, the local density of states is rather ‘unstable’. The latter is illustrated in Fig. 4 where the mesh employed is the one detailed in Ref. frequenzen, . Following the notation of this Ref. frequenzen, , the discretization parameters of the upper panel are given by , , and , and those of the lower panel read , , and .
- (17) The symbols in Fig. 3 show the Padé approximation for the derivative only at the physical Matsubara frequencies . If one seeks to calculate the linear-response conductance by virtue of the continued fraction expansion of the Fermi function (e.g., using Eq. (14) for the SIAM), it is necessary to know at the modified frequencies . As discussed in Ref. Ozaki, , the first 60 percent of the almost coincide with their Matsubara counterparts, and it turns out that can be computed in a stable way for arbitrary arguments except those close to the real axis ().
- (18) If one implements the FRG in Keldysh space, the finite-temperature linear-response (as well as the non-equilibrium) conductance of the SIAM can be obtained without carrying out an analytic continuation.severinsiam () The structure of the FRG flow equations, however, is far more complex in the Keldysh than in the Matsubara formalism. If one is interested in linear response only, the latter thus provides the most simple framework to study more complex geometries.
- (19) Standard NRG is a well-established numerical tool to compute low-energy equilibrium properties of quantum impurity systems. A detailed introduction to this method can be found in Ref. thomasnrg, .
- (20) R. Bulla, T. Costi, and Th. Pruschke, Rev. Mod. Phys. 80, 395 (2008).
- (21) R. Hedden, V. Meden, Th. Pruschke, and K. Schönhammer, J. Phys.: Condensed Matter 16, 5279 (2004).
- (22) More precisely, we employ the so-called ‘approximation 1’ in the notation of Ref. frequenzen, .
- (23) A. Hewson, The Kondo Problem to Heavy Fermions (Cambridge University Press, Cambridge, 1993).