Density Functional Theory of the Hubbard-Holstein Model

We present a density functional theory (DFT) for lattice models with local electron-electron (e-e) and electron-phonon (e-ph) interactions. Exchange-correlation potentials are derived via dynamical mean field theory for the infinite-dimensional Bethe lattice, and analytically for an isolated Hubbard-Holstein site. These potentials exhibit discontinuities as a function of the density, which depend on the relative strength of the e-e and e-ph interactions. By comparing to exact benchmarks, we show that the DFT formalism gives a good description of the linear conductance and real-time dynamics.

Similarly to what is done in (TD)DFT for electronnuclei systems in the continuum [49,50], or quantum electrodynamics [29], we describe the HH model via a twocomponent formulation for the electron and phonon subsystems, where for the electron (phonon) component the basic variables are the electron occupations {n i } (phonon coordinates {x i }) at each site i, and each component is governed by its own XC potential.
After presenting our approach, we explicitly determine the XC potentials for the analytically solvable onesite, zero-dimensional (D = 0) model, and the infinite dimensional (D = ∞) homogenous Bethe lattice (where n i = n, x i = x), via dynamical mean field theory (DMFT) [51][52][53]. These potentials are then used to compute the dynamics in a finite system. We find that: i) e-ph interactions screen the e-e interaction, and the behavior of the electronic XC potential is mainly determined by the screened interaction U . At the same time, because the e-ph coupling is linear in both n and x, as well as local, the XC phonon potential is always zero. ii) The electronic XC potential is discontinuous at half-filling density n = 1 for U > 0, and at n = 0, 2 for U < 0; for the Bethe lattice, the discontinuity appears above a nonzero value of |U |. iii) For an infinite chain with a HH impurity, (TD)DFT conductances have a smooth transition from the charge-to the spin-Kondo regime upon varying the e-ph coupling. iv) TDDFT dynamics in a test system subject to interaction quenches or external fields compares well with exact numerics in an appreciable range of interaction strengths. These results demonstrate that (TD)DFT is a promising formalism for the study of e-ph lattice systems.
The system -Aiming for a DFT description, we start with an inhomogeneous version of the HH Hamiltonian: where v i is a local site dependent electron potential, µ is the chemical potential, U the e-e interaction strength, J the hopping amplitude (set equal to 1, as the energy unit), and . . . denotes nearest neighbour sites. The operator c † iσ creates an electron on site i with spin σ, witĥ n iσ = c † iσ c iσ the corresponding density operator. The phonon frequency is ω, and g is the e-ph coupling parameter. We use λ = g 2 /ω as a measure of the e-ph interaction strength. Finally, the site-dependent external phonon po-arXiv:1903.04984v1 [cond-mat.str-el] 12 Mar 2019 tential η i is introduced to control the phonon coordinateŝ where b i destroys a phonon at site i. The form ofĤ in Eq. (1) allows to address formal aspects of the (TD)DFT description (see the Supplemental Material, SM), and to use a homogeneous HH reference system (∀i, v i = v and η i = η) in (adiabatic) local density approximations. To calculate the ground state energy of the reference homogenous HH model, we perform the Lang-Firsov transformation H →Ĥ = e iŜĤ e −iŜ [54] (see the SM). In H , the hopping amplitude is renormalized as 2 the phonon momentum, and the other parameters transform as v → v = v + (g 2 + 2gη)/ω, and U → U = U − 2g 2 /ω.
Density functional theory -For a DFT description, we consider the pair of sets of variables (n, x) ≡ ({n i }, {x i }) and the conjugated fields (v, η) ≡ ({v i }, {η i }), with n i = n i↑ +n i↓ the total electron density at site i. In the SM we prove that i) the total energy E = E[n, x] is a functional of n and x, with a minimum at the ground state values (n 0 , x 0 ) (i.e., the Hohenberg-Kohn theorem for the HH) and ii) where v 0 representability holds, (n 0 , x 0 ) is obtained by solving a two-component Kohn-Sham (KS) problem. In a homogeneous system (useful to derive a local density approximation), n i = n and x i = x for all i, and the KS Hamiltonian is The electronic KS potential can be written as v KS = v ext + v Hxc , with v ext ≡ v and the Hartree-exchangecorrelation (Hxc) part v Hxc = U n/2 + √ 2gx + δE xc /δn. To obtain E xc , we subtract from E[n, x] the e-e and e-ph Hartree interaction terms, and the energies of the corresponding non-interacting HH system. For the phonons, η KS = η ext +η Hxc , where η ext ≡ η and η Hxc ≡ η H +η xc = √ 2g(n − 1) + δE xc /δx. Using the Heisenberg equation ∂ tp = 0, we get x = − √ 2 [g(n − 1) + η] /ω. Inverting this relation and using η Hxc = η| g=0 − η , we find η Hxc = η H , i.e. η xc = 0. In the following, we set ω = 1 and η = 0.
A) v xc from a single Hubbard-Holstein site -We consider a single HH site exchanging energy and particles with a bath at chemical potential µ and temperature β −1 . This problem is analytically solvable, by determining the partition function Z = Tr e −βH(v,η) (see SM).
where δn = n − 1, R = [e −βU 1 − δn 2 + δn 2 ] 1/2 , and U defined above. Thus, v xc is independent of x and depends on the e-e interaction renormalized by phonons. For g = 0 we recover the expression for a single-site Hubbard system [23]. In Fig. 1(a) we display v xc as a function of n and λ for U = 1 and β = 20. For U > 0 (i.e. λ < 1/2) the potential is discontinuous at n = 1, as in the case of the purely electronic repulsive Hubbard model [7,23,57,58]. For U < 0 (i.e. λ > 1/2) the discontinuity is at n = 0 and n = 2, as in a negative-U Hubbard model [20,48]. Notably, in the present model, the transition from positive to negative U results from the phonons screening of the e-e interactions.
Eq. (3) shows that, save for the linear term g 2 δn/ω, the analytic expressions for v xc in a HH and a Hubbard single-site model only differ by the renormalization U → U , i.e e-ph interactions primarily affect the discontinuities at n = {0, 1, 2}. Phonon effects are instead explicitly manifest in the electronic spectral function A( ). Starting from the many-body Matsubara Green's func- can be extracted via analytic continuation to real energies and the fluctuation-dissipation theorem (see SM). Fig. 1(b,c) shows A( ) for µ = 0, β = 5, v = −U/2 and for four pairs (U, λ). For U = 3 and λ = 0.08 the two main peaks correspond to the electronic excitation energies. Instead, for λ = 1, phonon replicas spaced by ω are seen. A similar behavior occurs at U = −1: for small λ, A( ) has two main peaks. Here, the electronremoval/addition parts of A contribute to both peaks, since the e-ph interaction reorders the energy levels.
The zero bias conductance G is related to the spectral function [55]. Using v xc from Eq. (3), we calculate G at zero temperature for a HH impurity connected to two 1D semi-infinite noninteracting leads. In this case, G/G 0 = sin 2 (πn/2) [56], with n = (2/π) arctan(−v KS [n]/γπ) + 1, G 0 the unit of quantum conductance, and γ the level width in the wide-band limit. Results for G as a function of v and λ are in Fig. 1(d), where U = 1 and µ = 0. For U > 0, G has a plateau of width U [23], but for U < 0 we find a single narrow peak [48,56]. Overall, G behaves smoothly as a function of the e-ph coupling, while the system evolves from the spin to the charge Kondo regime.
For further insight into e-e and e-ph interactions, in the SM we compare the exact double occupancy n ↑ n ↓ to the DFT one (the latter is obtained via the single-site potential). While good agreement is found in many situations, the comparison in the SM clearly suggests that kineticenergy effects in v xc are important and thus a more general reference system than a single HH site should be used. This is taken into account in the next section.
B) v xc from the infinite-dimensional Bethe lattice -For the HH model on the D = ∞ Bethe lattice with bandwidth 4 (in units of the hopping parameter), we estimate the ground state energy E tot within DMFT at β = 200 [53]. The exchange-correlation potential is explicitly determined for U = 3, 8 and g = 0.2, 1, corresponding in all cases to a screened interaction U > 0 [59]. After the calculation of E xc as a function of n and x via DMFT (see the SM), we compute the electron and phonon potential as v xc = ∂E xc (n, x)/∂n and η xc = ∂E xc (n, x)/∂x. To perform the necessary derivatives, we fit the DMFT data in n ∈ [0, 1] with piecewise fourth-order polynomials. For n ∈ [1, 2], we employ the symmetry E xc (n) = E xc (2 − n).
In Fig. 2(a)-(d) we show the v xc obtained from DMFT (for E xc results, see the SM). At U = 8, v xc is discontinu-ous at n = 1, but not at U = 3. This is a DFT signature of the Mott-Hubbard transition in the HH model, in analogy with the purely electronic Hubbard model [58] (for the D = ∞ Bethe lattice, when β → ∞, U c1 ≈ 4.7, and U c2 ≈ 5.8 [60]). Interestingly, e-ph interactions not only renormalize the value of the XC discontinuity, but also "delay" its onset. For further insight, in Fig. 2 we also plot by yellow lines the single-site results from Eq. (3), with the value of β fitted to best reproduce the DMFT curves. For panels (a) and (b), this gives β = 5, and a smeared XC discontinuity. Also, due to the small U and large g values, the shape of the single-site solution in panel (b) is dominated by the linear term g 2 ω δn. In contrast, in panels (c) and (d) the fit gives β = 100 for the single-site potential, already close to the zero temperature limit (where the discontinuity exists for all nonzero interactions U ). Overall, single-site and DMFT potentials agree for n ∈ [0, 1/2] but significant differences appear at higher fillings, with important consequences for time-dependent simulations.
Bethe lattice mapping and real-time dynamics -We now use the single-site and DMFT potentials for the realtime dynamics of an L-site chain with a HH impurity at one end (so-called Anderson-Holstein chain, see Fig. 3, left). This test system was chosen because the local density of states of a homogeneous Bethe lattice of coordination Z and hopping term J is identical to the one at the first site (site 0) of a semi-infinite chain. The mapping is obtained via Lanczos recursion (see SM), and also holds with HH interactions and time-varying fields at a single site of the Bethe lattice: in this case the local Green's function at that site is the same as the one at site 0 of the chain. When Z → ∞, the chain Hamiltonian is where t labels time, v(t) is a local perturbation, and the rescaling J → J/ √ Z keeps the hopping probability finite. In the simulations, we use a finite chain of L = 8 sites, which allows for exact numerical solutions. By virtue of the mapping, we are actually dealing with a Z = ∞ Bethe lattice truncated after eight layers and with one HH impurity in the center (Fig. 3). We consider N ↑ = N ↓ = 3 electrons in the chain (as before, J = ω = 1). The system's time evolution is performed via exact diagonalization, as well as by TDDFT time propagation via the KS equations [61] within the adiabatic local density approximation (ALDA) [62]. By setting v xc to zero, we also consider the Hartree-Fock (HF) dynamics. Figure 3(a)-(c) shows the dynamics after a sudden interaction quench (U i , g i ) → (U f , g f ) at t = 0. This situation is within the scope of TDDFT, by freedom of choice of the initial state [4]. Further, quenches severely test the ALDA (typically employed within TDDFT, and used here). For the quench (0, 0) → (3, 1), panel (a) in Fig. 3, exact and TDDFT-DMFT results are in excellent agreement, while the single-site and HF solutions give a moderately good description. Instead, for (8, 1) → (3, 1) and (3, 1) → (8, 1), panels (b) and (c), the agreement worsens, due to stronger interactions. However, the DMFT potential still qualitatively performs well, while the HF solution fails to capture the main features. There is excellent agreement between exact and TDDFT results. In (e, f) v(t) is a soft square pulse of duration T = 4 and amplitude v = 1, switched on and off in a time T = 8. In (e), where U = 8 and g = 1, the DMFT potential initially gives a very good agreement but this worsens near n = 1. This is a known behavior, due to the discontinuity of v xc at n = 1 [26,58]. Finally, we briefly turn to the U < 0 region, using the single-site potential (the case of the DMFT potential is left for future work). The results of panel (f), where U = −1 and g = 0.5, suggest that (TD)DFT can also be used for the attractive regime, but better potentials than the single-site one are clearly needed.
Conclusions.-By means of a two-component density functional theory (DFT), we have introduced a novel approach to the Hubbard-Holstein (HH) model, a popular template to study e-e and e-ph interactions in lattice systems. We also explicitly determined and characterized electron and phonon exchange-correlation (XC) po-tentials, analytically for a simple one-site system, and via dynamical mean field theory for a Bethe lattice in D = ∞. Comparisons between DFT and exact results showed that the newly found potentials perform very well across an appreciable range of interaction strengths and electron densities.
Possible directions for immediate extensions include the analysis of the phonon overscreening regime, and a formulation for lattices in D ≤ 3 or for the linear response regime. Application-wise, an appealing option would be to explore how phonon-like degrees of freedom affect the physics of cold atoms in optical lattices (e.g. cloud expansion after trap removal, disorder effects, entanglement distillation, etc.). Finally, a key development would be the introduction of memory and non local effects in the XC potentials, by exploiting connections to many-body approximations within Green's function schemes.