A parameter-free version of the recently developed driven Liouville-von Neumann equation [T. Zelovich *et al.*, J. Chem. Theory Comput. **10**(8), 2927–2941 (2014)] for electronic transport calculations in molecular junctions is presented. The single driving rate, appearing as a fitting parameter in the original methodology, is replaced by a set of state-dependent broadening factors applied to the different single-particle lead levels. These broadening factors are extracted explicitly from the self-energy of the corresponding electronic reservoir and are fully transferable to any junction incorporating the same lead model. The performance of the method is demonstrated via tight-binding and extended Hückel calculations of simple junction models. Our analytic considerations and numerical results indicate that the developed methodology constitutes a rigorous framework for the design of “black-box” algorithms to simulate electron dynamics in open quantum systems out of equilibrium.

## INTRODUCTION

Over the past decade, the study of electron dynamics in open quantum systems out of equilibrium has gained growing attention within the molecular electronics community.^{1–9} Many important aspects of molecular junctions have been addressed including the characterization of transient current dynamics,^{10–12} dynamical response to time-dependent bias voltages,^{13,14} optically induced current variations,^{15–24} and non-equilibrium thermodynamics in externally driven systems.^{25–28} Despite the many developments made in the field, understanding dynamical effects in molecular junctions remains a major theoretical, computational, and experimental challenge. This challenge needs to be addressed at a fundamental level to enable the future design of molecular-based electronic components with fast response times.

In recent work we presented the driven Liouville-von Neumann (DLvN) equation of motion (EOM) for simulating time-dependent electron transport in molecular junctions.^{3,23,29–33} Within this approach the molecular junction is represented by a fully atomistic finite model system consisting of two sufficiently large lead sections bridged by an (extended) molecule. Open boundary conditions are enforced by augmenting the Liouville-von Neumann EOM with an appropriate non-unitary source/sink term. The latter continuously drives the leads’ state occupations toward the equilibrium Fermi-Dirac distribution of the (implicit) electronic reservoir to which each lead is coupled. With appropriate choices for the chemical potential and the electronic temperature of the various reservoirs, a non-equilibrium charge-polarized state, characterized by charge accumulation and depletion near the corresponding junction model edges, is achieved. This, in turn, results in well-defined voltage and electronic temperature gradients that induce dynamic current flow through the system. The performance of the DLvN approach was demonstrated for simple molecular junctions based on tight-binding (TB) Hamiltonian models^{29,30} as well as for explicit non-orthogonal basis-set representations based on extended Hückel (EH) theory.^{31} The dynamics obtained by the DLvN EOM were shown to conserve density matrix positivity and to obey Pauli’s exclusion principle.^{29–32} Furthermore, the method was shown to accurately describe dynamic currents in junctions subjected to time-dependent perturbations.^{23}

The success of the DLvN approach suggests that it offers an efficient and physically motivated methodology for time-dependent charge transport in molecular junctions. Nevertheless, the original theory incorporates a single fitting parameter that hinders its complete first-principles implementation.^{34} This parameter is a driving factor that dictates the rate at which the system is driven out of equilibrium. While the results were shown to be fairly insensitive to the exact value of the driving rate,^{29} in practice it has to be fitted to reproduce reference steady-state Landauer currents. Recently, a procedure to further reduce the sensitivity of the calculated currents to the driving rate was suggested within the framework of time-dependent density functional theory.^{34} Interestingly, for simple model systems it was shown that an appropriate value of the driving rate can be deduced from physical considerations involving wave packet reflection time scales at the boundaries of the finite junction model.^{29} This suggests that the driving rate should be generally attainable from first-principles considerations.

It is the purpose of the current study to present a rigorous methodology that replaces the single driving rate appearing in the original theory by a set of single-particle lead state broadening factors that are extracted explicitly from the self-energy of the corresponding reservoir. This eliminates the need for any fitting procedure, resulting in an autonomous methodology that can be readily implemented within existing packages in a “black-box” fashion.

## METHODOLOGY

Within the DLvN approach, a finite junction model is formally divided into lead sections that are bridged by an extended-molecule region (see Figure 1(a)). The latter incorporates the molecule augmented by lead sections that are sufficiently large such that the electronic properties at their far edges are converged to those of the infinite lead. Assuming a two-lead setup, the original EOM has the following block form:^{29}

where *H* and *𝝆* are the Hamiltonian matrix representation and density matrix of the entire finite junction, respectively, given in the eigenstate basis of the individual system sections (see Figure 1(b)); $\bm{\rho}L/R$ are the density matrix blocks of the left and right lead sections, respectively; $\bm{\rho}L/R,EM$ and $\bm{\rho}EM,L/R$ represent the lead/extended-molecule coherences; $\bm{\rho}L/R0$ are diagonal target density matrices, whose fixed diagonal elements follow the Fermi-Dirac occupation distribution sampled at the left/right lead eigenvalues with chemical potential and electronic temperature of the corresponding reservoir to which the lead section is implicitly coupled; $\u210f$ is the reduced Plank constant; and $\Gamma /\u210f$ is a real positive number determining the rate at which the leads are driven toward their target equilibrium states.

In the energy domain, the choice of a single driving rate corresponds to the application of a uniform Lorentzian broadening to the entire single-particle lead state manifold. This broadening represents the effect of coupling the finite lead section to an implicit semi-infinite electronic reservoir within the wide-band approximation. In practice, in order to obtain a continuous lead density of states (DOS) in the transport energy window,^{35} the discrete DOS should be slowly varying in this range such that the broadening factor can be chosen as the corresponding typical level spacing. Furthermore, only when the different leads present similar electronic structure can the same broadening factor be applied to all.

In order to eliminate the ambiguity involved in the choice of the driving rate and to represent the effect of coupling the leads to external reservoirs more accurately, we introduce state-dependent broadenings to each lead level. To this end, we employ Green’s function (GF) theory, where the effects of coupling a subsystem (in our case the finite lead section) to its environment (the semi-infinite reservoir) are described by augmenting the Hamiltonian of the subsystem with a self-energy operator, in which information as to the electronic structure of the environment and its coupling to the subsystem is encoded. The real part of the self-energy operator constitutes level shifting within the subsystem, whereas the imaginary part induces level broadening due to the finite life-time of the coupled states.

To obtain the self-energy operator, explicit reservoir models are constructed as duplicates of the finite lead sections (see Figure 1(c)).^{36} The energy dependent matrix representation of the retarded surface Green’s function of each isolated reservoir, $GResr,0(\epsilon )=[\epsilon SRes\u2212HRes0]\u22121$, is first calculated. Here, $HRes0$ is the Hamiltonian matrix representation of the uncoupled semi-infinite reservoir in a given basis-set, ** SRes** is the corresponding overlap matrix, and $\epsilon =E+i\eta $, where

*E*is a real variable representing the electron energy and $\eta \u21920+$ is a small imaginary part added to soften singularities (see the supplementary material for a discussion of the sensitivity of the calculated state-dependent broadenings to the value of

*η*). To this end, we use the iterative principal layer approach,

^{37–43}where the semi-infinite reservoir is sliced into finite sections (termed layers, see Figure 1(c)) that are sufficiently wide such that interactions beyond nearest-neighboring layers can be neglected. Under the above conditions, $\mathbf{H}Res0$ assumes a block tri-diagonal form, with the dimension of each block being that of the principal layer. An efficient iterative algorithm is then used to successively augment the surface layer by bulk layers until convergence of the surface electronic properties, to within the required accuracy, is reached. The resulting $\mathbf{G}Resr,0(\epsilon )$ has the dimension of one principal layer. Next, the reservoir’s retarded self-energy matrix is constructed via $\mathbf{\Sigma}Resr(\epsilon )=(E\mathbf{S}l,Res\u2212\mathbf{V}l,Res)\mathbf{G}Resr,0(\epsilon )(E\mathbf{S}Res,l\u2212\mathbf{V}Res,l)$, where $\mathbf{S}l,Res$ and $\mathbf{V}l,Res$ are the real-space lead/reservoir overlap and coupling matrix blocks, respectively, calculated between the finite lead section and its adjacent reservoir principal layer (see Figure 1(c)), with $\mathbf{S}Res,l=\mathbf{S}l,Res\u2020$ and $\mathbf{V}Res,l=\mathbf{V}l,Res\u2020$.

For each bare lead state, $|l\u27e9$, the reservoir’s self-energy, evaluated at the corresponding energy $\epsilon l$, is added to the finite lead model Hamiltonian $\mathbf{H}Lead+\mathbf{\Sigma}Resr(\epsilon l)$ in the real-space (site) representation.^{44,45} The resulting dressed Hamiltonian is then diagonalized and the complex eigenvalue corresponding to the original bare lead state $|l\u27e9$ is extracted, by comparing the eigenvectors of the dressed Hamiltonian to those of the bare lead Hamiltonian (more details regarding this procedure are provided below). The imaginary part (multiplied by −2)^{45} of this eigenvalue is chosen as the (energy independent) broadening factor of state $|l\u27e9$ (see Appendix A). We note that by evaluating the reservoir’s self-energy at the energy of the given lead eigenstate a “local wide-band approximation” is invoked, where $\mathbf{\Sigma}Resr(\epsilon )$ is assumed to be weakly dependent on energy in the vicinity of $\epsilon l$ (see the supplementary material for a discussion regarding the validity of this approximation).^{46,47} It should be emphasized that unlike the full wide-band approximation, which assumes constant (i.e., independent of the electron energy) and uniform broadening of all lead states, our approach results in constant but state-dependent broadening factors. The latter, which are explicitly derived from the reservoir’s self-energy, allow for a more reliable representation of the electronic structure of the semi-infinite lead by the finite lead model.

To identify which dressed lead state corresponds to the bare lead state $|l\u27e9$, we turn on the lead/reservoir coupling gradually.^{44} To this end, we multiply the lead/reservoir overlap and coupling matrices ($\mathbf{S}l,Res$ and $\mathbf{V}l,Res$) by a scaling factor that is progressively ramped up from 0 to 1. At each step, *k*, the (left or right) eigenfunction, $|lk\u27e9$, of the new dressed Hamiltonian that has the largest overlap, $|\u27e8lk|lk\u22121\u27e9|$, with the eigenfunction that was associated with the bare lead state $|l\u27e9$ in the previous step, $|lk\u22121\u27e9$, is identified and stored.^{45,48} This adiabatic connection assumption allows us to follow the evolution of the bare lead state continuously into the fully dressed one. As a complementary test we also follow the gradual shift in the real-part of the lead eigenvalue with increasing coupling to the reservoir, allowing us the association of the bare lead state with its fully dressed counterpart. These shifts are usually found to be small relative to the full bandwidth of lead states and are hence neglected throughout the rest of the calculation (see the supplementary material for a discussion regarding the validity of this approximation).

The procedure described above is repeated for all eigenstates of a given finite lead section, resulting in a set of state-specific broadening factors. Since the calculation requires just the bare lead overlap and Hamiltonian matrices as input, it is performed only once per given lead. The obtained broadening factors are then fully transferable to every molecular junction using the same lead model.

Given the set of state broadening factors, the DLvN EOM is written as follows (see Appendix B):

where $\mathbf{\Gamma}L$ and $\mathbf{\Gamma}R$ are diagonal matrix blocks of dimensions of the finite left and right lead model basis-sets, respectively, containing the calculated state-specific broadening factors on their diagonals. Notably, this EOM naturally reduces to the original DLvN EOM of Eq. (1) when a single uniform driving rate is used (see Appendix B).

## TESTING AND VALIDATION

To validate the developed parameter-free DLvN (PF-DLvN) methodology, we performed a set of time-dependent transport calculations on several test systems. We start from the simplest one-dimensional TB junction model, consisting of two TB chains representing the leads bridged by a third TB chain acting as the active molecular region (all model parameters appear in the caption of Figure 2). The calculated state-dependent broadening factors for lead models consisting of 200 and 300 sites are presented in Figure 2(a). The broadening factors form a band that reaches a maximum at the lead’s Fermi energy and vanishes near the band edges. This can be rationalized by invoking Fermi’s golden rule for a single lead level coupled to the entire reservoir state manifold.^{49} Here, the lead level broadening is given by $\Gamma i=2\pi \u210f|\u27e8f|V^|i\u27e9|2\delta (Ef\u2212Ei)$, where $\u27e8f|V^|i\u27e9$ is the coupling matrix element between lead state $|i\u27e9$ of energy $Ei$ and reservoir state $|f\u27e9$ of energy $Ef$, and $\delta (x\u2212x0)$ is the Dirac delta function. Therefore, within this approximation the broadenings follow the lead/reservoirs coupling scheme that forms a Newns-Anderson type band.^{30,50} Increasing the lead size by a factor of 1.5 without modifying its bandwidth results in a corresponding increase of the density of states. This is followed by a reduction of the calculated broadening factors (see Figure 2(a)) such that the 200 and 300 sites lead models present an overall similar continuous density of states. In both cases, it is clearly seen that the choice of a single uniform driving rate is justifiable in a small energy window around the Fermi energy but should be done with care when larger bias voltages, approaching the lead band edges, are applied.

The simulated time-dependent currents for a bias voltage of 0.3 V, applied by a symmetric shift of the chemical potentials of the two reservoirs in the target density matrices, are presented in Figures 2(b) and 2(c) for strong and weak lead/molecule coupling, respectively. Polarized initial conditions are employed, where the leads are populated according to their target level occupations and the molecule levels are filled up to the Fermi energy of the entire finite model system. Using the calculated broadening factors both the 200 and 300 lead site model systems present very similar current traces that reproduce the Landauer steady-state current and the plateau regions in the microcanonical^{1,2,4–7} simulation. The latter are performed by running unitary LvN dynamics starting from the same non-equilibrium initial polarized conditions. This provides a clear indication of the ability of the developed methodology to reliably simulate current dynamics in open quantum systems with no parameter fitting involved.

Next, we turn to examine the performance of the PF-DLvN method for non-orthogonal basis-set representations of molecular junctions, within the framework of EH theory.^{31} To this end, a Slater orbital basis-set is introduced to evaluate the overlap and Hamiltonian matrix representation elements using the procedure and parameters described in Ref. 31. First, we consider a hydrogen chain junction model equivalent to the TB chain model discussed above. Here, two hydrogen chain leads, 300 atoms in length each, are bridged by a third 110 atom hydrogen chain acting as the extended-molecule region. The time-dependent current obtained at a bias voltage of 0.5 V and a uniform inter-atomic distance of 2 Å is presented in Figure 3(a). Here, as well, the PF-DLvN current trace matches well the microcanonical simulation results up to the point where the latter exhibits current inversion due to wave-packet reflection from the edges of the finite lead models. Furthermore, the steady-state current corresponds well with the Landauer current (X mark) calculated using the procedure described in Ref. 31. The broadening factors, plotted in the inset of Figure 3(a), form a band similar to the one obtained in the TB calculation (Figure 2(a)), with slight asymmetry around the Fermi energy (−13.5 eV) due to the non-uniform density of states obtained within the EH model.

A somewhat more involved picture is obtained when considering finite-width lead models. Figure 3(b) presents the current dynamics obtained for a 110 atoms hydrogen chain bridging two finite-width hydrogen lead models constructed from five 100 atoms long rows each. Nullification of all broadening factors results in microcanonical dynamics that, in the present case, suffers from strong current oscillations. This indicates that the finite lead models used are insufficient to obtain a stable quasi-steady-state. To remedy this, much larger finite leads have to be modeled, thus hindering the practical applicability of the microcanonical approach to a reliable description of current dynamics in such molecular junctions. Using the PF-DLvN approach, after a typical initial rise the current trace gradually approaches the Landauer results, indicating the suitability of the calculated broadening factors for driving the system toward the required non-equilibrium state. Unlike the atomic chain leads discussed above, however, here these broadening factors form five bands (see inset of Figure 3(b)) around the Fermi energy (−13.09 eV), which correspond to the five single-particle states spanned along the direction perpendicular to the main lead axis by the minimal basis-set used.^{30}

As a last test case we consider a carbon chain junction model, where each atomic center contributes four (one 2s and three 2p) Slater orbital basis functions. The leads are modeled by 200 carbon atom chains each bridged by a 110 atom extended-molecule carbon chain. The lead broadening factors of this system form three bands (see Figure 4(a)). The two higher broadening bands, appearing in the ranges (−23.5)–(−18.6) eV and (−13.1)–(−8.5) eV, correspond to the valence and conduction bands of the system, respectively. The lower band appearing between −11.8 and −10.8 eV corresponds to the increased density of states in the vicinity of the Fermi energy at −11.6 eV due to the appearance of van Hove singularities.

The resulting PF-DLvN time-dependent current approaches the Landauer reference value at steady-state (Figure 4(b)), providing additional support for the physical validity of the developed methodology. This is furtheremphasized in light of the microcanonical simulation results that present a complex set of current plateaus with abrupt current variations (Figure 4(b)). This reflects the different time scales associated with the dynamics of various lead states as manifested by the three broadening factor bands in Figure 4(a). While one of the obtained plateaus indeed matches the correct steady-state value it is difficult to assess, without *a priori* knowledge, when the system attains the correct quasi-steady-state.

Finally, it should be noted that all systems considered herein consist of identical lead models that present relatively uniform DOS within the Fermi transport window. This suggests that the use of a single driving rate should be sufficient to sustain appropriate open boundary conditions. Indeed, if one uses the maximal broadening value calculated for each system using the $\mathbf{H}+\mathbf{\Sigma}$ procedure described above as the single driving rate, the resulting current traces and steady-state occupations are comparable to those obtained using the full PF-DLvN procedure (see the supplementary material). Nevertheless, more realistic model systems, like those considered in Ref. 31, often exhibit complex non-uniform DOS in the vicinity of the lead Fermi energy. For such systems, the use of a single driving rate per lead is insufficient for obtaining an appropriate description of the lead surface electronic properties with a finite lead model. Hence, the entire set of state-dependent broadening factors should be used to drive the system out of equilibrium. Notably, the computational overhead associated with using state-dependent broadenings over the single driving rate during the dynamics is negligible.

## SUMMARY AND CONCLUSIONS

A parameter-free version of the driven Liouville-von Neumann approach to time-dependent simulation of non-equilibrium electron dynamics in open quantum systems was presented. The computational overhead, with respect to the single fitting parameter version of the method, involves the calculation of the reservoir self-energy operator and a diagonalization of the corresponding dressed lead Hamiltonian block for the evaluation of each state-dependent broadening factor. The whole procedure is performed only once per lead/reservoir model, and the obtained set of broadening factors is fully transferable to any molecular junction using the same lead model. The performance of the method was demonstrated on tight-binding and non-orthogonal basis-set representations of simple molecular junction models, including hydrogen and carbon atomic chains bridging one-dimensional and finite-width lead models. The developed methodology eliminates the need for any fitting procedure. This potentially allows for a black-box implementation within the framework of advanced electronic structure methods, at least under the assumption that single-electron states can be interpreted as approximate quasi-particle states.^{52–54}

## SUPPLEMENTARY MATERIAL

See the supplementary material for a discussion regarding the validity of the assumption that the state-specific broadening factors can be approximated to be energy independent, a comparison of simulation results obtained using state-specific and uniform broadening factors, and an assessment of the sensitivity of the results to the value of the parameter *η* used in the reservoir’s surface Green’s function calculation.

## ACKNOWLEDGMENTS

Work at TAU was supported by the Israel Science Foundation under Grant No. 1740/13, the Raymond and Beverly Sackler Fund for Convergence Research in Biomedical, Physical, and Engineering Sciences, the Lise Meitner-Minerva Center for Computational Quantum Chemistry, and the Center for Nanoscience and Nanotechnology at Tel-Aviv University. Work at Weizmann was supported by the Israel Science Foundation and the Lise Meitner-Minerva Center for Computational Quantum Chemistry. Work at the University of Copenhagen was supported by the Lundbeck Foundation. Work in Berkeley was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering under Contract No. DE-AC02-05CH11231. Work performed at the Molecular Foundry was also supported by the Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under the same contract number.

### APPENDIX A:DERIVATION OF THE $\mathbf{H}+\Sigma $ APPROACH FOR CALCULATING LEAD STATE-DEPENDENT BROADENING FACTORS

The accurate description of quantum dynamics in large systems is usually a prohibitively demanding computational task. Nevertheless, often it is sufficient to focus on the dynamics of a small section of the entire system. In these cases, the equation of motion can be recast as a sum of terms describing the internal subsystem dynamics and the influence of its coupling to the rest of the system. In the context of the present study, we are interested in describing electron dynamics within a finite lead model that is coupled to an external semi-infinite reservoir. The contribution of the reservoir to the lead dynamics, $L(t)$, can be expressed in terms of the lesser Green’s function, $\bm{G}<(t,t)$, and Hamiltonian, ** H**, matrix representations of the entire lead/reservoir system as follows:

^{23,55}

where *l* and *m* are single-particle state indices within the finite lead model and the index *r* runs over all reservoir states. The lesser GFs can be factorized using Langreth rules to obtain^{23,55}

where $\bm{G}adv(\tau ,t)$, $\bm{G}ret(t,\tau )$, $\bm{G}(0),adv(\tau ,t)$, $\bm{G}(0),ret(t,\tau )$ are the advanced and retarded GFs of the entire system in the presence and absence of lead/reservoir coupling, respectively, and $\bm{G}(0),<(\tau ,t)$ is the lesser GF of the entire uncoupled lead/reservoir system. The indices *n* and *s* run over all lead and reservoir states, respectively.

We can now define the retarded, advanced, and lesser self-energies as follows:

and rewrite Eq. (A3) as

#### 1. Sink term

The first sum in Eq. (A5) can be identified as the sink term absorbing electrons traveling in the lead toward the reservoir,

The retarded and advanced self energies appearing in Eq. (A6) can be recast as Fourier integrals of the form

where we have used the definitions appearing in Eq. (A4) and the fact that we are working in the basis of eigenfunctions of the isolated lead and reservoir subsystems (state representation), where the isolated reservoir’s retarded and advanced GFs are diagonal.

Invoking the wide band approximation (WBA), we can neglect the real-part of the isolated reservoir’s retarded and advanced GFs, associated with level shifting. Notably, in all numerical examples appearing in the main text, lead level shifts due to coupling to the reservoir were found to be very small with respect to the full lead energy bandwidth and level crossing was rarely observed (see the supplementary material for further discussion regarding the effect of coupling to the reservoir on the lead level shifts). This indicates the validity of this approximation. The remaining imaginary part can be written as

Eq. (A9) can be rewritten in terms of the broadening matrix element,

as

The WBA further implies that $\Gamma lnlead(\epsilon )$ is energy independent (see the supplementary material for a discussion regarding the validity of this assumption) for all indices *l* and *n* so it can be taken out of the integration in Eq. (A11) to give

where we have introduced the reduced density matrix of the finite lead model $\rho lead(t)=\u2212iG<(t,t)$.

Eq. (A13) can be written in matrix form as

If we further assume that the lead states are not mixed by the coupling to the reservoir, $\mathbf{\Gamma}lead$ becomes a diagonal matrix holding the state-dependent (and energy independent) broadening factors and Eq. (A14) becomes identical to the lead block damping term appearing in Eq. (2) above.

#### 2. Source term

The second sum in Eq. (A5) can be identified as the source term injecting electrons from the reservoir traveling toward the extended-molecule through the lead,

The lesser self-energies appearing in Eq. (A15) contain information regarding the reservoir states populations. To obtain this information explicitly we first assume, as above, that the lead states are not mixed by the bath such that both the lesser Green’s function appearing in Eq. (A4) and the lesser self-energy itself are diagonal. Next, the diagonal elements of the lesser self-energy can be expressed in terms of a Fourier integral,

Since the lesser GF appearing in Equation (A16) relates to the uncoupled reservoir at thermal equilibrium, it can be expressed in terms of the advanced and retarded GFs as

where the spectral function of the uncoupled reservoir, $Arr(\epsilon )$, consists of a set of Dirac $\delta $ functions centered at the corresponding eigenvalues. Inserting (A17) in (A16) gives

Using this in the expression for the source term in Eq. (A15), one obtains

Next, we assume that the lead Hamiltonian is time-independent such that the retarded and advanced GFs appearing in Eq. (A19) have time-translational invariance and can be written as the following Fourier integral:

Inserting these expressions in the time integrals appearing in Eq. (A19) we obtain

and

where we have defined $\omega r\u2261\epsilon r\u210f$. Inserting Eqs. (A21) and (A22) in Eq. (A19) yields

Following the assumption that the lead state are not mixed by the coupling to the reservoir we focus on the diagonal elements of the source term appearing in Eq. (A23),

where we have introduced the spectral function $All(\epsilon r)=i[Gllret(\epsilon r)\u2212Glladv(\epsilon r)]$ and used the properties of the Dirac $\delta $ function in the last equality. Using the definition of the broadening matrix elements of Eq. (A10) above, Eq. (A24) can be written as

Since the spectral function, $All(\epsilon )$, is peaked around the lead eigenvalue $\epsilon l$, we can evaluate the broadening at this eigenvalue $\Gamma lllead(\epsilon )\u2248\Gamma lllead(\epsilon l)$ and take it out of the integral to obtain the equilibrium density matrix (see the supplementary material for a discussion regarding the validity of this approximation),

Defining the diagonal matrix $\mathbf{\Gamma}lead$ whose diagonal elements host the values of $\Gamma lllead(\epsilon l)$ and the diagonal matrix $\bm{\rho}0$ whose diagonal elements have the Fermi-Dirac distribution sampled at the isolated lead’s eigen-energies $f(\epsilon l)$ encoding the chemical potential and electronic temperature of the corresponding reservoir, we can rewrite Eq. (A26) as

This is identical to the lead block injection term appearing in Eq. (2).

### APPENDIX B: HEURISTIC DERIVATION OF THE DRIVEN LIOUVILLE-VON NEUMANN EQUATION WITH STATE-SPECIFIC BROADENING FACTORS

To obtain the DLvN equation in the case of state-dependent broadening factors, we follow the heuristic derivation of the original equation with a single driving rate, appearing in Ref. 29. First, the Hamiltonian of the entire finite junction model is written in the basis of eigenfunctions of the leads and the extended-molecule section. The procedure to transform the real-space (site) Hamiltonian to the energy (state) representation in the case of orthogonal and non-orthogonal basis-set representations is given in Refs. 29 and 31, respectively. Next, to account for absorption of outgoing electrons propagating through the leads towards the edges of the finite lead models, the Hamiltonian blocks representing the lead sections are augmented by diagonal imaginary matrices $\u2212i\bm{\gamma}L/R$ representing the broadening (or damping rate) of each single-particle lead state. For a two-lead setup this translates to the following form:

Inserting this expression in the Liouville-von Neumann Eq. (A8) of Ref. 29 and assuming that in the state representation $\bm{\gamma}L$ and $\bm{\gamma}R$ are diagonal blocks, the absorption term obtains the following anti-commutation form (see Eq. (A9) of Ref. 29):

The emission term is obtained in a similar manner by replacing the density matrix of the system with the density matrix of the electronic reservoirs and inverting the overall sign (see Eq. (A12) of Ref. 29),

where $\bm{\rho}L/R0$ are diagonal matrices holding equilibrium Fermi-Dirac single-particle state occupations on their diagonal with the chemical potential and electronic temperature of the left/right reservoirs.

Summing the absorption (Eq. (B2)) and emission (Eq. (B3)) contributions results in the driving term appearing in Eq. (2) of the main text,

Defining $\mathbf{\Gamma}L=2\bm{\gamma}L$ and $\mathbf{\Gamma}R=2\bm{\gamma}R$, we obtain

When using a uniform driving rate the broadening matrices obtain the structure $\mathbf{\Gamma}L=\Gamma \mathbf{I}L$ and $\mathbf{\Gamma}R=\Gamma \mathbf{I}R$, where $\mathbf{I}L/R$ are unit matrices of dimensions of the left and right lead model basis. In this case, the driving term of Eq. (B5) assumes the form

This is the driving term appearing in the original DLvN EOM of Eq. (1).

## REFERENCES

The energy range between the chemical potential of the two leads.

The choice of full lead model duplicates to construct the reservoir model is done for convenience. In many cases, a considerably smaller principal layer can be used without loss of accuracy. Note that the wider the principal layer used the less iterations are required to converge the electronic properties of the semi-infinite reservoir model but the whole process involves inversion and multiplication of larger matrices.

For the systems considered herein typically 10-20 steps were sufficient to follow the dressed lead levels progress. As the lead DOS becomes more complex and level-crossing occurs more steps will be required.