Merging Features from Green's Functions and Time Dependent Density Functional Theory: A Route to the Description of Correlated Materials out of Equilibrium?

We propose a description of nonequilibrium systems via a simple protocol that combines exchange-correlation potentials from density functional theory with self-energies of many-body perturbation theory. The approach, aimed to avoid double counting of interactions, is tested against exact results in Hubbard-type systems, with respect to interaction strength, perturbation speed and inhomogeneity, and system dimensionality and size. In many regimes, we find significant improvement over adiabatic time dependent density functional theory or second Born nonequilibrium Green's function approximations. We briefly discuss the reasons for the residual discrepancies, and directions for future work.

We propose a description of nonequilibrium systems via a simple protocol that combines exchangecorrelation potentials from density functional theory with self-energies of many-body perturbation theory. The approach, aimed to avoid double counting of interactions, is tested against exact results in Hubbardtype systems, with respect to interaction strength, perturbation speed and inhomogeneity, and system dimensionality and size. In many regimes, we find significant improvement over adiabatic time dependent density functional theory or second Born nonequilibrium Green's function approximations. We briefly discuss the reasons for the residual discrepancies, and directions for future work. DOI: 10.1103/PhysRevLett.116.236402 Hybrid methods are a valuable option in physics, to merge concepts and perspectives into a more general and effective level of description. This Letter adds an item from condensed matter physics to the list; we propose a hybrid method that combines nonperturbative exchangecorrelation (XC) potentials from time dependent density functional theory (TDDFT) [1][2][3] with many-body perturbative self-energy schemes from nonequilibrium Green's functions (NEGF) [4][5][6][7], to deal with systems with strong electronic correlations and out of equilibrium.
An accurate first-principles description of the real-time dynamics of systems with strong electronic correlations is an important, difficult, and basically unsolved problem of condensed matter research. General frameworks like TDDFT and NEGF do indeed allow for an in-principleexact treatment of strong electronic correlations. However, they both rely on key ingredients (the XC potential for TDDFT and the self-energy Σ for the NEGF) that in general are only approximately known. For TDDFT, a systematic and general way to include nonlocal, nonadiabatic effects in the XC potential is lacking, while for NEGF a main hindrance is that self-energies based on many-body perturbation theory, already computationally demanding, are usually inadequate for strong electronic correlations. While considerable progress has been made for model systems far away from equilibrium (see, e.g., [8][9][10][11][12][13][14]) or for the ab initio description of nearequilibrium situations (see, e.g., [15,16]), a reliable firstprinciples treatment of the far-from-equilibrium regime is still lacking.
Here, we suggest a step towards the solution of this problem, by a novel combination of TDDFT and NEGF, where perturbative (but systematic) memory-effect corrections augment a nonperturbative local adiabatic treatment of electronic correlations. The approach is fully conserving in the Kadanoff-Baym sense [17] and, using the so-called generalized Kadanoff-Baym ansatz [18] (see below), can be made viable for realistic systems.
Putting in practice our proposal at the ab initio level requires access to continuum nonperturbative XC potentials, and this point is addressed at the end of the paper. However, the scope of our method can already be illustrated here using simple lattice models. This has the advantage of avoiding complex implementations and technicalities that, indispensable to deal with real-world systems, are usually unnecessary (possibly even unwanted) for an explorative assessment of a new methodology. Our results show that in many situations (see also Supplemental Material [19]) the hybrid method provides significant progress over adiabatic-TDDFT and perturbative schemes for NEGF, thus holding promise for an improved treatment of the nonequilibrium dynamics of realistic correlated systems.
Model systems.-We consider small Hubbard-type onedimensional (1D) and three-dimensional (3D) clusters, isolated or coupled to two 1D semi-infinite noninteracting leads. In the latter case, the cluster consists of one site (single impurity). These systems are exposed to time dependent local perturbations and/or (where applicable) to electric biases in the leads. The Hamiltonian for the above setups isĤ which has contributions from the cluster, the leads, and the cluster-leads couplings. In standard notation, where hiji labels nearest-neighbor sites in the cluster C, V 0 > 0 is the tunneling amplitude, ϵ i ðtÞ are time dependent on-site energies in the cluster, and U i are contact-interaction strengths. Further,n i ¼n i↑ þn i↓ . For the lead Hamiltonian,Ĥ l ¼ P αĤα , where α ¼ RðLÞ refers to the right (left) lead, and Here, b α ðtÞ is the (site-independent) bias in lead α, V > 0 is the tunneling amplitude, andN α ¼ P i∈αni . The coupling between the leads and the cluster (impurity) is given bŷ All energy units are expressed in terms of the hopping parameter V 0 (for the one-site impurity cluster we use V instead), and time is measured in the units of the inverse hopping parameter (assuming atomic units). We now switch to continuum variables for generality and notational convenience, and provide some elements of TDDFT and NEGF relevant to our approach.
is the singleparticle Hamiltonian, with kinetic energy T, Hartree potential v H , and external potential v ext . Σ ¼ Σ emb þ Σ XC ½G is the self-energy, which introduces a memory dependence. We integrate over the Keldysh contour γ [4,5]. Σ emb is an embedding self-energy that accounts for the leads (if present), while Σ XC accounts for XC effects [20]. Standard approximations for Σ XC are the second Born (2BA), T matrix (TMA), and screened interaction (GW) [6,7,21]. For real time, the lesser part of G (denoted G < ) gives the density nðt; rÞ ¼ −iG < ðt; r; t; rÞ and the current.
TDDFT.-The time dependent density n KS is obtained in terms of the Kohn-Sham (KS) orbitals ϕ κ ðt; rÞ. These obey and v XC accounts for XC effects. Then, n KS ðt; rÞ ¼ P occ κ jϕ κ ðt; rÞj 2 . Within a NEGF treatment, the KS density can be obtained from In practical implementations, the functional dependence of v XC on n is often replaced by an adiabatic local density approximation (ALDA), i.e., v XC ð½n; t; rÞ ≈ v ref XC (nðt; rÞ). A hybrid TDDFT-NEGF approach.-Our proposal is to augment a perturbative self-energy Σ PT XC from a conserving many-body scheme with a nonperturbative XC potential v NP XC , local in space and time. Alternatively, this prescription can be seen as recasting an ALDA-TDDFT approach based on v NP XC in a NEGF approach, but augmenting it with a nonlocal, nonadiabatic perturbative self-energy Σ PT XC . To avoid double counting we subtract an ALDA potential v PT XC obtained from the same approximation as was used for Σ PT XC . The basic equation of our approach is To actually proceed with Eq. (5), at t ¼ 0 we solve for G; i.e., we find v NP XC ½G, v PT XC ½G, and Σ PT XC ½G self-consistently on the imaginary-time track; then we propagate G selfconsistently on the Keldysh contour, thus fulfilling the conservations laws of Kadanoff and Baym [17]. The hybrid scheme involves no additional computational costs compared to standard NEGF time propagation. Since the augmentation v NP XC ðtÞ − v PT XC ðtÞ is of the form of a timelocal potential, our scheme can similarly be implemented in a density matrix formalism. This means that a generalized Kadanoff-Baym ansatz (GKBA) [18,22,23] can be employed to reduce computational costs allowing for first-principles calculations of realistic systems.
The v PT XC correction.-For concreteness, in this Letter Σ PT XC and v PT XC are computed in the 2BA (some results in the TMA are also shown). The calculation and use of Σ 2B XC for Hubbard-type interactions in a NEGF time evolution has been discussed before (see, e.g., [30]) and is not repeated here. Rather, we provide additional details of the perturbative correction v 2B XC . For the homogeneous (Hubbard) reference system, we use v 2B , and the three terms on the rhs are respectively the total energy in the 2BA, the noninteracting kinetic energy, and the Hartree energy for the 1D (or 3D) homogeneous Hubbard model. We compute E 2B tot ðnÞ in ðω; qÞ space, with G R being the retarded propagator, f the statistical Fermi factor (we consider zero temperature), ϵ q the single-particle energies, and n ¼ ð−2=ð2πÞ Dþ1 Þ R ∞ −∞ dω × R BZ dqImG R ðω; qÞfðωÞ. In Fig. 1 we plot v 2B XC for the 1D and 3D Hubbard model, for different interaction values. We also show the nonperturbative potentials v BALDA XC , v DMFT XC used in Eq. (2). They exhibit a discontinuity at half filling, which is always present in 1D but only for large U values in 3D, reflecting the Mott-Hubbard metal-insulator transition [25]. The discontinuity is absent in the 2BA. Note that, at exactly half filling, v NP XC and v PT XC are both zero. Finally, v PT XC from the TMA is shown. The discontinuity is absent also in this case, and at low and high filling v TMA XC approaches v NP XC . Closed systems: the 3D case.-We start our analysis with a 3D cubic cluster with 5 3 sites, open boundary conditions, and a single interacting (and perturbed) site at the cluster center [ Fig. 2(c)]. We compare time dependent densities from the hybrid approach, 2BA, and ALDA against exact results. The system is highly inhomogeneous, and despite the local character of the interaction and external perturbation, nonlocal effects are important: the exact v XC (not shown) can have large nonzero components at all sites [25]. Using symmetry, we map the cluster to a 10-site one [ Fig. 2(d)], as in [25]. We consider both weak (U ¼ 8, panels a and b) and strong (U ¼ 24, panels c and d) correlations. The temporal shape of the external fields we use is Gaussian (bottom-row panels, red curves), with a slower or faster onset or offset (in the following referred to as fast or slow perturbations). For additional time profiles we refer the reader to Supplemental Material [19].
For the weakly correlated, slowly perturbed case (panel a), all approximations follow the exact solution. For the fast perturbation (panel b), nonadiabatic effects emerge, and this leads to the failure of the ALDA; the remaining approximations perform well, with the hybrid method being marginally better than the 2BA. In contrast, for the slow perturbation and stronger correlations (panel c), the agreement of the 2BA is poor, while the other treatments still follow the exact solution. For the most unfavorable and extreme regime of strong correlations and fast perturbations (panel d), the ALDA and 2BA are largely out of phase, and only the hybrid approximation reproduces the main structures of the exact solution with the correct phase. Overall, the hybrid approximation exhibits a fairly good agreement in all regimes, and is superior to the others in the most extreme regime.
Closed systems: the 1D case.-We next consider when all sites are interacting and exposed to a space and time dependent perturbation. A 3D system for this situation that is also an exactly solvable benchmark is not easily accessible, due to the unfavorable scaling of the configuration space. We thus turn to a numerically more convenient 1D test case (this also makes it possible to assess the hybrid approach at low dimensionality), choosing a 1D ring with eight interacting sites (Fig. 3). To explore the role of space inhomogeneity, we resort to a (rather artificial) perturbation sinusoidally modulated in space: V ext ðl; tÞ ¼ sin½ð2π= and F is the temporal profile. The phase ϕ k guarantees that the sine nodes are between sites and the amplitude at site l ¼ 3 has always the same sign. For the time profile, FðtÞ ≡ θðtÞ (step, s), FðtÞ ≡ 1=ð1 þ e −t Þ (ramp, r) or FðtÞ ≡ expð−t 2 Þ (Gaussian, g). Results are shown in Fig. 3 (for a more systematic study see Supplemental Material [19]).
With highly inhomogeneous fields (λ 1 , λ 2 ) no approximation reproduces the exact dynamics. Moreover, for r k¼1 the hybrid method shows artificial density oscillations. The latter, also present in the TDDFT-ALDA approach based on v BALDA XC , are induced by the sharp discontinuity in v BALDA XC and are not removed by the 2BA self-energy (thus, nonlocal, nonadiabatic effects beyond the 2BA should also be taken into account). For more homogeneous fields (λ 3 ) the different approximations compare more favorably to the exact dynamics with superiority of the hybrid method. Looking at s k¼3 , the hybrid approximation is in phase with the exact curve but, for densities changing across half filling, it still exhibits the artificial oscillations (see Supplemental Material [19]). Further, the ALDA does not perform well, and the 2BA tends to be out of phase with the exact solution. Finally, for a slowly varying-inspace perturbation (λ 4 ) the hybrid approach (in contrast to the other approximations) is in excellent agreement with exact results. This applies for all time profiles g, s, r.
Open systems.-Finally, we test the hybrid method in open systems (Fig. 4). Specifically, using a single-orbital Anderson impurity coupled to two 1D semi-infinite leads [31] [system shown in Fig. 4(e)], we consider (i) the conductance G in the wide-band limit (WBL), Fig. 4(a); and (ii) the finite-bias, finite-lead-width regime, Figs. 4(b)-4(e). Starting with (i), we find the exact density (and thereby the exact linear conductance via the Friedel sum rule) in the WBL [32][33][34][35]. Figure 4(a) displays for U ¼ 2 the absolute deviation Δ from the exact G as a function of V gate and for different approximate treatments. We consider stronger correlations (Γ WBL ¼ V 2 link =V ¼ 0.09; see the plateau in the conductance); here, except for 0.15 < n=2 < 0.28, the hybrid method performs as the best compared to the 2BA or ALDA, and it is significantly better in the range 0.28 < n=2 < 0.42 (symmetrical considerations apply above half filling). (ii) Next, we consider 1D tight-binding leads (of bandwidth 4V). We fix a static V gate to be away from the particle-hole symmetric ground state (where v NP XC ¼ v PT XC ¼ 0). As a benchmark, we use open-ended, Anderson-impurity finite chains with up to L ¼ 96 sites treated with time dependent density matrix renormalization group (tDMRG) [36,37]. When V link ¼ 0.5 (panels b and d), the agreement between hybrid and tDMRG densities or currents is fairly good, especially in the transients (n and j from tDMRG never fully reach a steady state within the simulation time, in contrast to the hybrid, 2BA, and ALDA ones). However, for stronger correlations and lower transparency U=V link ¼ 2=0.3 (panels c and e), the impurity density from the hybrid scheme is closest to the tDMRG one than other schemes, while for the currents the ALDA performs best. The unconvincing performance of the hybrid approximation for U=V link ¼ 2=0.3 comes probably from multiplescattering processes, neglected by the 2BA. To corroborate this conjecture we have tested the hybrid method also using the TMA, which includes such processes. In Figs. 4(c) and 4(e) the TMA hybrid method shows an improvement over the ALDA and the pure TMA calculation and thus supports the conjecture (for an expanded discussion and additional results, see the Supplemental Material [19]).
Conclusions and outlook.-By merging elements of TDDFT and NEGF, we proposed a simple, easy to implement, nonequilibrium scheme aimed to improve the treatment of local nonperturbative correlation effects and, at the same time, incorporate nonlocal, nonadiabatic effects. Results from Hubbard-type systems are quite encouraging. Taking a mildly optimistic stand, we can argue that our approach extends the applicability of the ALDA-TDDFT approach and NEGF based on perturbation theory, thus providing a way forward to merge (strong) correlations and memory effects in general. On the other hand, one can certainly envisage situations where nonperturbative and nonlocal correlations are very important, and this is where perhaps corrections beyond the 2BA (e.g., GWA or TMA, mixed, or other) could be employed. We note that Hubbardtype systems usually are challenging benchmarks to perturbative approximations such as the 2BA, TMA or GWA. The latter generally perform much better for continuum systems with long-range interactions. Thus, we speculatively suggest that our hybrid method could perform even better for realistic systems. This is where the real merits of our proposal could possibly be: Using continuum XC potentials tailored for strong correlations (obtained from, FIG. 4. Single-impurity, one-orbital Anderson model with U ¼ 2 (shown in panel e). a): Linear conductance G in the wide-band-limit for Γ WBL ¼ 0.09 (strong correlations). The exact G is displayed, together with the Hartree-Fock (HF), 2BA, ALDA and hybrid-method results. The density/spin-channel n=2 at the impurity is also shown (dashed line). n=2 and G share the same vertical scale (in different units). b-e): Time dependent density n (b, c) and average current ¼ ðj L þ j R Þ=2 (d, e) for the Anderson impurity with constant impurity gate voltage V gate ¼ ϵ 0 ¼ 0.25 and bias b L ðtÞ ¼ 0.5θðtÞ. The hopping parameter in the leads is V ¼ 1, the impurity-lead coupling is V link ¼ 0.5 (b,d) and V link ¼ 0.3 (c,e).