Nonadiabatic Electron Dynamics in Tunneling Junctions: Lattice Exchange-Correlation Potential

The search for exchange-correlation functionals going beyond the adiabatic approximation has always been a challenging task for time-dependent density-functional theory. Starting from known results and using symmetry properties, we put forward a nonadiabatic exchange-correlation functional for lattice models describing a generic transport setup. We show that this functional reduces to known results for a single quantum dot connected to one or two reservoirs and furthermore yields the adiabatic local-density approximation in the static limit. Finally, we analyze the features of the exchange-correlation potential and the physics it describes in a linear chain connected to two reservoirs where the transport is induced by a bias voltage applied to the reservoirs. We find that the Coulomb blockade is correctly described for a half-filled chain, while additional effects arise as the doping of the chain changes.


INTRODUCTION
The miniaturization of electrical circuits is of great importance for future technologies. Fundamental building blocks for circuits can implement their functionality within single molecules, e.g., by trapping them in mechanically controllable break junctions, and transport can be studied experimentally. 1−9 At these length scales, the understanding of such devices requires a quantum mechanical description. A wealth of approaches can be used to study said devices, and through the years, a lot of effort has been put into this matter. One of the most successful theories for describing charge transport through nanoscale devices is the Landauer−Buẗtiker (LB) approach, 10−12 where transport is expressed as a scattering problem. In the LB framework, the electrons traverse the device ballistically; i.e., electron−electron scattering is ignored. However, the LB approach can be generalized for interacting electrons within nonequilibrium Green's (NEGF) function theory. 13,14 A numerically more efficient alternative to NEGF for the description of interacting quantum systems is density-functional theory (DFT). 15,16 Since transport is by definition a nonequilibrium phenomenon, timedependent DFT (TD-DFT) 17−19 has to be employed for studying molecular transport. 20 This leads to the remarkable result that, in principle, the steady-state flow of interacting electrons though a transport device can be computed from the LB formula, provided it is evaluated using non-interacting Kohn−Sham (KS) electrons of TD-DFT moving in an effective local potential. In practice, however, this approach is limited by the quality of the available approximations to the effective potential. A common approach is to use approximations derived in the context of equilibrium or ground state DFT in order to evaluate the effective potential in the device region. This procedure is referred to as adiabatic approximation within TD-DFT, which ignores dynamical effects (memory). 21−23 Embedding approaches, like dynamical mean-field theory 24 or density matrix embedding theory, 25 have received a lot of attention recently due to their ability to describe strongly correlated electrons. In the context of density-functional theories, impurity models, such as the Anderson impurity model, 26 have been studied extensively over the past years in the Kondo and Coulomb blockade regime. 27−39 These studies provided important insight in how to describe strongly correlated physics within DFT approaches by reverse engineering the (nearly) exact effective potentials for simple model systems. A big remaining challenge is to generalize these effective potentials for generic transport setups.
In this work, we put forward an effective potential based on its analytical derivation for the case of a single impurity coupled to one lead. 38 A recent study comparing results using this effective potential to results from numerically exact, time-dependent density-matrix renormalization group demonstrated the quality of this effective potential for an impurity site coupled to a single lead. 40 We generalize the effective potential to any tight-binding transport setup. The generalization to any system is carried out by exploiting a simple gauge symmetry argument: in a system composed of an impurity connected to a lead, the application of a potential to the impurity is equivalent to the application of the opposite potential to the lead. We first show how this argument can be used to generalize the effective potential of the one site− one lead setup 38 to a simple transport setup composed of one site coupled to two leads. The resulting effective potential yields the known exact results, which reproduce Coulomb blockade binding chains, giving the dispersion ϵ αk = −2t cos k. This means that the electrons have the same hopping amplitude in the leads and in the device region. In the following, we consider the spin unpolarized case; i.e., observables do not depend on spin, e.g., n i↑ = n i↓ ≡ n i /2. We use a TD-DFT approach to numerically describe the dynamics of the system. This means that the evolution of the charge density for the system governed by Hamiltonian (1) is computed from an effectively noninteracting system. This is governed by a Hamiltonian similar to eq 1, where the interaction term is replaced by an effective local potential, the so-called KS potential. Once an approximation for the effective potential is given, the charge flow can be computed from the non-interacting system. In the steady state (SS), this charge flow is completely determined by the spin-summed density matrix of the device region, D SS , the currents between sites i and j of the device, I ij SS , and the currents from the device to lead α, I α

SS
. For the non-interacting system, they are given by Here, f(ε + b α ) is the Fermi function of lead α with associated voltage bias b α , and the matrices D α and T αα′ are defined as D α is the density of states of the device projected on lead α, and the off-diagonal matrix elements of T αα′ (ε) are the transmission coefficients between leads. Both are defined in terms of the advanced/retarded Green's function of the device region, G A/R , and the decay rate associated with lead α, Γ α . The trace is performed over the device region. In principle, Γ α is energy dependent, however, here we take the wide-band limit approximation. This disregards the energy dependence of the leads' density of states, taking it constant and equal to its value at the Fermi energy. 41−47 This is justified when the applied bias is small compared to the leads' bandwidth. One major advantage of taking the wide-band limit is that it allows for an analytical evaluation of the LB formula at finite temperatures. 48−50 3. GENERALIZED NONADIABATIC EXCHANGE-CORRELATION FUNCTIONAL The effective potential contains two contributions: the scattering induced by the nuclei or the boundaries of the device and the effect of electron−electron collisions. The former is commonly referred to as "external potential" and characterizes the device, while the latter is the so-called (Hartree-)exchangecorrelation (Hxc) potential. The Hxc potential derived in ref 38 for a quantum dot coupled to a single reservoir is given by Ä where Γ is the decay rate associated with the only reservoir, U is the interaction strength, and β is the inverse temperature. We point out that there is a typo in the corresponding formula (equation 8b in ref 38) that has been corrected here. Note the dependence of the functional on the normalized current c: this will be important in the generalization to generic lattice geometries. The potential (4) can be rewritten as a functional of the density n and current to the dot I by using

Journal of Chemical Theory and Computation
Article the continuity equation I ≡ −n. The two forms of the functional given in eq 4 are equivalent but numerically stable in different regimes: the first one is stable for n − I/(2Γ) > 1, while the second one is stable for n − I/(2Γ) < 1. These expressions differ from the one given in ref 38 by a constant of U/2. This constant shift ensures half-filling at vanishing chemical potential (μ = 0) for any interaction strength. Starting from this functional, we first present a generalization to a reservoir-dot-reservoir system, then to a dot connected to any number of leads, and finally to any transport setup in the tight-binding model. The functional (4) provides a form for the dot potential in the dot-reservoir system. By gauge symmetry, the dynamics are not altered if one applies the opposite potential to the reservoir instead. In general, any pair of dot (v Hxc ) and reservoir (V xc ) potentials v n c ( , ) Hxc Hxc with 0 ≤ α ≤ 1 gives rise to equivalent dynamics by gauge symmetry.
The crucial aspect of the previous argument is that we can subdivide the system into two parts. When we consider a transport scenario, where an additional reservoir is attached to the device, there are two possible ways of subdividing the system by grouping the dot either with the left or the right lead. For both configurations, we can apply the procedure leading to eq 5a, and summing the resulting potentials leads to v n c n c ( , ) where we use the fact that in the steady state the currents flowing from the dot to the two leads are opposite to each other by Kirchhoff's law and that both leads are characterized by the same decay rate. By imposing the additional condition α 1 + α 2 = 1, one can ensure that in the limit of vanishing current the potentials in the reservoirs are identically zero. This means that the density of the dot does not induce a potential anywhere else but locally when no currents flow (LDA). The generalization to any number of leads connected to a dot follows immediately from the previous considerations, and the potentials read where c k = I k /Γ k is the normalized current flowing to reservoir k. Again, the condition 1 k k α ∑ = leads to a local approximation in the adiabatic limit (no currents), and V xc k is identically zero for all k's. This condition is still not enough to fix an arbitrary number of α's. We make the choice α k = α motivated by the fact that only the dimensionless currents, c k , enter the functional. As a direct consequence, the first order nonadiabatic correction for the potential v Hxc , i.e., the correction due to having finite but small currents flowing, vanishes by virtue of Kirchhoff's law (charge conservation in the steady state) provided all decay rates are identical. With the aforementioned considerations, the potentials can be rewritten as where N L is the number of reservoirs connected to the dot. From the expressions (8a and 8b), it is evident how the electron density of the dot induces a potential locally and on the reservoirs, i.e., the elements connected to it. This is important for the generalization to any tight-binding model, where the dot is connected not only to leads but also to other sites. For each site, we have the local potential v Hxc and the nonlocal contributions V xc coming from other sites. When only nonlocal contributions due to directly connected sites are considered, the exchange-correlation potential of a general site i can be written as where the first term is the local contribution and the remainder is the sum of the nonlocal ones. The local component v Hxc i is eq 8a calculated with the density n i and the normalized currents c i k going out of site i. The nonlocal part, V xc ik ≡ v Hxc k − ϵ Hxc (n k ,c k i ), is the nonlocal contribution (8b) evaluated at the density of the nearest neighbors and the associated currents going to site i. Finally, ⟨k|i⟩ indicates the restriction of the sums to the nearest neighbors only. In equilibrium, when no currents are flowing in the system, the nonlocal contribution vanishes and the effective potential reduces to the LDA for lattice systems. For the reservoirs, the local component, v Hxc i , cannot be defined as no density can be associated with them in the wide-band limit and, like in the expression for the reservoir-dot-reservoir, only the nonlocal contributions appear. It is important here to notice the dependence of the functional on c i k : this is the normalized (dimensionless) current flowing from site i to site k. The normalization factor, Γ k , can only be evaluated for sites connected to leads. For currents between sites in the device region, eq 3c can be evaluated for V α → t, yielding 2t for the normalization factor. Equation 9 represents the main result of this paper. It is applicable to any lattice system of any dimensionality and connectivity with local interactions. In the following section, an analysis of the effect of this nonadiabatic functional is presented.

RESULTS
In this section, we investigate the influence of the derived exchange-correlation functional on the transport properties of nanojunctions. In particular, we analyze the steady state of tightbinding chains connected to two leads, one at each end of the chain. The development of a steady state is a consequence of the application of a potential difference to the leads. The xcpotential (9) is then found self-consistently as it depends on the steady state densities and currents and vice versa. To compare results coming from different functionals or different interaction strengths, we tune the gate voltage to have a fixed number of electrons in the chain. The reason for choosing this approach stems from the strong dependence of the number particles on the temperature of the leads (for a fixed gate), which is reflected, e.g., in the resistance, and this can be different for different functionals/interaction strengths.
The tight-binding parameters chosen for the simulations expressed in units of the hopping amplitude, which is set to t = 1 eV, are temperature 0 < k B T = β −1 ≤ t and bias voltage

Journal of Chemical Theory and Computation
Article b L,R = ±t/8. The other parameters that determine the junction are not constant for all calculations and will be specified when needed. All resistances shown in the following are normalized by the resistance (R 0 ) of a non-interacting infinite chain at zero temperature at the particle-hole symmetric point with hybridization Γ = 0.5.

Half-Filling.
When the gate voltage is set to 0, the junction is at the particle-hole symmetric point, meaning there are N electrons in it, with N being the length of the chain. In Figure 2, the resistance of a junction of length N = 80 at different interaction strengths is shown as a function of the temperature of the system. As expected, the non-interacting system is metallic; i.e., its resistance is a monotonically increasing function of the temperature and becomes linear for high temperatures. In this regime, the effect of the interaction is negligible as the energy scale associated with temperature is much larger than interaction strength. The role of the interactions becomes more important at lower temperatures, where it hinders transport, resulting in an increase of the resistance due to the Coulomb blockade. Furthermore, Figure 2 shows that at half-filling a larger interaction strength results in higher resistance. This effect is due to the potential induced in the leads, which suppresses the bias voltage, effectively shrinking the transport window, resulting in a smaller current flowing through the junction. The adiabatic functional, in turn, gives rise to only a minimal (negligible) increase in the resistance and does not describe the Coulomb blockade physics, as it misses this nonlocal potential induced in the leads. As stated before, the correct Coulomb blockade physics of the single impurity model is caught by the functional; 39 see the inset of Figure 2.
When considering lower temperatures, a peculiar feature appears in the resistance of the junction shown as a softening of the Coulomb blockade effect. In the upper panel of Figure 3, we report the resistance of junctions of length N = 80 with different interaction and hybridization strengths: U ∈ {0.3, 0.4, 0.5} and Γ ∈ {0.32, 0.5, 1.125}. The value of Γ determines the hightemperature behavior of the resistance, since the interaction becomes irrelevant when T ≫ U, as observed above. We observe, as expected, a higher resistance for lower values of Γ, since this is related to the hopping amplitude to the leads, facilitating tunneling through the junction at higher values. For fixed Γ, the interaction strength dictates the temperature at which the system transitions to the Coulomb blockade regime and also the position of the softening feature. In particular, the feature appears at higher temperatures when there is a higher interaction. Moreover, for fixed U, the softening appears amplified as the hybridization strength deviates from the hopping t, linking this peculiar feature to the dynamics of the edge regions. In fact, even though finite size and odd−even effects affect the resistance of the junction, its qualitative behavior persists up to convergence with the system size; see the lower panel of Figure 3. In all the previous calculations, the shown results do not go down to T = 0; they stop in fact at βU = 20, as for example in Figure 3. This is due to practical difficulties in the convergence of the self-consistent exchangecorrelation potential. Its determination is, in fact, a fixed point iteration with a function (ϵ Hxc ) highly nonlinear in the density and current. Moreover, the exchange-correlation potential is an oscillating function in a chain junction, stemming from the underlying Friedel's oscillation of the non-interacting case. These complications are worsened at lower temperatures, where the Fermi−Dirac distribution becomes step-like, sharpening all the nontrivial features of the functional. An undamped selfconsistent iteration, in fact, fails to converge as temperature goes down and cannot reach the regime where the Coulomb blockade is established. Advanced mixing schemes are needed to converge the functional at low temperatures: in this work, we use a Pulay mixing scheme. 51 At each iteration, this method exploits the history of guesses for the quantity needed to establish a better guess for minimizing the error in the next iteration. In the presented results, the potential has been mixed with a history length of 9. This is found to give faster and more stable convergence for this problem.
4.2. Doped Regime. By applying a gate voltage to the chain junction, we can tune the number of particles occupying the system and move it away from the particle-hole symmetric point. This affects the transport properties of the system itself. Different functionals and/or different interaction strengths yield a different number of particles in the junction for a fixed gate potential. For this reason, we choose to compare results by

Journal of Chemical Theory and Computation
Article tuning the gate voltage to have a fixed number of electrons in the chain for all temperatures. In the top panel of Figure 4, we compare the resistances of a chain of length N = 80 with a fixed total number of particles N el = 60 resulting from the adiabatic and nonadiabatic functionals at different interaction strengths. In accordance with ref 23, the nonadiabatic approximation yields a higher resistance. An additional distinction is given by the qualitative behavior of the resistances. While the adiabatic functional becomes less resistive as the interaction strength increases (behavior in agreement with mean-field results, not shown), the nonadiabatic approximation exhibits a crossover. In fact, the resistances for different interaction strengths cross at a temperature k B T/t ≃ 0.18. The decrease in the resistance for increasing interaction strength in the adiabatic functional is explained simply by the change induced in the transmission functions T αα′ due to the exchange-correlation potential. For the nonadiabatic case, this effect competes with the suppression of the transport window caused by the nonlocal potential induced in the leads. At high temperatures, the transmission functions are unaltered by the interaction; hence, the effect of the bias renormalization, although small in amplitude, dominates. At a fixed temperature, |V xc | is larger for larger interaction strengths; this explains why larger U's produce larger resistances. As temperature decreases, the two effects have comparable impacts and compete with one another. |V xc | starts decreasing roughly at a constant k B T/U, meaning different temperatures for different interaction strengths; see the lower panel of Figure 4. This results in a stronger drop in the resistance for, e.g., U = 0.70 while the transport window suppression for weaker interaction is still growing, resulting in the crossing of the resistances. This effect is amplified as we move away from the particle-hole symmetric point, since by doing so the magnitude of the bias renormalization drops. At a fixed number of particles N el = 76, in fact, the resistances do not cross and the results are similar to the ones at particle-hole symmetry.

CONCLUSIONS AND OUTLOOK
In this work, we have derived a novel nonadiabatic, nonlocal exchange-correlation potential for studying transport in devices modeled by lattice Hamiltonians within TD-DFT. Our exchange-correlation potential is based on the solution for a single quantum dot coupled to a reservoir, which exhibits a steplike feature important for the description of Coulomb blockade. 38 The derivation is built by exploiting gauge symmetry: the dynamics of a reservoir-dot device stay unchanged if a potential is applied to the dot or to the reservoir with opposite sign. For a generic lattice model, we apply this idea to every link through which a current flows, independent of whether it connects two sites or a site to a lead. The proposed functional reduces to the LDA 27 in the static limit; i.e., nonadiabaticity and nonlocality of the potential are linked closely together in our construction. In the nonequilibrium situation, our approximation reduces to the exact result for a single dot transport setup. 39 We have analyzed the effect of this exchange-correlation potential in chain systems connected by the ends to two reservoirs. We found that the "step" in the approximation for the effective potential makes it challenging to find a self-consistent solution for low temperatures/strong interactions, since the step sharpens as βU increases. First results for two-dimensional tightbinding models indicate that this is amplified for lattices with higher coordination. To this end, it would be interesting to develop alternative strategies for converging self-consistent approaches based on new insights. 52 For one-dimensional tightbinding chains, we have shown how Coulomb blockade physics is well described by this functional as an effect of the suppression of the transport window due to an induced potential in the leads opposite to the voltage bias applied. At the particle-hole symmetric point, an additional feature appears at roughly constant βU, showing a softening of the Coulomb blockade effect. The amplitude of this softening depends on the decay rate of the leads Γ. Away from particle-hole symmetry, the nonadiabatic functional differs qualitatively from the LDA as the on-site attraction (repulsion) competes with the renormalization of the transport window produced by the nonlocality of the functional. We believe that TD-DFT for tight-binding models offers a numerically efficient scheme for studying transport in quantum systems approaching mesoscopic scales. The aim of this work was to propose an approximation for the effective potential taking electron−electron scattering effects into account, which are believed to play an important role in socalled superballistic transport and has recently received a lot of attention. 53−55 An interesting route to explore is to study thermoelectric transport using a novel DFT approach proposed recently, 50,56−59 by generalizing the exchange-correlation potential for describing the combined charge and energy flow through the transport device.

Notes
The authors declare no competing financial interest.