Driven Liouville-von Neumann Equation for Quantum Transport and Multiple-Probe Green's Functions
Francisco Ramírez, Daniel Dundas, Cristián G. Sánchez, Damian A. Scherlis, Tchavdar N. Todorov
TThis document is the unedited authors version of a Submitted Work that was subsequently acceptedfor publication in The Journal of Physical Chemistry C, copyright c (cid:13) ACS after peer review. To accessthe final edited and published work see: https://pubs.acs.org/doi/abs/10.1021/acs.jpcc.8b12319
Driven Liouville-von Neumann Equation forQuantum Transport and Multiple-Probe Green’s
Functions
Francisco Ram´ırez, ∗ , † Daniel Dundas, ‡ Cristi´an G. S´anchez, ¶ Damian A.Scherlis, † and Tchavdar N. Todorov ‡ † Departamento de Qu´ımica Inorg´anica, Anal´ıtica y Qu´ımica F´ısica/INQUIMAE, Facultadde Ciencias Exactas y Naturales, Universidad de Buenos Aires, Ciudad Universitaria,Buenos Aires (C1428EHA) Argentina ‡ Atomistic Simulation Centre, School of Mathematics and Physics, Queen’s UniversityBelfast, Belfast BT7 1NN, UK ¶ Departamento de Qu´ımica Te´orica y Computacional, Facultad de Ciencias Qu´ımicas,Universidad Nacional de C´ordoba, Ciudad Universitaria X5000HUA, C´ordoba Argentina,and CONICET & Facultad de Ciencias Exactas y Naturales, Universidad Nacional deCuyo, Mendoza, CP5500, Argentina
E-mail: [email protected] a r X i v : . [ c ond - m a t . m e s - h a ll ] M a y bstract The so called Driven Liouville-von Neumann equation is a dynamical formulationto simulate a voltage bias across a molecular system and to model a time-dependentcurrent in a grand-canonical framework. This approach introduces a damping term inthe equation of motion that drives the charge to a reference, out of equilibrium density.Originally proposed by Horsfield and co-workers, further work on this scheme has ledto different coexisting versions of this equation. On the other hand, the multiple-probescheme devised by Todorov and collaborators, known as the hairy-probes method, isa formal treatment based on Green’s functions that allows to fix the electrochemicalpotentials in two regions of an open quantum system. In this article, the equationsof motion of the hairy probes formalism are rewritten to show that, under certainconditions, they can assume the same algebraic structure as the Driven Liouville-vonNeumann equation in the form proposed by Morzan et al. [J. Chem. Phys. ,146, 044110]. In this way, a new formal ground is provided for the latter, identifyingthe origin of every term. The performance of the different methods are explored us-ing tight-binding time-dependent simulations in three trial structures, designated asballistic, disordered, and resonant models. In the context of first-principles Hamiltoni-ans the Driven Liouville-von Neumann approach is of special interest, because it doesnot require the calculation of Green’s functions. Hence, the effects of replacing thereference density based on the Green’s function by one obtained from an applied fieldare investigated, to gain a deeper understanding of the limitations and the range ofapplicability of the Driven Liouville-von Neumann equation.
The interest in molecular conductance and electronic transport across nanostructures hasinspired the development of a multiplicity of theoretical treatments to compute the currentunder an applied bias. The proposed schemes vary in complexity and computational cost,going from the original Landauer-B¨uttiker method and different static models meant to2escribe steady state transport, to dynamical methodologies that take into account thetime evolution of the charge density.
The present article is concerned with the laterclass of approaches, and in particular with the subset of methods based on the so-calledDriven Liouville-von Neumann (DLvN) equation. In this framework, the dynamics at theelectrodes is modulated by an additional driving term that augments the standard Liouville-von Neumann equation of motion. The role of these terms is to enforce part of the densitymatrix associated with the electrodes to remain close to a reference density matrix, thusintroducing a charge imbalance between a “source” or left lead ( L ) and a “drain” or rightlead ( R ). This driving term originally takes the form of a difference between the time-dependent density matrix ˆ ρ ( t ), and a reference density matrix ˆ ρ ,˙ˆ ρ ( t ) = − i (cid:126) [ ˆ H, ˆ ρ ( t )] − Γ( ˆ ρ ( t ) − ˆ ρ ) (1)where ˆ H is the electron Hamiltonian, Γ is the driving rate parameter, and the matrix ˆ ρ canbe defined as follows: ρ ij = ρ ij ( t ) if i, j ∈ L ∪ Rρ ij ( t ) if i, j / ∈ L ∪ R (2)This type of approach was first introduced by Horsfield and co-workers as an intuitiveway of including at the level of the density matrix the circulation of charge between theelectrodes, and it was further enriched by the work of Nitzan and Mazziotti. These ideaswere later taken on by Hod and co-workers and reelaborated by deriving the drivingterms from the theory of complex absorbing potentials. In doing so, they arrived at amodified form in which coherences were damped to zero,˙ˆ ρ ( t ) = − i (cid:126) [ ˆ H, ˆ ρ ( t )] − Γ2
2( ˆ ρ LL − ˆ ρ LL ) ˆ ρ LC ρ LR ˆ ρ CL ρ CR ρ RL ˆ ρ RC
2( ˆ ρ RR − ˆ ρ RR ) , (3)3hat improved the stability of the calculations and the steady-state convergence, and thatwas shown to satisfy Pauli’s principle regardless of the initial conditions. Franco and col-laborators presented a formal derivation of this formula from the theory of non-equilibriumGreen’s functions, demonstrating that it can accurately capture time-dependent transportphenomena. More recently, Zelovich et al. proposed a strategy to replace the single rateparameter Γ by diagonal matrices containing the broadening factors corresponding to thelead states. These factors can be computed from the self-energies of the electrodes, thusproviding a parameter-free version of the Driven Liouville-von Neumann approach. Because of its conceptual simplicity and good compromise between computational costand physical accuracy, the DLvN method has attracted significant attention, and severalfurther refinements of its implementation and analysis of its theoretical foundation weremade. Among these, of particular relevance is the adaptation to a first-principles real-timeTDDFT framework carried out by Morzan and co-authors, where an observed imbalancebetween injection and absorption of charge during the dynamics prompted a reformulationof the driving term. In an orthonormal basis, this equation of motion assumes the followingstructure: ˙ˆ ρ ( t ) = − i (cid:126) [ ˆ H, ˆ ρ ( t )] − Γ2
2( ˆ ρ LL − ˆ ρ LL ) ˆ ρ LC − ˆ ρ LC
2( ˆ ρ LR − ˆ ρ LR )ˆ ρ CL − ˆ ρ CL ρ CR − ˆ ρ CR
2( ˆ ρ RL − ˆ ρ RL ) ˆ ρ RC − ˆ ρ RC
2( ˆ ρ RR − ˆ ρ RR ) (4)where the subscripts indicate the corresponding blocks of the time-dependent density matrix,and ˆ ρ denotes a reference, time-independent density matrix. A major change with respect toprevious expressions is that here all off-diagonal contributions are damped to their referencevalues (obtained e.g. from the polarized density matrix). This proved to be important forcharge conservation and an appropriate balance between incoming and outgoing currents inthe relatively small models tractable in TDDFT simulations. Despite intuitively justifyingthis modification both from a numerical standpoint and interpreting it as a change in the4oundary conditions, no formal theoretical derivation was provided at the time.In parallel with these developments, a different driving term for the Liouville-von Neu-mann equation was put forward by Todorov and co-workers. This new embedding methodcalled multiple-probe or “hairy probes” (HP), formally derived a different structure of thedriving term via the application of Green’s functions and the Lippmann-Schwinger equa-tion.
In spite of its good results and wide range of applicability, its relationship withthe reference density driving method remained undetermined.The objective of the present study is to bridge the gap between the two methodologies:the heuristic formulation of the Driven Liouville-von Neumann equation, and the hairy-probes scheme. For that purpose, we have been able to rewrite the driving terms of the HPtheory into two terms, one of which resembles the damping term of the reference densityapproach. In doing so we have found that the latest addition of Morzan and co-workers isthe most compatible with HP, thus finally providing a formal framework for their correctionof the non-diagonal blocks of the driving term. Through the application of these schemes toa series of model systems, we better characterize these methodologies for different situationsand shed light on the reliability of the equation of motion including the driving term withthe reference density matrix ˆ ρ .In the next two sections we introduce the theoretical and methodological frameworkfor this work: we first show how the HP theory can be rewritten to arrive to the DrivenLiouville-von Neumann equation plus an additional term (section 2), and then provide a briefdescription of the models employed in the simulations (section 3). In section 4 we examinethe impact that the progressive simplifications of the HP equation have on the accuracy ofthe physical description of these systems. After that we examine the effect that some ofthe relevant parameters of the models have on the baseline performance of these differentapproximations. In particular, we will focus on the Γ parameter (section 5) and the shape ofthe field used to generate the reference matrix in the DLvN scheme (section 6). Finally, insection 7, we discuss the differences between the two forms of the DLvN equation reported5n the literature: the one emerging from the truncation of the HP method (eq. 4), and theone proposed by Hod and co-workers (eq. 3). We consider in this work a system formed by a central molecule or device, identified hereafterwith the symbol C , coupled to a right and a left electrode (or lead), denoted as regions R and L respectively. In the multiple-probes framework this system is embedded in an implicitlyrepresented environment via the coupling of each atom of the left and right leads to anexternal probe, which in isolation has a retarded and an advanced surface Green’s function g ± j ( E ) and a surface local density of states d j ( E ) = − π − (cid:61) g + j ( E ). These probes have fixedelectrochemical potentials µ L and µ R , depending on whether they are connected to the leftor right electrode, with Fermi-Dirac distributions f L ( E ) and f R ( E ). For this model, andassuming that the Hamiltonian of the system is time-independent, the HP theory providesthe following equation of motion for the electronic density (see eq. 38 in reference ): i (cid:126) ˙ˆ ρ S ( t ) = [ ˆ H S , ˆ ρ S ( t )] + ˆΣ + ˆ ρ S ( t ) − ˆ ρ S ( t ) ˆΣ − + (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] dE (5)where ˆ ρ is the time-dependent density matrix, ˆ H is the electron Hamiltonian, subscript S refers to the full system ( L , C and R sections), and the ˆΣ and ˆ G matrices are the system’s self-energy and Green’s function. We also assume that (i) all probes interact with the electrodesthrough the same coupling term γ , and (ii) the wideband limit applies to the probes so thattheir density of states becomes a constant, d j ( E ) = d . In these conditions, these matricesadopt the following form:ˆΣ ± = ˆΣ ± L + ˆΣ ± R = − i Γ2 · ˆ P L − i Γ2 · ˆ P R = (cid:16) ˆΣ ∓ (cid:17) † , (6)6Σ < ( E ) = Γ2 π · f L ( E ) · ˆ P L + Γ2 π · f R ( E ) · ˆ P R , (7)ˆ G ± S ( E ) = (cid:16) E · ˆ P S − ˆ H S − ˆΣ ± (cid:17) − . (8)In the above equations, Γ = 2 πγ d , and f L ( E ) and f R ( E ) correspond to the Fermi-Diracdistributions for the left and right probes. The matrices ˆ P L and ˆ P R are the correspondingprojector operators, with ˆ P S being the projection over the whole explicit system (i.e. theidentity).We will now take equation 5 as the starting point for our derivation. Our goal in thissection is to express it in terms of a difference between weighted density matrices—onecorresponding to the current state of the system and the other to a reference system—thatassumes the form of the driving term in the DLvN equations.For this, we focus on the second ( ˆΣ + ˆ ρ S ( t ) − ˆ ρ S ( t ) ˆΣ − ) and third ( (cid:82) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] dE ) terms on the right hand side of equation 5, and examine separately eachof the sub-blocks corresponding to the orbitals of the electrodes ( L , R ) and of the device( C ). The first thing to note is that, since all ˆΣ matrices contain projections on the leadsonly, the second and third terms will vanish in the central block: i (cid:126) ˙ˆ ρ CC ( t ) = (cid:16) [ ˆ H S , ˆ ρ S ( t )] (cid:17) CC . (9)In turn, the off-diagonal blocks of the second term can be expanded as:ˆΣ + ˆ ρ S − ˆ ρ S ˆΣ − = − i Γ2 (cid:16) ˆ P L ˆ ρ S + ˆ P R ˆ ρ S + ˆ ρ S ˆ P L + ˆ ρ S ˆ P R (cid:17) = − i Γ2 (2 ˆ ρ LL + 2 ˆ ρ LR + 2 ˆ ρ RL + 2 ˆ ρ RR + ˆ ρ LC + ˆ ρ CL + ˆ ρ RC + ˆ ρ CR ) (10)where for conciseness we have omitted the time-dependence of ˆ ρ . In matrix form, this can7e written as: ˆΣ + ˆ ρ S − ˆ ρ S ˆΣ − = − i Γ2 ρ LL ˆ ρ LC ρ LR ˆ ρ CL ρ CR ρ RL ˆ ρ RC ρ RR . (11)We now turn our attention to the third term on the right hand side of equation 5,considering separately each off-diagonal block. For the upper left one, we have: (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LL dE = Γ2 π (cid:90) (cid:16) f L ( E ) ˆ P L ˆ G − S ( E ) ˆ P L − f L ( E ) ˆ P L ˆ G + S ( E ) ˆ P L (cid:17) dE = Γ2 π ˆ P L · (cid:90) f L ( E ) (cid:104) ˆ G − S ( E ) − ˆ G + S ( E ) (cid:105) dE · ˆ P L = i Γ ˆ P L · (cid:90) f L ( E ) ˆ D S ( E ) dE · ˆ P L (12)where we used the relation between the Green’s function and the density of states matrixˆ D S ( E ) = (2 πi ) − (cid:16) ˆ G − S ( E ) − ˆ G + S ( E ) (cid:17) . The integral in the last term of equation 12 definesa fictitious equilibrium density matrix with an electronic distribution f L corresponding tothe left probes, that we will denote ˆ ρ L . After applying on the left and on the right theprojection operator associated with the upper left block, we arrive at the final result, (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LL dE = i Γ · ˆ P L · ˆ ρ L · ˆ P L = i Γ · ˆ ρ LLL . (13)For the RR sub-matrix we get an analogous expression, but with the R projectors andthe Fermi-Dirac distribution for the right probes, which defines a different reference densitymatrix ˆ ρ R , 8 [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] RR dE = i Γ ˆ P R · (cid:90) f R ( E ) ˆ D S ( E ) dE · ˆ P R = i Γ · ˆ ρ RRR . (14)To evaluate the off-diagonal terms involving the central region (blocks LC, CL, RC , and CR ), we first observe that in all these cases part of the integrand vanishes when operatingon ˆΣ < due to the contiguous application of projectors belonging to different regions, (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LC dE = (cid:90) (cid:16) ˆ P L ˆΣ < ( E ) ˆ G − S ( E ) ˆ P C − ˆ P L ˆ G + S ( E ) ˆΣ < ( E ) ˆ P C (cid:17) dE = Γ2 π (cid:90) f L ( E ) ˆ P L ˆ G − S ( E ) ˆ P C dE = Γ2 π ˆ P L · (cid:90) f L ( E ) ˆ G − S ( E ) dE · ˆ P C . (15)Similarly, (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] CL dE = − Γ2 π ˆ P C · (cid:90) f L ( E ) ˆ G + S ( E ) dE · ˆ P L . (16)To recover the density of states from just one of the Green’s matrices, we rewrite the latterin the following way:ˆ G − S ( E ) = (cid:32) ˆ G − S ( E )2 − ˆ G + S ( E )2 (cid:33) + (cid:32) ˆ G − S ( E )2 + ˆ G + S ( E )2 (cid:33) = iπ ˆ D S ( E ) + ˆ G (cid:60) S ( E ) (17) − ˆ G + S ( E ) = (cid:32) ˆ G − S ( E )2 − ˆ G + S ( E )2 (cid:33) + (cid:32) − ˆ G − S ( E )2 − ˆ G + S ( E )2 (cid:33) = iπ ˆ D S ( E ) − ˆ G (cid:60) S ( E ) , (18)for which we have defined the matrix ˆ G (cid:60) S ( E ) = ( ˆ G − S ( E ) + ˆ G + S ( E )). Inserting equations 179nd 18 in equations 15 and 16, we arrive at (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LC dE = i Γ2 ˆ P L · (cid:90) f L ( E ) ˆ D S ( E ) dE · ˆ P C + Γ2 π ˆ P L · (cid:90) f L ( E ) ˆ G (cid:60) S ( E ) dE · ˆ P C = i Γ2 ˆ ρ LLC + Γ2 π ˆΩ LLC . (19)Similarly (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] CL dE = i Γ2 ˆ ρ LCL − Γ2 π ˆΩ LCL (20)where we have introduced the matrix ˆΩ
L/R = (cid:82) f L/R ( E ) ˆ G (cid:60) S ( E ) dE . The terms involving theright lead can be treated on the same footing, to obtain an analogous result but incorporatingthe reference density and ˆΩ corresponding to the right region, (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] RC dE = i Γ2 ˆ ρ RRC + Γ2 π ˆΩ RRC (21) (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] CR dE = i Γ2 ˆ ρ RCR − Γ2 π ˆΩ RCR . (22)To complete the derivation, we consider the two remaining blocks that couple both leadstogether ( LR , RL ): (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LR dE = Γ2 π (cid:90) (cid:16) f L ( E ) ˆ P L ˆ G − S ( E ) ˆ P R − f R ( E ) ˆ P L ˆ G + S ( E ) ˆ P R (cid:17) dE = Γ2 π ˆ P L · (cid:90) (cid:16) f L ( E ) ˆ G − S ( E ) − f R ( E ) ˆ G + S ( E ) (cid:17) dE · ˆ P R . (23)The expression within the brackets can be rewritten introducing relations 17 and 18 to10licit the reference density and omega matrices, (cid:90) (cid:16) f L ( E ) ˆ G − S ( E ) − f R ( E ) ˆ G + S ( E ) (cid:17) dE = (cid:90) (cid:16) iπf L ( E ) ˆ D S ( E ) + f L ( E ) ˆ G (cid:60) S ( E ) + iπf R ( E ) ˆ D S − f R ( E ) ˆ G (cid:60) S (cid:17) dE = iπ ˆ ρ LS + ˆΩ LS + iπ ˆ ρ RS − ˆΩ RS , (24)leading to the following result: (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] LR dE = i Γ2 (cid:0) ˆ ρ LLR + ˆ ρ RLR (cid:1) + Γ2 π (cid:16) ˆΩ LLR − ˆΩ RLR (cid:17) (25)which combines the left and right probes populated matrices. The RL block is treatedanalogously, to obtain: (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] RL dE = i Γ2 (cid:0) ˆ ρ LRL + ˆ ρ RRL (cid:1) + Γ2 π (cid:16) ˆΩ RRL − ˆΩ LRL (cid:17) . (26)Thus, the third term in equation 5 can in matrix form be written as: (cid:90) [ ˆΣ < ( E ) ˆ G − S ( E ) − ˆ G + S ( E ) ˆΣ < ( E )] dE = i Γ2 ρ LLL ˆ ρ LLC ˆ ρ LLR + ˆ ρ RLR ˆ ρ LCL ρ RCR ˆ ρ LRL + ˆ ρ RRL ˆ ρ RRC ρ RRR + Γ2 π LLC ˆΩ LLR − ˆΩ RLR − ˆΩ LCL − ˆΩ RCR ˆΩ RRL − ˆΩ LRL ˆΩ RRC (27)where the first matrix on the right hand side will be referred to as the “region-weightedreference density matrix”, or simply reference density matrix.11inally, collecting all the pieces together, the HP master equation can be rewritten inthe following way:˙ˆ ρ ( t ) = − i (cid:126) [ ˆ H S , ˆ ρ S ( t )] − Γ2 (cid:126)
2( ˆ ρ LL − ˆ ρ LLL ) ˆ ρ LC − ˆ ρ LLC (cid:16) ˆ ρ LR − ˆ ρ LLR +ˆ ρ RLR (cid:17) ˆ ρ CL − ˆ ρ LCL ρ CR − ˆ ρ RCR (cid:16) ˆ ρ RL − ˆ ρ LRL +ˆ ρ RRL (cid:17) ˆ ρ RC − ˆ ρ RRC
2( ˆ ρ RR − ˆ ρ RRR ) − i Γ2 π (cid:126) LLC ˆΩ LLR − ˆΩ RLR − ˆΩ LCL − ˆΩ RCR ˆΩ RRL − ˆΩ LRL ˆΩ RRC (28)We have thus arrived at a mathematical structure that is very similar to that of theDLvN scheme as presented in equation 4. By this we mean that the second term has theform of a block-weighted difference between the current time dependent density matrix anda reference density matrix. The implications of the definition of this reference matrix andthe impact of including or excluding the remaining Ω-term are analyzed in the remainder ofthis work. In what follows we assess and compare the performance of three different implementations ofthe HP equation of motion: (i) the full formula 28; (ii) a version excluding the Ω-term; and(iii) same as (ii), but where the region-weighted reference density matrix is computed usinga step-shaped field encoding the bias potential. In the latter, notice that ˆ ρ L and ˆ ρ R are thesame matrix ˆ ρ , corresponding to the polarized density in the presence of a step-field, andthus the equation of motion becomes identical to equation 4. These three schemes (i), (ii)and (iii), will be referred to as “full-HP”, “partial-HP”, and “step-potential” methods, or, for12hort, F-HP, P-HP, and ST-P, respectively. As mentioned above, the mathematical structureof the P-HP approach resembles closely the one of the DLvN equation 4, implemented inthis study as the ST-P method. In particular, the driving term arising in the HP formulareinstates the damping of the coherences to the equilibrium values. In this sense, the ST-Pscheme can also be thought of as an approximation to the P-HP method with an alternativereference density matrix. The motivation to explore it is to bypass the calculation of Green’sfunctions, which becomes especially appealing in the context of ab-initio Hamiltonians.For comparison purposes, we also consider in this work the alternative form of DLvNgiven by equation 3, which hereafter will be referred to as DLvN-z, whereas the version inequation 4 will be denoted as DLvN-e (since these two equations damp the lead-moleculecoherences to zero and to their equilibrium values, respectively). It is important to notefrom this is that “ST-P” and “DLvN-e” refer to the exact same method: we will in generaluse the first notation when comparing only with other HP-derived methods, and the secondone when also including the other version of the DLvN equation.To examine the performance of these different equations of motion, a series of modelsystems were chosen as case studies in numerical simulations. All these systems consist oflinear chains of atoms arranged in different spatial configurations (see below). The electronicstructure is represented by a tight-binding Hamiltonian with an orthonormal one-electronbasis, using the model presented by Sutton and co-workers (other possible, more elaborateelectronic models are discussed in ). This scheme adopts a single orbital | j (cid:105) per atom (orsite), coupled with each other by nearest-neighbour hopping terms that only depend on thedistance between sites. The onsite energies will in general be identical for all atoms, exceptin the “disordered” model (see below), or in the presence of an external field used to obtainthe reference matrix.For the HP method, the potential difference ∆ V between the left and right probes de-termines the Fermi-Dirac distributions according to the corresponding electrochemical po-tentials µ L = +∆ V / µ R = − ∆ V /
2. The presence of an external field (to generate the13eference density to be used either in the modified HP or the DLvN schemes) is simulatedthrough the modification of the onsite energies. The shape of this field is a matter of analy-sis in section 6, but for the simulations in sections 4 and 5, constant values of +∆
V / − ∆ V / C ). In all systems the electrodes L and R are made of 200 atoms,except for some results shown in section 7 in which shorter leads of 50 atoms were explored.In all cases the systems are metallic, exhibiting a half-filled conduction band. These threesystems are: • Ballistic: all atoms in the system ( S ) are equally spaced by 2.5 ˚ A , involving a hoppingintegral of -3.88 eV (which corresponds to the tight-binding reference parameters forgold ). The on-site energy is the same for all atoms. There are 21 atoms in the centralregion C . • Disordered: this system has the exact same geometry as the ballistic (same number ofatoms all equally spaced), but with a different on-site energy for every atom belongingto the central region C . These different on-site energies were randomly generated usinga homogeneous distribution with values between plus and minus the absolute value ofthe hopping term above. • Resonant: all sites are equally spaced as in the previous cases, with the exception of14wo distances separating a group of 15 atoms in the center of the device, from twogroups of 9 atoms on each end, connected to the electrodes. All of these belong to the C region, as can be seen in Figure 1. These two bonds were 0.5 ˚A longer than the rest.The on-site energies are the same for all atoms, but the hopping terms of the longerbonds are -1.88 eV.Figure 1: Scheme of the resonant system. All on-site energies are the same.Modified versions of the resonant system were later used in section 6. Instead of having a9-15-9 configuration of the central region (with the dash standing for the longer separation),these are 9-45-9 and 30-15-30. Figure 2 depicts the current as a function of time for transport simulations based on theF-HP, P-HP, ST-P (or DLvN-e), and DLvN-z methods. It is seen that the currents reacha stable steady state, regardless of the method, model system, and the value of Γ. Theupper, middle and lower panels of Figure 2, correspond to the ballistic, the resonant, andthe disordered models, respectively, for an applied bias of 1.5 V (the same stable behavioris found for other bias). Left and right panels show the currents for Γ equal to 0.1 eV and0.6 eV. It can be seen that for large couplings, the steady state is reached faster, and thatin none of the cases is there any ambiguity in the identification of the final current, since ittypically converges to a precise value within the first 10 - 40 fs of dynamics, depending onΓ. While the figures present the current for the initial 400 fs, most of the simulations havebeen evolved for up to 1 ps, without the detection of instabilities.15igure 2: Current as a function of time obtained for a bias of 1.5 V with the three differentimplementations of the hairy-probes method (including the so called DLvN-e, given by equa-tion 4), and with the Driven Liouville-von Neumann formula given in equation 3 (DLvN-z).The upper, middle and lower panels, correspond to the ballistic, the disordered, and theresonant models, respectively. Left and right panels show the currents for Γ equal to 0.1 eVand 0.6 eV. 16igure 3 presents the steady state current-voltage curves obtained from the three methodsderiving from HP. For the ballistic system all three approaches show only marginal differencesbetween each other, which become negligible in the low bias region. The same is true forthe disordered case. In both kinds of model system, the differences in the steady-statecurrents resulting from the three methods are below 5%, and typically much less than thisat low biases. Somehow curiously, the scheme based on the step-potential reproduces moreaccurately the response of the full HP method, despite the difference in their referencedensities. This suggests the existence of a compensation mechanism through which the step-generated reference density matrix balances the absence of the ˆΩ contribution. A possibleexplanation can be found in its higher polarization, as discussed below.The I − V curves in the resonant system display a more complex behavior. The P-HPand ST-P equations of motion still reproduce the currents of the full HP method, but withslightly larger deviations. Here, these appear to be comparable in the P-HP and in theST-P method. The latter scheme produces in all three systems, for a given bias, highercurrents than the first one. This trend was confirmed also in exploratory calculations forother resonant systems with variable separation between atoms: the step-potential yieldssystematically larger steady states currents than the P-HP method for a fixed potentialdifference. This might be related to the fact that in the case of a reference density generatedfrom an abrupt step-field, the device is subject to a larger effective bias than in the caseof a smoother reference density constructed from the Green’s functions. In other words,the action of the field, through the direct modification of the on-site energies, affects moreseverely the density matrix of the electrodes when compared to the HP method, in whichthe electrochemical potential at the leads is controlled via the coupling with external probesaccording to γ . This is consistent with the eigenstates populations plotted in Figure 4. Thefull and the P-HP simulations produce a manifold of partially occupied states, both methodsexhibiting essentially no differences from each other. When the reference density derivesfrom the step-wise potential, instead, the eigenvalues of the density matrix are 0 or 1, with a17igure 3: Current as a function of the applied bias, computed in the steady state fromquantum dynamics simulations based on the three different implementations of the hairy-probes method. Results are shown for the ballistic (top), disordered (center), and resonant(bottom) model systems. 18mall fraction of states associated with intermediate occupations. This circumstance reflectsa lower electronic temperature and a more drastic polarization in the ST-P approach arisingfrom a highly perturbed ˆ ρ , in comparison with the other two treatments. On the other hand,in the former there are a few states that violate the Pauli’s principle, bearing negative, orlarger than unity occupations. The results in Figure 4 correspond to the resonant system,but analogous behavior is found for the ballistic and resonant models. The behavior of theeigenvalues for the different schemes, including the one in equation 3, is further discussed inthe final section. Γ parameter The question of the proper value of Γ is an interesting problem that has been investigatedrecently, but has not been entirely settled. In the context of the DLvN approach, it has beenargued in the literature that the magnitude of the driving rate should be somewhere betweenthe energy spacing separating the lead levels, and the effective device-lead coupling. If Γis too high, the target density is enforced in the lead sections too tightly to allow for anappropriate response and relaxation at the interface. As a consequence, in this regime thecurrent is a decreasing function of Γ (see e.g. refs. , , ). Moreover, when Γ is too low,the relaxation process in the leads is not sufficiently fast to allow for the electronic densityin the leads to reach ˆ ρ . The net effect in this case can be assimilated to a situation in whichthe leads are disconnected or weakly connected to the reservoirs and the system approachesa microcanonical regime.On the other hand, in the context of hairy probes Γ depends on the system parameters d and γ , where d is the probe surface density of states and γ the matrix element coupling thesites at the electrodes with the external probes. While the value of d can be defined to fit thedensity of states of the device, there is no particular criterion to uniquely assign the value of γ . In principle, in the wideband limit it is expected to satisfy the second requirement listed19igure 4: Temporal evolution of the eigenvalues of the density matrix for the resonant modelsystem, computed using the full HP (top), P-HP (center), and ST-P (bottom) methods.20bove for Γ in DLvN: it must be larger than the mean energy level spacing in the leads, toensure an effective broadening of the electronic states and hence a metallic behavior of theelectrodes. In the present case this implies a lower bound for Γ of around 0.1 eV.Figure 5 illustrates the effect of Γ on the I − V plots in the case of the resonant system,also showing the derivatives of the current (conductance) in the insets. We chose this systembecause in its case the I − V curve displays a distinctive physical pattern, and the discrepan-cies between the three approaches are more significant than in any other model, both factsmaking improvements easier to identify. Indeed, these discrepancies tend to fade away as Γdecreases from 0.6 to 0.1 eV (within this range Γ remains larger than the mean energy levelspacing in the leads, of 0.078 eV). Interestingly, the insets reveal that conductances derivedfrom the ST-P methodology turn out to be more rugged than in the other methods. Thiscan be attributed to the presence of the self-energy in the Green’s functions, that, providedΓ is larger than the energy level spacing in the electrodes, screens the far ends of the systemin the case of the F-HP and P-HP methods. In the absence of the self-energy, finite sizeeffects manifest in the form of interference oscillations. As a matter of fact, for a small Γ thewiggles are visible in all three curves (see the inset on the right panel of Figure 5), becausethe screening effect dilutes as the Green’s function carries a dependence on this parameter inthe denominator. The insets also reveal that a large Γ value tends to damp the resonancesin the ST-P case, which are otherwise intact for Γ = 0.1 eV.The currents obtained for a bias of 1.5 V are depicted as a function of Γ in Figure 6,where it can be seen that the agreement between methods is progressive as the parametergets smaller. Additionally, this Figure shows that all approaches are relatively robust againstthe variation of Γ, in particular the full HP scheme, for which the change in the steady statecurrent is just 2% when this parameter is reduced by a factor of six.In the F-HP and P-HP dynamics, the effect of Γ arises not only from the multiplicationof the driving terms, but also because the reference matrices ˆ ρ and ˆΩ are all functionsof Γ through the self-energy and the Green’s function (see equations 6 and 8). Therefore,21igure 5: Steady state current as a function of the applied bias computed with the threedifferent implementations of the hairy-probes method, for two different values of the Γ pa-rameter: 0.6 eV (left panel) and 0.1 eV (right panel). The results correspond to the resonantsystem. For a better comparison, the derivatives are shown in the insets.Figure 6: Steady state current as a function of the Γ parameter, computed with the threedifferent implementations of the hairy-probes method. The shown results correspond to theresonant system with an applied bias of 1.5 V.22t becomes relevant to check how the structure of these matrices change when varying thevalue of Γ. Figures 7 and 8 show that the elements of these matrices exhibit a slight inversedependence on Γ, which, as mentioned above, limits the spatial range of the Green’s functionsand reference density matrices in the electrode region, in a way that has no analogue in ST-P.This spatial damping is visible on Figures 7 and 8.Figure 7: Color maps representing the structure of different blocks of the ˆΩ matrix fordifferent values of the Γ parameter, computed for the resonant system. Top panels: LCblocks. Bottom panels: LR blocks. Left panels: Γ=0.1 eV. Right panels: Γ=0.6 eV. Thecolors reflect the absolute values of the matrix elements, in a natural logarithmic scale. Thismatrix determines the difference between the F-HP and P-HP schemes.In particular, the magnitude of the ˆΩ contribution is responsible for the difference betweenthe F-HP and P-HP methods. All diagonal blocks of the ˆΩ matrix are identically zero, andhence only the LC and LR blocks are depicted in Figure 7 (the remaining blocks, CL, RL,RC and CR are similar). The off-diagonal elements of ˆΩ are of the same order of magnitudeas those of the reference matrix ˆ ρ . Yet, the impact of the former matrix on the dynamics is23igure 8: Color maps representing the LL block of the reference density matrix ˆ ρ in theF-HP and P-HP methods, for different values of the Γ parameter, in the case of the resonantsystem. Left panel: Γ=0.1 eV. Right panel: Γ=0.6 eV. The colors reflect the absolute valuesof the matrix elements, in a natural logarithmic scale.only secondary because it has zeros throughout its diagonal blocks, which is where ˆ ρ givesthe highest contribution. In spite of its seemingly small importance in numerical terms,though, the presence of the ˆΩ matrix appears to balance the effect of ˆ ρ in the current,furnishing the full HP method with a reduced sensitivity with respect to Γ.Figure 9 displays the difference between the reference density matrices in the P-HP andST-P methods. Given that in the latter this matrix is independent of Γ, Figure 9 just reflectsthe effect of this parameter on the hairy probes reference density. Interestingly, since thedecrease of Γ tends to relax the spatial damping of the P-HP reference density, which is inany case absent from the step-field generated ˆ ρ , both matrices become with such a decreasemore similar to each other. The agreement between these two approaches at low Γ is manifestin the behavior of the currents in Figures 5 and 6.As mentioned in the introduction, Zelovich and co-workers proposed a method for thecomputation of a set of broadening factors that are applied to the lead states. In this way, therate parameter Γ is replaced by diagonal matrices ˆΓ L and ˆΓ R , with dimensions given by thenumber of basis sets associated with the leads. The calculation of these matrices involves aself-consistent procedure to extract from the self-energy of the isolated lead, the broadeningfactors that afterwards must be assigned to the corresponding levels of the lead coupled to24igure 9: Color maps representing the LL block for the difference between the referencedensity matrices of the P-HP and ST-P methods, for the resonant system and differentvalues of the Γ parameter. Left panel: Γ=0.1 eV. Right panel: Γ=0.6 eV. The colors reflectthe absolute values of the matrix elements.the reservoir. Whereas this procedure is somehow involved, the ˆΓ matrices are transferable toany calculation using the same lead model, and the propagation of the dynamics representsa negligible additional computational cost. The authors rationalized the magnitude of theresulting broadening factors Γ i by invoking Fermi’s golden rule for a single lead level coupledto a reservoir, which gives a value determined essentially by the coupling matrix elementbetween lead and reservoir states. This is precisely the meaning of γ in the context of hairyprobes. In any case, this treatment was shown to be of importance in particular when theDOS of the leads is inhomogeneous in the vicinity of the Fermi energy. However, the effectof considering a single rate parameter instead of one per level was shown negligible in simpletight-binding models as the present one, where the DOS is uniform around the Fermi energy.As a matter of fact, the authors reported that the adoption of the maximal broadening valuecalculated for such systems using their procedure as the single driving rate, produces currenttraces and steady-state occupations almost indistinguishable from those obtained using theparameter-free DLvN method. 25 Generation of the reference density
The encouraging results obtained with the ST-P implementation raise the question of whetherit is possible to better reproduce the F-HP dynamics with a truncated equation of motion,through the optimization of the reference density. This pathway is considered in the presentsection, using the shape of the electric field as a tool to prepare ˆ ρ , although other strategiescould be adopted as for example the state representation of Hod and co-workers. We initially compare the results obtained with a step-like potential, with those corre-sponding to a ramp, a sigmoidal decay, and a double sigmoidal decay at the start and theend of the central region. These different patterns are illustrated in Figure 10. The analysisof the reference matrices generated in this way did not reveal any significant improvement inthe description of the full hairy probes density, with the exception of a very slight increasein similarity for those matrix elements lying in the proximity of the central region. Thisimprovement was most noticeable in the case of the sigmoid field (denoted as “SI-P”), ascan be visualized in Figure 11.As a matter of fact, the adoption of the sigmoidal field proved to be at least as good asthe step potential, when not better, to reproduce the I − V plots. Figure 12 presents thesecurves for three resonant systems exhibiting different morphologies. Aside from the originalstructure displayed in Figure 1, two other models were examined in which the number ofatoms in either the central or the lateral segments was extended. The value of Γ was fixed to0.6 eV, for which the disagreement between methods was most noticeable. For the standardresonant system the performances obtained from the step or the sigmoidal potentials arecomparable. The same is true for the alternative resonant model with a longer intermediateregion, for which no significant differences are found between the results yielded by eithermethod. However, in this case the description provided by these two approaches is manifestlyworse. In this sense, it is noteworthy that the P-HP scheme is still able to capture the mainfeatures of the full hairy probes curve. On the other hand, in the third model where thedevice bears longer lateral segments, the reference calculated with the sigmoidal potential26igure 10: Shapes of the different electric fields applied to generate the reference density. L,C, and R, represent the left electrode, central, and right electrode regions, respectively.Figure 11: Color maps representing the LL block for the difference between the referencedensity matrix of the full HP method with the one generated with a step field (left), or theone generated with a sigmoidal field (right). Data corresponding to the resonant system forΓ=0.6 eV. The colors reflect the absolute values of the matrix elements.27utperforms the one calculated with the step-like field. Equations 3 and 4 are two expressions of the DLvN approach, both of which had an heuristicorigin, and for which different formal derivation routes have been proposed. From a mathe-matical point of view, the difference between these two formulas is found in the off-diagonalelements of the driving term: while in equation 3 these damp the lead-device coherences tozero, in equation 4 coherences are pushed towards their equilibrium values. Their derivationsinvolve different approximations and assumptions, and it is of particular interest to comparethe paths that lead to one or the other. This is the goal of the present section, where theperformances of the two schemes are also confronted.It is possible to identify three assumptions or approximations specific to one or the otherformula, that explain their different mathematical structures:1. In equation 3, the explicit lead levels are in contact with an implicit fermionic reservoir,for which coupling the wideband limit is assumed. Equation 4 is in turn a truncationof the hairy probes formula, where the leads are coupled to a set of probes in whichthe wideband limit holds. The latter procedure broadens but preserves the electronicstructure of the (finite) leads adjoining the central region.2. To arrive to DLVN-z, it is assumed that the relaxation dynamics in the leads is inde-pendent of the presence of the device. This amounts to the zeroth-order approximationof the Green’s functions in which G adv and G ret commute with the leads subspace pro-jector Q . This step eventually leads to the disappearence of ρ eq in the last term ofequation 36 in the work by Franco et al. Interestingly, the same result would also beobtained e.g. from equations 15 or 16 of our manuscript, if ˆ P L or ˆ P R were assumed28igure 12: Steady state currents as a function of the applied potential, computed withdifferent implementations of the hairy-probes method with Γ = 0.6 eV. Results are shownfor the resonant system with different arrangements of the atoms in the central region: 9-15-9(top), 9-45-9 (center), and 30-15-30 (bottom).29o commute with the Green’s functions ˆ G ± S . In such a case, these contributions wouldvanish, suppressing ρ L and ρ R from the off-diagonal elements of the equation of mo-tion, and thus giving rise to a driving term analogous to DLvN-z. Thus, the dampingof the coherences either to zero or to equilibrium is related to this approximation inthe Green’s functions, which physically translates into an independency of the electrondynamics in the leads with respect to the device.3. Finally, the ˆΩ terms are dropped from the HP expression to arrive at DLvN-e. Theseterms arise from the Green’s functions keeping track of the forward time-propagationof carriers coming from the probes. Thus, they introduce a propagation sense to theinjected lead-device charge. Its suppression must imply a drop in the current, as iscertainly observed in the simulations.In ref. it has been reasoned that while DLvN-z represents a device between electronicreservoirs at equilibrium with well-defined chemical potentials, DLvN-e resembles the stateof a charged capacitor, where the target density represents the equilibrium state of the entirefinite junction model and not just the leads. This seems to be consistent with the origins ofeach of these schemes, that have now come to light. Figure 13 depicts the current-voltagecurves for the two DLvN approaches, together with the results of the HP method. Atsmall couplings it can be seen that the performance of any of these methods is essentiallyindistinguishable from each other. As the Γ parameter is increased, discrepancies start toemerge.DLvN-z has been shown to observe Pauli’s principle regardless of the initial conditions. It is interesting to note that, while the hairy probes scheme also fulfills Pauli’s principle atlong times (by construction), it does not necessarily obey it in the transient. This can be seenin Figure 14, which displays the occupations for F-HP, P-HP, DLvN-e, and DLvN-z in theresonant system. It must be recalled that in all these schemes, the dynamics is switched onsmoothly, to avoid sudden jumps which could lead to numerical discontinuities. Specifically,the driving term and the Ω term are introduced in the first part of the simulation using a30igure 13: Steady state current as a function of the applied bias computed with the twoforms of the DLvN equation and with the hairy-probes method. The top panel presents thedata corresponding to the ballistic (left) and disordered (right) systems, for Γ equal to 0.6eV. The bottom panel shows the results for the resonant system with Γ equal to 0.6 eV (left)and 0.1 eV (right). 31inear ramp. For F-HP, the positive and negative deviations from 1 and 0 respectively—which become smaller with a smoother ramp—disappear in the long term. When the Ωterm is suppressed from the F-HP scheme, the dynamics and in particular the occupationsare affected but the exclusion principle is still obeyed in the steady state. However, ifthe reference density matrix is replaced by the one obtained from an electric field to givethe DLvN-e approach, the violation to the exclusion principle persists even in the steadystate. This is not observed with the DLvN-z method, for which the exclusion principle issatisfied throughout the full dynamics, despite the adoption of the same ˆ ρ as in DLvN-e. Interestingly, this makes manifest that the observance of Pauli’s principle is determinedneither by Ω nor by the reference density, but by their combination.Figure 14: Temporal evolution of the eigenvalues of the density matrix for the resonantmodel system. The top row corresponds to the full HP method where the driving term wasappplied gradually using ramps of different durations: 5.0 fs, 12.5 fs and 25.0 fs (from left toright). The bottom row compares the behavior found with the other methods using a rampof 5.0 fs: the P-HP (left), the DLvN-e (center) and the DLvN-z (right) schemes. Part of thisdata has been already given in a different size scale in Figure 4.Finally, Figure 15 explores the effect of Γ and of the electrode size on the currents, using32he resonant model. The bottom left panel shows that the HP method, consistently withits physical origin, is particularly robust with respect to electrode size and Γ. In the case ofan electrode of 50 atoms, significant deviations in the current are observed only for Γ < , which is a consequence of the robustnessof the HP method from which it is descended. In this article, it was shown that the Green’s functions based multiple-probes—or hairyprobes—formalism, adopts a form equivalent to the heuristic Driven Liouville-von Neumannmethod as proposed in reference (equation 4), plus an additional term involving a matrix( ˆΩ) with null diagonal blocks. A distinctive feature of this form of the DLvN equation33igure 15: Effect of electrode size and of the Γ parameter on the current-voltage curvesobtained for the resonant model. Top panels present in black the curve obtained from F-HPwith 200 atoms in each lead and Γ=0.6 eV (labelled ref), in comparison with the curvescomputed using 50 atoms leads and different Γ values for the DLvN-e (left) and DLvN-z (right) methods. We also include the behaviour of the F-HP method with the smallelectrodes and varying Γ (bottom left). The bottom right panel displays, for DLvN-e andDLvN-z (Γ = 0.3 eV), the difference with respect to the F-HP current (obtained with Γ =0.6 eV). 34f motion is that, at variance with the previous versions introduced in references or , the coherences are damped to the equilibrium density. It has been argued that in thisapproach the electrodes are not meant to represent infinite reservoirs with homogeneousand well defined chemical potentials. Instead, through the action of the driving operator,the leads are driven close to, rather than exactly to, the target density. This is parallelto the physics in the HP model, where the electrodes are connected to multi-probes withelectrochemical potentials µ L/R . The electrochemical potential of the lead is not necessarilythat of the probes but depends on the strength of their coupling (and on position down thelead). These particular boundary conditions of the HP method are reminiscent of those inthe DLvN implementation of equation 4 and demonstrated to be well suited for small sizeelectrodes. To neglect the ˆΩ matrix in the driving term of the HP formula has minor effects onthe dynamics and the steady state currents obtained for a variety of model systems. Theseeffects are even less significant when the coupling between the probes and the leads is reducedby decreasing the Γ parameter. More specifically, our results show that the P-HP methodcan reproduce the behavior described by the full HP scheme for ballistic and disorderedsystems, and with some care it can also be tuned for more complicated resonant systems. Ingeneral, at least in the context of tight-binding Hamiltonians, it is possible to conclude thatthe DLvN method, incarnated here in the P-HP or ST-P equations of motion, convergesto the hairy probes description in the limit of small couplings between the probes and theleads (providing Γ is still larger than the energy level spacing). The P-HP method canbe thought of as a version of the DLvN approach in which the calculation of the referencedensity is based on Green’s functions. The ST-P method, on the other hand, reproduces theformulation presented in reference . It has proved to be a very good approximation to thefull HP description, but the finite size effects which manifest in the absence of a self-energyimply a limitation in comparison with the P-HP approach that seems difficult to overcome.The possibility to avoid the calculation of Green’s functions acquires special interest for35he applications in the context of TDDFT or other first-principles schemes. In that respect,we explored the substitution of the reference density computed from Green’s functions inthe P-HP method, by the one emerging from an electric field applied to the system. Thisstrategy systematically produced higher currents than the P-HP method, presumably be-cause the application of a field, in the present setting, amounts to an additive constant inthe density matrix elements, whereas in the HP scheme the chemical potential is fixed at anexternal probe and not directly on the leads, whose polarization is mediated by a couplingparameter. The effective bias arising between the electrodes in operation conditions maythus not be as large as the one imposed by the field. Further sources of improvement for theregion-weighted field-generated reference density method proved hard to find. Our resultshere seem to suggest that the step and sigmoideal shapes fields can better fit those parts ofthe reference matrices at the boundary with the device, providing the most accurate repre-sentations of the current in comparison with the other fields tested. Whereas the ballisticand disordered models could be described reasonably well with the ST-P and SI-P schemes,the resonant systems proved to be challenging for those methods using electric-field basedreference densities. This limitation becomes particularly relevant for TDDFT applications,where the complexity of the chemical structures, far from the simple tight-binding modelsexamined in the present work, may lead to stronger discrepancies between the results ob-tained using a reference density generated from Green’s functions or from an electric field.It would be interesting to establish how the different reference densities and the omission ofthe ˆΩ matrix affect the dynamics in the case of first-principles simulations. That questionwill be the subject of future work. Acknowledgement
This work was supported by a research grant from Science Foundation Ireland (SFI) and theDepartment for the Economy Northern Ireland under the SFI-DfE Investigators Programme36artnership, Grant Number 15/IA/3160, by the Argentinean Agency for Scientific and Tech-nological Promotion (ANPCYT) through PICT 2015-2761, and by the University of BuenosAires, UBACYT 20020160100124BA. We are grateful to Uriel Morzan for very appreciateddiscussions.
References (1) Landauer, R. Spatial Variation of Currents and Fields Due to Localized Scatterers inMetallic Conduction.
IBM J. Res. Dev. , , 223–231.(2) Buttiker, M. Four-Terminal Phase-Coherent Conductance. Phys. Rev. Lett. , ,1761.(3) Brandbyge, M.; Mozos, J.-L.; Ordej´on, P.; Taylor, J.; Stokbro, K. Density-FunctionalMethod for Nonequilibrium Electron Transport. Phys. Rev. B , , 165401.(4) Kosov, D. S. Lagrange Multiplier Based Transport Theory for Quantum Wires. J.Chem. Phys. , , 7165–7168.(5) Goyer, F.; Ernzerhof, M.; Zhuang, M. Source and Sink Potentials for the Descriptionof Open Systems with a Stationary Current Passing Through. J. Chem. Phys. , , 144104.(6) Koentopp, M.; Chang, C.; Burke, K.; Car, R. Density Functional Calculations ofNanoscale Conductance. J. Phys-Condens. Mat. , , 083203.(7) Thoss, M.; Evers, F. Perspective: Theory of Quantum Transport in Molecular Junc-tions. J. Chem. Phys. , , 030901.(8) Kurth, S.; Stefanucci, G.; Almbladh, C.-O.; Rubio, A.; Gross, E. K. U. Time-DependentQuantum Transport: A Practical Scheme Using Density Functional Theory. Phys. Rev.B , , 035308. 379) Burke, K.; Car, R.; Gebauer, R. Density Functional Theory of the Electrical Conduc-tivity of Molecular Devices. Phys. Rev. Lett. , , 146803.(10) McEniry, E. J.; Bowler, D. R.; Dundas, D.; Horsfield, A. P.; S´anchez, C. G.;Todorov, T. N. Dynamical Simulation of Inelastic Quantum Transport. J. Phys-Condens. Mat. , , 196201.(11) M¨uhlbacher, L.; Rabani, E. Real-Time Path Integral Approach to NonequilibriumMany-Body Quantum Systems. Phys. Rev. Lett. , , 176403.(12) Horsfield, A. P.; Bowler, D.; Fisher, A. Open-boundary Ehrenfest Molecular Dynamics:Towards a Model of Current Induced Heating in Nanowires. J. Phys-Condens. Mat. , , L65.(13) S´anchez, C. G.; Stamenova, M.; Sanvito, S.; Bowler, D.; Horsfield, A. P.; Todorov, T. N.Molecular Conduction: Do Time-Dependent Simulations Tell You More than the Lan-dauer Approach? J. Chem. Phys. , , 214708.(14) Subotnik, J. E.; Hansen, T.; Ratner, M. A.; Nitzan, A. Nonequilibrium Steady StateTransport Via the Reduced Density Matrix Operator. J. Chem. Phys. , ,144105.(15) Rothman, A. E.; Mazziotti, D. A. Nonequilibrium, Steady-State Electron Trans-port With N-Representable Density Matrices from the Anti-Hermitian ContractedSchr¨odinger Equation. J. Chem. Phys. , , 104112.(16) Zelovich, T.; Kronik, L.; Hod, O. State Representation Approach for Atomistic Time-Dependent Transport Calculations in Molecular Junctions. J. Chem. Theory Comput. , , 2927–2941.(17) Zelovich, T.; Kronik, L.; Hod, O. Molecule-Lead Coupling at Molecular Junctions:38elation Between the Real- and State-Space Perspectives. J. Chem. Theory Comput. , , 4861–4869.(18) Hod, O.; Rodrguez-Rosario, C. A.; Zelovich, T.; Frauenheim, T. Driven Liouville vonNeumann Equation in Lindblad Form. J. Phys. Chem. A. , , 3278–3285.(19) Zelovich, T.; Kronik, L.; Hod, O. Driven Liouville von Neumann Approach for Time-Dependent Electronic Transport Calculations in a Nonorthogonal Basis-Set Represen-tation. J. Phys. Chem. C. , , 15052–15062.(20) Chen, L.; Hansen, T.; Franco, I. Simple and Accurate Method for Time-DependentTransport along Nanoscale Junctions. J. Phys. Chem. C , , 20009–20017.(21) Zelovich, T.; Hansen, T.; Liu, Z.-F.; Neaton, J. B.; Kronik, L.; Hod, O. Parameter-FreeDriven Liouville-von Neumann Approach for Time-Dependent Electronic TransportSimulations in Open Quantum Systems. J. Chem. Phys. , , 092331.(22) Morzan, U. N.; Ram´ırez, F. F.; Gonz´alez Lebrero, M. C.; Scherlis, D. A. ElectronTransport in Real Time from First-Principles. J. Chem. Phys. , , 044110.(23) Todorov, T. N. Tight-Binding Simulation of Current-Carrying Nanostructures. J. Phys-Condens. Mat. , , 3049.(24) Dundas, D.; McEniry, E. J.; Todorov, T. N. Current-Driven Atomic Waterwheels. Nat.Nanotechnol. , , 99.(25) McEniry, E.; Wang, Y.; Dundas, D.; Todorov, T.; Stella, L.; Miranda, R.; Fisher, A.;Horsfield, A.; Race, C.; Mason, D. et al. Modelling Non-Adiabatic Processes UsingCorrelated Electron-Ion Dynamics. Eur. Phys. J. B , , 305–329.(26) Horsfield, A. P.; Boleininger, M.; D’Agosta, R.; Iyer, V.; Thong, A.; Todorov, T. N.;White, C. Efficient Simulations with Electronic Open Boundaries. Phys. Rev. B , , 075118. 3927) Sutton, A. P.; Todorov, T. N.; Cawkwell, M. J.; Hoekstra, J. A Simple Model of AtomicInteractions in Noble Metals Based Explicitly on Electronic Structure. Philos. Mag. A ,81