Open-system approaches are gaining traction in the simulation of charge transport in nanoscale and molecular electronic devices. In particular, “extended reservoir” simulations, where explicit reservoir degrees of freedom are present, allow for the computation of both real-time and steady-state properties but require relaxation of the extended reservoirs. The strength of this relaxation, γ, influences the conductance, giving rise to a “turnover” behavior analogous to Kramers turnover in chemical reaction rates. We derive explicit, general expressions for the weak and strong relaxation limits. For weak relaxation, the conductance increases linearly with γ and every electronic state of the total explicit system contributes to the electronic current according to its “reduced” weight in the two extended reservoir regions. Essentially, this represents two conductors in series—one at each interface with the implicit reservoirs that provide the relaxation. For strong relaxation, a “dual” expression-one with the same functional form-results, except now proportional to 1/γ and dependent on the system of interest’s electronic states, reflecting that the strong relaxation is localizing electrons in the extended reservoirs. Higher order behavior (e.g., γ2 or 1/γ2) can occur when there is a gap in the frequency spectrum. Moreover, inhomogeneity in the frequency spacing can give rise to a pseudo-plateau regime. These findings yield a physically motivated approach to diagnosing numerical simulations and understanding the influence of relaxation, and we examine their occurrence in both simple models and a realistic, fluctuating graphene nanoribbon.

“Extended reservoir” simulation recognizes that there is a hierarchy of length and time scales in transport. Particles (electrons, phonons, etc.) flow from very large reservoirs—essentially, external sources or sinks—into smaller, more confined regions before flowing through some “system of interest,” a molecular junction, a nanotube, a graphene nanoribbon, etc. This concept gives rise to provably accurate simulation approaches that incorporate part of the reservoir explicitly into the simulation.1–3 For electron transport, this will give exactly the general results of Jauho, Meir, and Wingreen4–6 for interacting and non-interacting systems alike (or, in the case of non-interacting systems, the Landauer formula7,8).

This approach requires the simulation of a larger system overall—not only does one simulate the “system of interest,” S, and possibly some small metallic leads, but rather that together with many extra degrees of freedom in the left and right “extended reservoirs,” denoted L and R—yet it yields a vast simplification for time-dependent phenomena. The two-time Green’s functions can be replaced with a time-local (Markovian) master equation,

ρ̇=ı[HLSR,ρ]+kLRγk+ckρck12ckck,ρ+kLRγkckρck12ckck,ρ,
(1)

so long as the relaxation is weak enough (but-in order to extract properties of S-not too weak, i.e., not in the small-γ regime we examine here2,3). The density matrix, ρ, is for the LSR system, and ck(ck) are the creation (annihilation) operators for state k. The Hamiltonian is

HLSR=HS+HL+HR+HI.
(2)

where HS is for S (including, potentially, many-body interactions), HL(R)=kL(R)ωkckck are for the “extended reservoirs,” and HI=kLRiS(vkickci+ h.c.) is their interaction. The index k includes all degrees of freedom (electronic state, spin, reservoir), while ωk and vki denote the level and hopping frequencies, respectively.

The first term in Eq. (1) describes Hamiltonian (unitary) dynamics of the LSR system. The terms outside of the commutator reflect particle injection (depletion) into the state k at a rate γk+ (γk). These will relax the extended reservoirs to equilibrium—a pseudo-equilibrium in terms of the isolated extended reservoir states2,3—if HI is removed when γk+γkfα(ωk) and γkγk[1fα(ωk)], where fα(ωk) is the Fermi-Dirac distribution in the α{L,R} reservoir.

Generically, one wants to simulate structures such as the one shown in Fig. 1(a), some molecule or nanoscale device (here, a graphene sensor) between two electronic reservoirs [see Fig. 1(b) for a schematic]. These types of structures are of interest to, e.g., sensing9–14 and sequencing.15–24 Within the extended reservoir approach, one models this setup by dividing the total system into three parts, the junction S (molecule/structure, possibly including some part of the metallic leads) and the electronic reservoirs. The latter will be further split into the explicit degrees of freedom L/R and implicit reservoirs that ensure proper sources/sinks are present. Considering one extended reservoir mode k, its relaxation—its connection to the implicit reservoir Ek—is given by the Hamiltonian

Hk=ωkckck+αEkωαcαcα+αEktαcαck+ckcα,
(3)

with tα the coupling to the reservoir mode (which is related to the relaxation rate γ2,3). The implicit reservoir is in equilibrium at some temperature and chemical potential (which are different in the left and right regions), thus providing an infinite environment that will relax the mode k to equilibrium in the absence of S. Including the environments Ek gives a non-Markovian generalization of Eq. (1). The Markovian approximation will be valid in the wide-band limit for Ek and when γkBT/, where kB is Boltzmann’s constant and T the temperature.2,3

FIG. 1.

(a) A fluctuating graphene nanoribbon in aqueous solution between two gold substrates (water is not shown for clarity). (b) Schematic of the electronic modes of a LSR system, where the couplings and onsite energies are arbitrary. Relaxation occurs only in the L and R regions, signified by the γk+ and γk, where the superscript L(R) is to explicitly indicate that the mode k is in the left or right region (which are at different chemical potentials). (c) For the small-γ limit, the whole LSR system is diagonalized to get a set of global modes. These determine the current, as the interface between LSR and the implicit reservoirs is the rate limiting process. (d) For the large-γ limit, each of the regions L, S, and R are separately diagonalized, as the interface between L(R) and S determines the current.

FIG. 1.

(a) A fluctuating graphene nanoribbon in aqueous solution between two gold substrates (water is not shown for clarity). (b) Schematic of the electronic modes of a LSR system, where the couplings and onsite energies are arbitrary. Relaxation occurs only in the L and R regions, signified by the γk+ and γk, where the superscript L(R) is to explicitly indicate that the mode k is in the left or right region (which are at different chemical potentials). (c) For the small-γ limit, the whole LSR system is diagonalized to get a set of global modes. These determine the current, as the interface between LSR and the implicit reservoirs is the rate limiting process. (d) For the large-γ limit, each of the regions L, S, and R are separately diagonalized, as the interface between L(R) and S determines the current.

Close modal

Both Eq. (1) and its full non-Markovian generalization have exact, closed form solutions for their steady-state currents, akin to the Meir-Wingreen formula.2 Our purpose here is to derive the small- and large-γ limits of the steady-state for both the non-Markovian and Markovian cases. We deal solely with non-interacting electrons since this yields compact, illustrative expressions.

For the non-Markovian case, the total current is found by integrating out the environment in Eq. (3),2,5 giving

I=e2πdωfL(ω)fR(ω)  ×trΓL(ω)Gr(ω)ΓR(ω)Ga(ω),
(4)

where fL(R) is the equilibrium distribution in the left (right) reservoir, ΓL(R) are the spectral densities, and Gr(a) are the retarded (advanced) Green’s functions (with bold capital letters indicating matrices). The environment being integrated out can either be the implicit reservoirs only or the implicit reservoirs plus L and R. Both will give the same result, but approximations will be easier using different forms in the different regimes, a fact that we will use below. The current, Eq. (4), in the presence of relaxation has a turnover behavior1–3 I1  smallγI2I3  largeγ, where it first increases (regime 1), then plateaus (2), and then decreases (3).

Weak relaxation—When γ is small, we will use the Green’s functions for LSR in Eq. (4). For HL and HR separately diagonalized, the self-energies are

ΣijL=ı2kLγkδijδik,  ΣijR=ı2kRγkδijδik,
(5)

where γ is nonuniform.2,25,26 These are matrices on the whole LSR system but are zero when i or j are outside the respective reservoir region. The spectral and Green’s functions are ΓL(R)=2ImΣL(R) and

Gr=1ωHLSRΣLΣR,
(6)

respectively. We rotate these expressions into the eigenbasis for LSR through the unitary transform U, i.e.,

UHLSRUij=ω̃iδij,
(7)

where ω̃i are the eigenvalues of HLSR. The transformation on the spectral functions and self-energies gives

Γ̃L(R)=UΓL(R)U,  Σ̃L(R)=UΣL(R)U.
(8)

The quantity ωHLSR is now diagonal and the remaining terms in the denominator of Gr are controlled by γ. The dominant terms are thus the diagonal components in this basis, as these diverge as 1/γ, yielding

(Gr)ij1ωω̃iΣ̃iiLΣ̃iiRδij,
(9)

with i,jLSR. The off-diagonal components of the self-energies yield a higher order correction to this. This creates a sharply peaked Lorentzian to the lowest order, which can be integrated directly.27 Using Eq. (9) in Eq. (4) yields a double sum over eigenstates of LSR. The cross terms, dω(ωω̃iıγ/2)1(ωω̃j+ıγ/2)1=(πγ2/2)/(γıω̃i+ıω̃j), in this sum behave as γ2 when ω̃iω̃j and γ is small compared to the spacing ω̃iω̃j. In this limit, the total current, I1, is

I1e4iLSRuiL2uiR2uiL2+uiR2fL(ω̃i)fR(ω̃i)γ,
(10)

where

uiL(R)2=kL(R)γk|Uik|2
(11)

is the γ-weighted extent of the electronic modes in the extended reservoir regions. Equation (10) is a more general result than that derived in classical thermal transport to understand the influence of topological edge modes on the conductance.28,29 This small-γ limit is the same for the non-Markovian and Markovian [Eq. (1)] cases.

For low temperatures (on the electronic scale) and uniform reduced couplings (for i in the bias window), Eq. (10) is proportional to the number of modes in the bias window. Thus, except for finite size effects [e.g., that change the mode structure and therefore change the form of the matrix elements in Eq. (11)], the current for small γ grows linearly with the number of modes in the extended reservoir. This means that the plateau region—the flat region that follows the small-γ regime, I2—grows with the size of the extended reservoirs [and, indeed, this is why the Markovian approximation becomes accurate for simulating the Meir-Wingreen/Landauer conductance of S,2,3 i.e., the Markovian approximation is valid for γ larger than that needed for the validity of Eq. (10)]. However, as we will see, inhomogeneity in the frequency spacing of LSR modes can give rise to undesirable oscillatory features in between the small- and intermediate-γ regimes (essentially, a symptom of finite sizes).

Strong relaxation—When γ is large, it is advantageous to work with Gr of S only, giving the spectral functions

ΓijL(R)(ω)=kL(R)vikvkjγk(ωωk)2+γk2/4,
(12)

where the v’s are the couplings in HI. The highest order (in 1/γ) term comes from the approximation

ΓijL(R)(ω)kL(R)4vikvkjγk.
(13)

That is, the relaxation must dominate all other energy scales. In practice, this means γkW, where W is the bandwidth of the extended reservoirs (e.g., all |ωk|<W), where we assume that the bias and temperature are small (so that the integration over ω is much smaller than γ30).

For just S, the Green’s function is

Gr=1ωHSΣLΣR,
(14)

where ΣL(R) are not given by Eq. (5) but rather are the self-energies of L(R) on S. Once again, we diagonalize the Hamiltonian component, i.e., ωHS. The remaining terms are controlled by 1/γ. As for small γ, the off-diagonal terms give a higher order contribution; hence,

(Gr)ij1ωωiΣiiLΣiiRδij,
(15)

where ωi are the frequencies of S in isolation (as compared to the entire LSR system previously). This is again a Lorentzian but now sharply peaked for large γ. The current, I3, in the strong relaxation limit is then

I3e4iSviL2viR2viL2+viR2fL(ωi)fR(ωi)1γ,
(16)

where

viL(R)2=kL(R)4|vik|2γk.
(17)

The current has the same form as the small-γ limit except proportional to 1/γ. As well, for small γ, the extent of the LSR modes across a structure determines the current [Fig. 1(c)], whereas for large γ, the coupling between S and L(R) determines the current [Fig. 1(d)]. The renormalized form of the relevant coupling, Eq. (17), is a reflection of the Zeno effect,31 where electrons are attempting to hop out of the reservoir into the system, but they are being strongly “measured” by the relaxation.

The Markovian case, though, is different than Eq. (16),

Ī3e4iSv¯iL2v¯iR2v¯iL2+v¯iR21γ,
(18)

where

v¯iL(R)2=kL(R)4|vik|2γkfL(R)(ωk).
(19)

This is due to the fact that the Markovian equation of motion populates the reservoir modes according to their isolated frequencies, i.e., it is a pseudo-equilibrium that does not account for the broadening of the density of states. Indeed, the Markovian approximation is not valid in this limit for this very reason.2,3

Equations (10), (16), and (18) are our main results. Figure 2(a) shows the current, I, versus γ for a simple example: A system with only a single mode symmetrically coupled to identically sized linear (1D) reservoirs in L and R all with strength v0. When the system mode has zero onsite energy, the small- and large-γ limits display γ and 1/γ behavior, respectively. However, if the system mode is shifted outside the bias window, the leading term for large γ becomes 1/γ2. Similarly, Fig. 2(b) shows a total LSR system that has no modes in the bias window, giving γ2 dependence for small γ. The error inherent in numerical integration of a Lorentzian can also yield incorrect limits for small γ [Fig. 2(c)].

FIG. 2.

(a) The current, I, computed via the non-Markovian result, Eq. (4), versus relaxation strength, γk = γ, for 1D reservoirs and S being a single electronic mode either within (blue line) or outside (green line) the bias window. Here, the small-γ behavior [blue dots, Eq. (10)] is nearly identical—the shift in the energies of the LSR system is negligible since the number of reservoir modes is large. For large γ, however, the behavior changes from 1/γ [blue dashes, Eq. (16)] to 1/γ2 (green dashes, a fit) when the system’s mode falls outside the bias window. (b) The current in the same 1D model with a bias such that no modes in LSR or S are in the bias window. The dotted and dashed blue lines show fits to γ2 and 1/γ2, respectively. (c) Numerical integration errors for a single Lorentzian. In most cases, an insufficient error tolerance results in the integral behaving as γ2 due to the tendency to exclude the bulk of the peak. Red dashed and cyan dotted lines show a Gaussian quadrature with two different error constraints. This systematic error is seen for most types of weighted integration. Hence, γ2 behavior can be either real or due to numerical issues (easily identified/solved with smaller tolerances or by using methods with a predetermined error, such as Monte Carlo, shown as blue circles).

FIG. 2.

(a) The current, I, computed via the non-Markovian result, Eq. (4), versus relaxation strength, γk = γ, for 1D reservoirs and S being a single electronic mode either within (blue line) or outside (green line) the bias window. Here, the small-γ behavior [blue dots, Eq. (10)] is nearly identical—the shift in the energies of the LSR system is negligible since the number of reservoir modes is large. For large γ, however, the behavior changes from 1/γ [blue dashes, Eq. (16)] to 1/γ2 (green dashes, a fit) when the system’s mode falls outside the bias window. (b) The current in the same 1D model with a bias such that no modes in LSR or S are in the bias window. The dotted and dashed blue lines show fits to γ2 and 1/γ2, respectively. (c) Numerical integration errors for a single Lorentzian. In most cases, an insufficient error tolerance results in the integral behaving as γ2 due to the tendency to exclude the bulk of the peak. Red dashed and cyan dotted lines show a Gaussian quadrature with two different error constraints. This systematic error is seen for most types of weighted integration. Hence, γ2 behavior can be either real or due to numerical issues (easily identified/solved with smaller tolerances or by using methods with a predetermined error, such as Monte Carlo, shown as blue circles).

Close modal

We also apply this method to a suspended graphene nanoribbon [Fig. 1(a)]. In this example, L and R are the gold contacts and S is the ribbon itself. The graphene Hamiltonian is

H=ijv0el/λcicj,
(20)

where v03 eV is the hopping energy between pz orbitals in pristine graphene, l is the carbon-carbon bond stretching, and λ0.047 nm (parameters are a rearrangement of those in Ref. 32). For the substrate, the gold atoms have a −3 eV hopping energy and an onsite energy of −1.6 eV, representing the 6s-conduction band.33 The carbon-gold coupling is of similar form to Eq. (20), with the same v0 and λ0.11 nm.

Figure 3 shows the current as a function of γ for a 50 meV bias and 200 meV Fermi level. The small-γ prefactor generally increases with reservoir size, but the exact form depends on the LSR mode structure. For the smallest reservoir size (blue data points), there are only four LSR modes in the bias window. These “turnover” out of the small-γ regime at about γ3×103 fs−1. Almost simultaneously, a γ2 contribution turns on from a mode just outside the bias window, giving the quasi-linear regime before the plateau. This can be made concrete by defining projection operators onto the subspace of these five modes (or even separate projectors onto the bias window modes and the one just outside the bias window). The blue solid line shows the contribution by only these handful of modes, clearly showing that they give the features in between the small- and intermediate-γ regimes. In addition, the defining role of LSR modes in the transition to the plateau regime allows one to determine when the mode spectrum (i.e., coupling across the whole structure versus frequency) is uniformly converging or still displaying finite-size effects. The green and red data points have a more rapid rise to the plateau (due to the larger number of LSR modes), but inhomogeneity of the mode spacings gives rise to a pseudo-plateau region, where oscillations occur before transitioning into the real plateau.

FIG. 3.

The electronic current, I [from Eq. (4)], through a graphene nanoribbon suspended in aqueous solution between two gold substrates, plotted with an increasing number of Au unit cells (i.e., the contact between the gold and the graphene is kept fixed but the reservoir regions are made deeper). The dashed lines are the exact results for the small-γ limit [Eq. (10)], which generally increase with the additional number of gold layers. However, within this small-γ regime, finite-size effects prevent a simple linear scaling of the current with the extended reservoir size, i.e., the total current is not dependent on solely the size of the extended reservoir but also on the structure of the electronic states. All cases go to the same large-γ limit (dash-dot line) since S is identical. In the intermediate regime, a larger plateau forms for the larger extended reservoirs. However, inhomogeneous mode spacing will create inhomogeneity in the turnover points, which in turn will create features in the small- to intermediate-γ region. That is, even outside the formally small-γ region, the results reflect the distribution in important quantities in the small-γ regime (solid lines show the contribution of the modes in and around the bias window in the LSR basis).

FIG. 3.

The electronic current, I [from Eq. (4)], through a graphene nanoribbon suspended in aqueous solution between two gold substrates, plotted with an increasing number of Au unit cells (i.e., the contact between the gold and the graphene is kept fixed but the reservoir regions are made deeper). The dashed lines are the exact results for the small-γ limit [Eq. (10)], which generally increase with the additional number of gold layers. However, within this small-γ regime, finite-size effects prevent a simple linear scaling of the current with the extended reservoir size, i.e., the total current is not dependent on solely the size of the extended reservoir but also on the structure of the electronic states. All cases go to the same large-γ limit (dash-dot line) since S is identical. In the intermediate regime, a larger plateau forms for the larger extended reservoirs. However, inhomogeneous mode spacing will create inhomogeneity in the turnover points, which in turn will create features in the small- to intermediate-γ region. That is, even outside the formally small-γ region, the results reflect the distribution in important quantities in the small-γ regime (solid lines show the contribution of the modes in and around the bias window in the LSR basis).

Close modal

For both weak and strong relaxation, the current depends on the amplitude of the eigenmodes at the boundaries, but different boundaries and therefore different eigenmodes. Going from small to large γ changes the relevant contact resistance from that between LSR and the implicit reservoirs to that between S and L and R. However, the LSR modes—and hence their spacing and uniform convergence—play the dominant role in the transition from the small-γ to the plateau regime, where the current reflects the intrinsic conductance of S.2,3 We note, of course, that the system of interest does not have a unique conductance, independent of how it is contacted. Thus, the intrinsic conductance reflects the natural conductance of the setup with all γ dependence eliminated.

Open-system simulations of the type we address here are increasingly in use for simulating nanoelectronics.26,34–38 These also have a relation to closed system simulations of real-time dynamics39–41 (including of particle transport in cold-atom setups42–45). To successfully employ this simulation approach, one should be in between the two limits and, in the case of Eq. (1), with γ still weak enough that a Markovian approximation is valid.2,3 The compact analytic forms we derive allow one to benchmark these simulations, as well as make simple predictions for other scenarios (e.g., quantum dot systems or molecules with variable contacts) where the relaxation/tunneling rate can be tuned. Moreover, the physics associated with this whole turnover process is analogous to Kramers turnover in chemical relaxation rates,1,2 thus giving a physically intuitive conceptual paradigm for simulating various transport processes, from thermal/energy to electronic.

See supplementary material for the complete atomic coordinates of the nanoribbon structure used in Fig. 3.

Daniel Gruss acknowledges support under the Cooperative Research Agreement between the University of Maryland and the National Institute of Standards and Technology Center for Nanoscale Science and Technology, Award No. 70NANB14H209, through the University of Maryland. Alex Smolyanitsky gratefully acknowledges support from the Materials Genome Initiative.

1.
K. A.
Velizhanin
,
S.
Sahu
,
C.-C.
Chien
,
Y.
Dubi
, and
M.
Zwolak
,
Sci. Rep.
5
,
17506
(
2015
).
2.
D.
Gruss
,
K. A.
Velizhanin
, and
M.
Zwolak
,
Sci. Rep.
6
,
24514
(
2016
).
3.
J. E.
Elenewski
,
D.
Gruss
, and
M.
Zwolak
, e-print arXiv:1705.00566 (
2017
).
4.
Y.
Meir
and
N. S.
Wingreen
,
Phys. Rev. Lett.
68
,
2512
(
1992
).
5.
A.
Jauho
,
N.
Wingreen
, and
Y.
Meir
,
Phys. Rev. B
50
,
5528
(
1994
).
6.
H.
Haug
and
A.-P.
Jauho
,
Quantum Kinetics in Transport and Optics of Semiconductors
(
Springer
,
1996
).
7.
R.
Landauer
,
IBM J. Res. Dev.
1
,
223
(
1957
).
8.
M.
Di Ventra
,
Electrical Transport in Nanoscale Systems
(
Cambridge University Press
,
Cambridge
,
2008
).
9.
D. M.
Willard
,
T.
Mutschler
,
M.
Yu
,
J.
Jung
, and
A.
Van Orden
,
Anal. Bioanal. Chem.
384
,
564
(
2006
).
10.
Y.
Choi
,
I. S.
Moody
,
P. C.
Sims
,
S. R.
Hunt
,
B. L.
Corso
,
I.
Perez
,
G. A.
Weiss
, and
P. G.
Collins
,
Science
335
,
319
(
2012
).
11.
B. R.
Goldsmith
,
J. G.
Coroneus
,
V. R.
Khalap
,
A. A.
Kane
,
G. A.
Weiss
, and
P. G.
Collins
,
Science
315
,
77
(
2007
).
12.
B. R.
Goldsmith
,
J. G.
Coroneus
,
A. A.
Kane
,
G. A.
Weiss
, and
P. G.
Collins
,
Nano Lett.
8
,
189
(
2008
).
13.
S.
Sorgenfrei
,
C.-y.
Chiu
,
R. L.
Gonzalez
, Jr.
,
Y.-J.
Yu
,
P.
Kim
,
C.
Nuckolls
, and
K. L.
Shepard
,
Nat. Nanotechnol.
6
,
126
(
2011
).
14.
S.
Sorgenfrei
,
C.-y.
Chiu
,
M.
Johnston
,
C.
Nuckolls
, and
K. L.
Shepard
,
Nano Lett.
11
,
3739
(
2011
).
15.
M.
Zwolak
and
M.
Di Ventra
,
Nano Lett.
5
,
421
(
2005
).
16.
M.
Zwolak
and
M.
Di Ventra
,
Rev. Mod. Phys.
80
,
141
(
2008
).
17.
J.
Lagerqvist
,
M.
Zwolak
, and
M.
Di Ventra
,
Nano Lett.
6
,
779
(
2006
).
18.
J.
Lagerqvist
,
M.
Zwolak
, and
M.
Di Ventra
,
Phys. Rev. E
76
,
013901
(
2007
).
19.
J.
Lagerqvist
,
M.
Zwolak
, and
M.
Di Ventra
,
Biophys. J.
93
,
2384
(
2007
).
20.
M.
Krems
,
M.
Zwolak
,
Y. V.
Pershin
, and
M.
Di Ventra
,
Biophys. J.
97
,
1990
(
2009
).
21.
S.
Chang
,
S.
Huang
,
J.
He
,
F.
Liang
,
P.
Zhang
,
S.
Li
,
X.
Chen
,
O.
Sankey
, and
S.
Lindsay
,
Nano Lett.
10
,
1070
(
2010
).
22.
M.
Tsutsui
,
M.
Taniguchi
,
K.
Yokota
, and
T.
Kawai
,
Nat. Nanotechnol.
5
,
286
(
2010
).
23.
E.
Paulechka
,
T. A.
Wassenaar
,
K.
Kroenlein
,
A.
Kazakov
, and
A.
Smolyanitsky
,
Nanoscale
8
,
1861
(
2016
).
24.
A.
Smolyanitsky
,
B. I.
Yakobson
,
T. A.
Wassenaar
,
E.
Paulechka
, and
K.
Kroenlein
,
ACS Nano
10
,
9009
(
2016
).
25.
T. M.
Henderson
,
G.
Fagas
,
E.
Hyde
, and
J. C.
Greer
,
J. Chem. Phys.
125
,
244104
(
2006
).
26.
T.
Zelovich
,
T.
Hansen
,
Z.-F.
Liu
,
J. B.
Neaton
,
L.
Kronik
, and
O.
Hod
,
J. Chem. Phys.
146
,
092331
(
2017
).
27.

While the approximations can be done with the self-energies for both small and large γ, the entire integrand cannot be expanded due to the nature of the integral [Eq. (4)] over the Lorentzian. Some terms must be kept until after the integration is complete.

28.
C.-C.
Chien
,
S.
Kouachi
,
K. A.
Velizhanin
,
Y.
Dubi
, and
M.
Zwolak
,
Phys. Rev. E
95
,
012137
(
2017
).
29.
C.-C.
Chien
,
K. A.
Velizhanin
,
Y.
Dubi
,
B. R.
Ilic
, and
M.
Zwolak
, e-print arXiv:1707.06669 (
2017
).
30.

This is necessary since the Ek broaden the density of states and there is not a hard cutoff to the integration in Eq. (4).

31.
W. M.
Itano
,
D. J.
Heinzen
,
J.
Bollinger
, and
D.
Wineland
,
Phys. Rev. A
41
,
2295
(
1990
).
32.
D. A.
Cosma
,
M.
Mucha-Kruczyński
,
H.
Schomerus
, and
V. I.
Fal’ko
,
Phys. Rev. B
90
,
245409
(
2014
).
33.
M.
Wahiduzzaman
,
A. F.
Oliveira
,
P.
Philipsen
,
L.
Zhechkov
,
E.
van Lenthe
,
H. A.
Witek
, and
T.
Heine
,
J. Chem. Theory Comput.
9
,
4006
(
2013
).
34.
C. G.
Sánchez
,
M.
Stamenova
,
S.
Sanvito
,
D.
Bowler
,
A. P.
Horsfield
, and
T. N.
Todorov
,
J. Chem. Phys.
124
,
214708
(
2006
).
35.
J. E.
Subotnik
,
T.
Hansen
,
M. A.
Ratner
, and
A.
Nitzan
,
J. Chem. Phys.
130
,
144105
(
2009
).
36.
S.
Ajisaka
,
F.
Barra
,
C.
Mejía-Monasterio
, and
T.
Prosen
,
Phys. Rev. B
86
,
125111
(
2012
).
37.
O.
Hod
,
C. A.
Rodríguez-Rosario
,
T.
Zelovich
, and
T.
Frauenheim
,
J. Phys. Chem. A
120
,
3278
(
2016
).
38.
U. N.
Morzan
,
F. F.
Ramírez
,
M. C.
González Lebrero
, and
D. A.
Scherlis
,
J. Chem. Phys.
146
,
044110
(
2017
).
39.
M.
Di Ventra
and
T. N.
Todorov
,
J. Phys.: Condens. Matter
16
,
8025
(
2004
).
40.
G.
Stefanucci
and
C.-O.
Almbladh
,
Europhys. Lett.
67
,
14
(
2004
).
41.
N.
Bushong
,
N.
Sai
, and
M.
Di Ventra
,
Nano Lett.
5
,
2569
(
2005
).
42.
C.-C.
Chien
,
M.
Zwolak
, and
M.
Di Ventra
,
Phys. Rev. A
85
,
041601
(
2012
).
43.
C.-C.
Chien
,
D.
Gruss
,
M.
Di Ventra
, and
M.
Zwolak
,
New J. Phys.
15
,
063026
(
2013
).
44.
C.-C.
Chien
,
M.
Di Ventra
, and
M.
Zwolak
,
Phys. Rev. A
90
,
023624
(
2014
).
45.
D.
Gruss
,
C.-C.
Chien
,
J.
Barreiro
,
M.
Di Ventra
, and
M.
Zwolak
, e-print arXiv:1610.01903 (
2016
).

Supplementary Material