Master equations are increasingly popular for the simulation of time-dependent electronic transport in nanoscale devices. Several recent Markovian approaches use “extended reservoirs”—explicit degrees of freedom associated with the electrodes—distinguishing them from many previous classes of master equations. Starting from a Lindblad equation, we develop a common foundation for these approaches. Due to the incorporation of explicit electrode states, these methods do not require a large bias or even “true Markovianity” of the reservoirs. Nonetheless, their predictions are only physically relevant when the Markovian relaxation is weaker than the thermal broadening and when the extended reservoirs are “sufficiently large,” in a sense that we quantify. These considerations hold despite complete positivity and respect for Pauli exclusion at any relaxation strength.

Nanoscale electronics have made inroads into a diverse range of applications, from tunneling-based DNA sequencing1–8 to high-performance microelectronics.9–18 The theoretical description of these devices is complicated by strong environmental effects, which profoundly influence electronic transport and lead to behavior beyond the static Landauer formalism. While a formally exact solution for such time-dependent transport exists, it requires the use of computationally demanding two-time Green’s functions,19,20 which are impractical for many applications. The description of sensing devices also necessitates an accounting of atomic fluctuations and unknown structural details, complicating their simulation.

One major goal is the development of a theoretical framework that can circumvent these limitations while remaining versatile enough to augment contemporary electronic structure methods. Quantum master equations and density-matrix propagation afford such an approach21–23 and encompass a diversity of well-established schemes that lie largely in the Markovian limit,24–30 with notable non-Markovian extensions.31–33 While these methods have recently come to the forefront, their conceptual history dates to the early work of Kohn and Luttinger,34 with subsequent developments that defined a device and its contacts as an open quantum system.21,22

We focus on a specific Markovian master equation,

ρ¯̇=ı[H¯,ρ¯]γ(ρ¯ρ¯0)Single Particle Density MatrixĊ=ı[H¯,C]γ(CC0)Single Particle Correlation Matrix,
(1)

in a particular context where both explicit and implicit extended reservoirs are present (introduced below).35 This master equation corresponds to a “relaxation approximation” where the system relaxes to ρ¯0 (C0)36 at some rate γ. We have written this expression in terms of both the single-particle density matrix, ρ¯, and the correlation matrix, Ckl=tr[ckclρ], since many recent studies have focused on noninteracting systems. The single-particle Hamiltonian H¯ is defined through H=k,lH¯lkckcl, where ck (ck) are the creation (annihilation) operators for the state k.37 We reserve ρ and H for the full, many-body density matrix and Hamiltonian.38 Strictly speaking, ρ, C, and H are defined for an arbitrary number of particles, allowing for fluctuations during time evolution. A current will flow when the reservoir component of ρ¯0 is “polarized” by different chemical potentials. This approach was applied to mean field electrons in Refs. 39 and 40 by casting ρ¯0 in a specific form, while other studies employ an alternative ρ¯0 that includes coherences between the device and reservoirs.41–47 

While related relaxation-type approximations have a lengthy history,21,22,24–26,34 this specific dual-reservoir setup is new and foundational to a family of promising real-time simulation methods.35,39–54 Here, we provide a rigorous justification to this setup, leading to both a well-defined domain of applicability and a connection between different variants of the formalism. More explicitly, this yields a mathematical rationale for their use in arbitrary systems (e.g., in terms of reservoir sizes and many-body interactions) and identifies relevant physical limitations, laying the foundation for future applications and implementations.

One issue with Eq. (1) is that—while it is Markovian—it is not in the standard Lindblad form.55,56 As such, it is not obviously positive and may yield both unphysical results and negative probabilities under certain conditions. For noninteracting electrons, the use of Eq. (1) has been shown to be positive for asymptotically large reservoirs.45 However, rather than starting from Eq. (1), we would like an expression that is already in Lindblad form, which will allow us to guarantee complete positivity.

We begin by examining the model depicted in Fig. 1 and analyzed in Ref. 35, where two electronic reservoirs, left (L) and right (R), drive current through a device S that contains the system of interest (for instance, a nanoscale junction and its electronic leads). The reservoir regions are finite and explicitly part of the simulation. We term these “extended reservoirs” to distinguish from the typical assumption that they are infinite and implicit.19,20 In order to have a true steady state, implicit environments EL(R) are introduced to relax L (R) to their equilibrium distributions—the notion of equilibrium is central to the use of these master equations. The Hamiltonian for this setup is

H=HS+HL+HR+HI,
(2)

where HS is the Hamiltonian for S, potentially including many-body interactions, HL(R)=kϵL(R)ωkckck are the “extended reservoir” Hamiltonians, and HI=kϵLRiϵS(vkickci+h.c.) is the interaction that couples them. The index k includes all labels (electronic state, spin, and reservoir), while ωk and vki denote the level and hopping frequencies.

FIG. 1.

Schematic of electronic transport. (a) “Extended reservoir” regions L and R drive a current I through the system S. Since the extended reservoirs are finite, external environments EL and ER act to relax them back to equilibrium, maintaining a true steady state. (b) Each extended reservoir state k is dressed by an infinite environment Ek in equilibrium (yielding non-Markovian dynamics) or with properly balanced injection, γk+, and depletion, γk, rates directing the state back to equilibrium (yielding Markovian dynamics).

FIG. 1.

Schematic of electronic transport. (a) “Extended reservoir” regions L and R drive a current I through the system S. Since the extended reservoirs are finite, external environments EL and ER act to relax them back to equilibrium, maintaining a true steady state. (b) Each extended reservoir state k is dressed by an infinite environment Ek in equilibrium (yielding non-Markovian dynamics) or with properly balanced injection, γk+, and depletion, γk, rates directing the state back to equilibrium (yielding Markovian dynamics).

Close modal

The LSR system is open. Under the influence of EL(R), its dynamics is given by the Markovian master equation

ρ̇=ı[H,ρ]+kγk+ckρck12ckck,ρ+kγkckρck12ckck,ρ
(3)

for the full, many-body density matrix ρ. The first term on the right is the Hamiltonian evolution of ρ under H and the second (third) term reflects particle injection (depletion) into the state k at a rate γk+ (γk). To ensure that the reservoirs relax to equilibrium—a pseudo-equilibrium, as we will see—in the absence of S, γk+γfα(ωk) and γkγ[1fα(ωk)], where fα(ωk) is the Fermi-Dirac distribution in the αϵ{L,R} reservoir and with γ nonzero only for reservoir states. We assume a general case where each reservoir may be at a different chemical potential or temperature. This specific master equation has appeared in previous efforts.35,48–54 In particular, Ref. 35 derives the closed-form solution for both the interacting and noninteracting cases, as well as those for the related non-Markovian problem.

To connect Eq. (3) to noninteracting approaches,39–47 we first differentiate the single-particle correlation matrix C, employ Eq. (3), and use that trckcl[H,ρ]=[H¯,C]kl, yielding

Ċ=ı[H¯,C]/+R[C].
(4)

The quantity R[C] is the relaxation,

(R[C])kl=γk+δklCkl2γk++γk+γl++γl=γfkαδklδkϵLRCkl2(δkϵLR+δlϵLR),
(5)

where γi±=0 when iϵS, δkϵα is 1 when kϵα (and zero otherwise), and δkl is the typical Kronecker delta.

Taking the block form

C=CL,LCL,SCL,RCS,LCS,SCS,RCR,LCR,SCR,R,
(6)

where Cα,α are for a subset of states, i.e., in the regions α,αϵ{L,S,R}, the relaxation component becomes

R[C]=γ(CL,LC0L)12CL,SCL,R12CS,L012CS,RCR,L12CR,S(CR,RC0R).
(7)

The “relaxed” distributions are C0α=diag[fα(ωk)]. For simplicity, Eqs. (5) and (7) are written in the single-particle eigenbasis of the decoupled L, S, and R regions.

The equation of motion defined by Eqs. (4) and (7) is exactly that of Refs. 41–47. Since the starting expression is in Lindblad form, this demonstrates that these prior approaches use a completely positive, trace-preserving master equation for the single-particle matrices. Our derivation shows that these properties always hold, including for finite reservoirs as well as those that are asymptotically large.57 Moreover, by virtue of the use of creation/annihilation operators in Eq. (3), Pauli exclusion is obeyed even though the particle number is not conserved. Furthermore, if we make an approximation where the off-diagonal coherences are negligible, the phenomenological expression in Eq. (1) using C0 from Refs. 39 and 40 is recovered. This is not, however, guaranteed to preserve positivity.

Too V, or not too V? In the preceding discussion, we adopted a Lindblad equation from the outset. To take a more foundational perspective, we can use the Born–Markov approach58 to derive this equation. In doing so, we see that there must be two implicit reservoirs with a high voltage between them. While this might appear to nullify the use of Eq. (3), we demonstrate that our approach is physically applicable. The following derivation is presented for a single extended reservoir state, which is sufficient for completeness as these states are separately relaxed.

Each extended reservoir state k is connected to an implicit reservoir, the environment Ek [EL(R) from Fig. 1(a) are composed of all Ek for kϵL(R)]. The Hamiltonian for this dressed, extended reservoir state [Fig. 1(b)] is

Hk=ωkckck+lϵEkωlclcl+lϵEkνlclck+ckclωkckck+lϵEkωlclcl+H,
(8)

with H=ckE+Eck and E=lϵEkνlcl. In the following, we work in the interaction picture, E(t)=UEkEUEk with UEk=exp(ıHEkt/), so that E(t)=lνlclexp(ıωlt) [similarly, ck(t)=ckexp(ıωkt)].

We begin with the Born–Markov master equation59 

ρ̇k(t)=120dttrEkH(t),H(tt),ρk(t)ρEk,
(9)

where ρEk denotes the initial state of Ek. This approximation requires a weak coupling between k and Ek (essentially, second order perturbation theory) and the assumption of an uncorrelated, time-local composite state of the system and environment. The latter needs justification, which Eq. (11) will provide. Expanding the commutators and taking the trace in Eq. (9) give

ρ̇k(t)=(γk/2)ckρkckckckρk(t)+h.c.+(γk+/2)ckρkckckckρk(t)+h.c.,
(10)

where the Hamiltonian component of the evolution will be recovered after returning from the interaction picture. The notation [](t) indicates that all operators within the brackets are in the interaction picture. The relaxation γk±=(2/2)0dtJ±(t)exp(ıωkt) is given in terms of the correlation functions J+(t)=trEkE(t)EρEk and J(t)=trEkE(t)EρEk. An ideal Markovian environment exhibits only time-local correlations,

J+(t)=dωJ(ω)f(ω)eıωtδ(t)2γk+,
(11)

where a similar expression holds for J(t) [but with f(1f) and a change of sign in the exponent]. To obtain a δ-function requires that the product J(ω)f(ω) is constant, which can only be physically satisfied if the spectral function is flat, J(ω)2γk+/2π, and f(ω)1 for all ω.60 In other words, the implicit reservoir is completely full. This is evident from the definition of a Markovian reservoir—an environment that couples to the system equally at all frequency scales. The presence of the Fermi level breaks this symmetry. Thus, this level must lie at ±, as adopted in other efforts.32 

Considering J(t), we find that J(ω)2γk/2π and 1f(ω)1. This implies that two distinct sets of states are required to obtain Eq. (3): In one set, the states are completely empty, acting only to deplete particles from k [the first line of Eq. (10)]. In the other set, the states are completely full and thus they only inject particles into k [the second line of Eq. (10)]. References 61–64 address this process when the implicit reservoirs connect directly to what we would call S, concluding that the equation of motion corresponds to a high bias V. In our approach, however, the implicit reservoirs are not connected directly to S, but rather indirectly through an intermediary—the extended reservoirs. The bias is thus wrapped into the simulation and we do not require a high value of V. We will see this more explicitly below where we show that, when quantifying the errors of the Markovian equation (3) for steady states, nowhere does the bias show up, but rather only the temperature and extended reservoir size. The requirement for Markovianity may also be relaxed since explicit reservoir states retain a memory up to time 1/γ of the dynamics.

True limitations. Equation (3) is completely positive, is trace preserving, respects Pauli exclusion, and does not require a high V for its use. Nonetheless, these properties are not sufficient to ensure physically meaningful behavior. To quantify this statement, we make use of the exact, closed-form solution of Eq. (3) and its non-Markovian counterpart, both given in Ref. 35. The latter also uses the model Hamiltonian from Eq. (8); however, it does not require two distinct full and empty components of Ek (nor a weak coupling between k and Ek, a flat spectral function, or the wide-band limit). Rather, Ek need only be infinite and in an equilibrium state described by the Fermi-Dirac distribution fL(R).

For this non-Markovian case, the current is

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

where Gr(a)(ω) are the—potentially many-body—retarded (advanced) Green’s functions for S [see Ref. 35 for the closed-form solution to the Markovian case, Eq. (3)]. The spectral densities of the couplings between S and L(R) are ΓijL(R)(ω)=ılϵL(R)vjlvli[glr(ω)gla(ω)], defined in terms of the “unperturbed”—but dressed—extended reservoir state Green’s functions, gkr(a)(ω)=[ωωk±ıγ/2]1. One may also obtain the lesser Green’s function, gk<(ω)=fL(R)(ω)[gkr(ω)gka(ω)]. For simplicity, we only address the wide-band limit (the more general case is in Ref. 35). Notable in these expressions is the term γ, which accounts for relaxation and is key to subsequent discussion. These results diverge from the Meir–Wingreen formula19,20 in the use of extended reservoirs as a finite-sized intermediary with relaxation.

While relaxation processes occur in real materials, these are not necessarily of the form in Eq. (8). Moreover, by taking Markovian relaxation as a further approximation, we cannot conclude that Eq. (3) will be physical for a given γ. We thus interpret the relaxation as a control parameter, diligently chosen to obtain meaningful results from Eq. (3). Numerical simulations [Fig. 2(a)] illustrate how the current in the full non-Markovian model behaves versus γ. There are three distinct regimes: a regime linear in γ, a plateau regime, and a 1/γ regime. In the intermediate plateau regime, the “intrinsic conductance” of the setup determines the current (for non-interacting systems, this would be the Landauer current). We note that the physics of the turnover versus γ is analogous to Kramers’ turnover for reaction rates in solution65 and holds for thermal66–68 as well as electronic transport.35,69

FIG. 2.

(a) Steady-state current, I, as a function of the reservoir relaxation parameter, γ, Eq. (12), with an applied bias of VL(R)=±v0/4 and an increasingly large number of explicit reservoir sites Nrϵ{32,128,512}. Here, we use a 1D model with an asymmetric shift: wkL(R)=2v0coskπ/(Nr+1)±v0/10 and vki=v02/(Nr+1)sinkπ/(Nr+1) (with kϵ{1,,Nr} for both L and R), kBT=v0/40, and a single site at zero energy in S. In the limit of Nr (black line), the plateau extends to γ = 0, yielding the Landauer/Meir-Wingreen current exactly (red dashed line). The dotted line demarcates where the Markovian master equation is valid, which can be reached numerically by increasing the number of reservoir sites. [(b) and (c)] Anomalous zero-bias (VL(R)=0) currents versus time t for (b) a weak relaxation (γ=0.1v0) and (c) a strong relaxation (γ=v0). The red dashed line indicates their limit as t. In the non-Markovian case, however, there is zero steady-state current. (d) Anomalous broadening of gk< in L for strong relaxation (γ=v0) as shown by Eq. (13). The solid blue line is the correct “broadening before occupying,” while the dashed green line is the Markovian counterpart, which extends the occupation well beyond the Fermi level (black dotted line).

FIG. 2.

(a) Steady-state current, I, as a function of the reservoir relaxation parameter, γ, Eq. (12), with an applied bias of VL(R)=±v0/4 and an increasingly large number of explicit reservoir sites Nrϵ{32,128,512}. Here, we use a 1D model with an asymmetric shift: wkL(R)=2v0coskπ/(Nr+1)±v0/10 and vki=v02/(Nr+1)sinkπ/(Nr+1) (with kϵ{1,,Nr} for both L and R), kBT=v0/40, and a single site at zero energy in S. In the limit of Nr (black line), the plateau extends to γ = 0, yielding the Landauer/Meir-Wingreen current exactly (red dashed line). The dotted line demarcates where the Markovian master equation is valid, which can be reached numerically by increasing the number of reservoir sites. [(b) and (c)] Anomalous zero-bias (VL(R)=0) currents versus time t for (b) a weak relaxation (γ=0.1v0) and (c) a strong relaxation (γ=v0). The red dashed line indicates their limit as t. In the non-Markovian case, however, there is zero steady-state current. (d) Anomalous broadening of gk< in L for strong relaxation (γ=v0) as shown by Eq. (13). The solid blue line is the correct “broadening before occupying,” while the dashed green line is the Markovian counterpart, which extends the occupation well beyond the Fermi level (black dotted line).

Close modal

When simulating Markovian dynamics using Eq. (3), instead of the non-Markovian problem, the three regimes are still present, but a large γ can result in non-zero currents at zero bias [Figs. 2(b) and 2(c)]. The origin of these currents is due to improper occupation of the extended reservoir levels. Calculating the real-time correlation functions from Eq. (3), we see that the advanced and retarded Green’s functions have a functional form that is identical to the non-Markovian case.35 In fact, the difference between the Markovian and non-Markovian limits is encapsulated by the replacement

gk<(ω)=ıγfL(R)(ω)(ωωk)2+γ2/4Non-MarkovianıγfL(R)(ωk)(ωωk)2+γ2/4Markovian.
(13)

In the Markovian case, the extended reservoir state is being relaxed to the occupation fL(R)(ωk) and then broadened by γ. For the non-Markovian case, the γ-dressed state has the proper occupation fL(R)(ω). In other words, the non-Markovian case dresses and then occupies and the Markovian case occupies and then dresses. Figure 2(d) demonstrates that the Markovian limit gives an additional occupancy above the Fermi level, leading to zero-bias currents under certain conditions. Thus, the Markovian master equation, Eq. (3) or Eq. (4), can yield unphysical behavior despite the fact that it is always completely positive and obeys Pauli exclusion. Another way to state the origin of this behavior is that the Markovian master equation relaxes the extended reservoir into pseudo-equilibrium—an equilibrium defined in terms of isolated extended reservoir states rather than those in the presence of the environment that provides the relaxation.

We can define precisely when the replacement in Eq. (13) yields a reasonable approximation, thus providing a satisfying quantification of the validity of Eq. (3): So long as the γ-induced broadening is less than the thermal broadening γkBT/, with temperature T and kB as the Boltzmann constant (or, in terms of time scales, γ125 fs at room temperature), the Markovian limit accurately gives the steady-state solutions.35 This is independent of any details of S—it may be interacting, may be non-interacting, may have electron-phonon coupling, etc. The validity hinges on the replacement made in Eq. (13), which, in turn, relies only on the fact that the reservoir states are non-interacting. This is generally a good approximation. While we do not quantify it here, the dynamics of interest will be correctly captured by Eq. (3) as long as they are faster than the relaxation, as the latter only cuts off behavior after a time γ−1. These limits must be carefully enforced in order to ensure physically meaningful dynamics and steady-state currents.

The considerations above lead to a natural estimate for the required number of extended reservoir states, Nr. So long as the extended reservoirs are sizable, one can take γkBT/ and be within the plateau regime. A simple estimate is given by the turnover point between the linear and plateau regions. This generically occurs when γW/Nr, where W is the bandwidth of the reservoirs, i.e., when γ is on the order of the mode spacing in the extended reservoirs35,69 (of course, inhomogeneity in the mode spacing can change this69). Such behavior was recognized in earlier studies that employed γ>W/Nr.39,40,48

Putting these conditions together gives NrW/kBT (hence, Nr should be very large at low temperature). A less stringent condition on γ would only require the current be on the plateau, which often extends to relatively large values of γ. It is a mistake, however, to conclude that an arbitrary, large value of γ will be acceptable, even if this condition holds. A sufficiently large γ will improperly occupy high energy states, which—when asymmetric reservoirs are present—will give rise to unphysical zero-bias currents.70 Even though γ plays the role of relaxation in the extended reservoirs, a small (or, in a sense, “intermediary”) value is still necessary.

We see that Markovian master equations can be a powerful tool for the simulation of electronic transport. Furthermore, the various approaches employed in the literature can be unified as equivalent expressions of Eq. (3) or some approximations thereof. These extended reservoir-based Markovian master equations do not require a large bias or even Markovianity. The true limit of the Markovian limit, Eq. (3), is the requirement that γkBT/ with Nr that is large enough to accommodate this slow relaxation and still yield the intrinsic conductance.

J.E.E. and D.G. acknowledge support under the Cooperative Research Agreement between the University of Maryland and the National Institute for Standards and Technology Center for Nanoscale Science and Technology, Award No. 70NANB14H209, through the University of Maryland.

1.
M.
Zwolak
and
M.
Di Ventra
,
Nano Lett.
5
,
421
(
2005
).
2.
J.
Lagerqvist
,
M.
Zwolak
, and
M.
Di Ventra
,
Nano Lett.
6
,
779
(
2006
).
3.
M.
Zwolak
and
M.
Di Ventra
,
Rev. Mod. Phys.
80
,
141
(
2008
).
4.
D.
Branton
,
D. W.
Deamer
,
A.
Marziali
,
H.
Bayley
,
S. A.
Benner
,
T.
Butler
,
M.
Di Ventra
,
S.
Garaj
,
A.
Hibbs
,
X.
Huang
 et al,
Nat. Biotechnol.
26
,
1146
(
2008
).
5.
S.
Chang
,
S.
Huang
,
J.
He
,
F.
Liang
,
P.
Zhang
,
S.
Li
,
X.
Chen
,
O.
Sankey
, and
S.
Lindsay
,
Nano Lett.
10
,
1070
(
2010
).
6.
M.
Tsutsui
,
M.
Taniguchi
,
K.
Yokota
, and
T.
Kawai
,
Nat. Nanotechnol.
5
,
286
(
2010
).
7.
S.
Huang
,
J.
He
,
S.
Chang
,
P.
Zhang
,
F.
Liang
,
S.
Li
,
M.
Tuchband
,
A.
Fuhrmann
,
R.
Ros
, and
S.
Lindsay
,
Nat. Nanotechnol.
5
,
868
(
2010
).
8.
T.
Ohshiro
,
K.
Matsubara
,
M.
Tsutsui
,
M.
Furuhashi
,
M.
Taniguchi
, and
T.
Kawai
,
Sci. Rep.
2
,
501
(
2012
).
9.
M.
Büttiker
and
S. E.
Nigg
,
Nanotechnology
18
,
044029
(
2007
).
10.
J.
Gabelli
,
G.
Fève
,
J.-M.
Berroir
,
B.
Plaçais
,
A.
Cavanna
,
B.
Etienne
,
Y.
Jin
, and
D.
Glattli
,
Science
313
,
499
(
2006
).
11.
G. P.
Lansbergen
,
R.
Rahman
,
C. J.
Wellard
,
I.
Woo
,
J.
Caro
,
N.
Collaert
,
S.
Biesemans
,
G.
Klimeck
,
L. C. L.
Hollenberg
, and
S.
Rogge
,
Nat. Phys.
4
,
656
(
2008
).
12.
A.
Morello
,
J. J.
Pla
,
F. A.
Zwanenburg
,
K. W.
Chan
,
K. Y.
Tan
,
H.
Huebl
,
M.
Möttönen
,
C. D.
Nugroho
,
C.
Yang
,
J. A.
van Donkelaar
,
A. D. C.
Alves
,
D. N.
Jamieson
,
C. C.
Escott
,
L. C. L.
Hollenberg
,
R. G.
Clark
, and
A. S.
Dzurak
,
Nature
467
,
687
(
2010
).
13.
K. Y.
Tan
,
K. W.
Chan
,
M.
Möttönen
,
A.
Morello
,
C.
Yang
,
J.
van Donkelaar
,
A.
Alves
,
J.-M.
Pirkkalainen
,
D. J.
Jamieson
,
R. G.
Clark
, and
A. S.
Dzurak
,
Nano Lett.
10
,
11
(
2010
).
14.
M.
Pierre
,
R.
Wacquez
,
X.
Jehl
,
M.
Sanquer
,
M.
Vinet
, and
O.
Cueto
,
Nat. Nanotechnol.
5
,
133
(
2010
).
15.
M.
Fuechsle
,
J. A.
Miwa
,
S.
Mahapatra
,
H.
Ryu
,
S.
Lee
,
O.
Warschkow
,
L. C. L.
Hollenberg
,
G.
Klimeck
, and
M. Y.
Simmons
,
Nat. Nanotechnol.
7
,
242
(
2012
).
16.
M. M.
Shulaker
,
G.
Hills
,
N.
Patil
,
H.
Wei
,
H.-Y.
Chen
,
H.-S. P.
Wong
, and
S.
Mitra
,
Nature
501
,
526
(
2013
).
17.
M. L.
Perrin
,
E.
Burzurí
, and
H. S. J.
van der Zant
,
Chem. Soc. Rev.
44
,
902
(
2015
).
18.
J.
Trasobares
,
D.
Vuillaume
,
D.
Théron
, and
N.
Clément
,
Nat. Commun.
7
,
12850
(
2016
).
19.
Y.
Meir
and
N. S.
Wingreen
,
Phys. Rev. Lett.
68
,
2512
(
1992
).
20.
A.-P.
Jauho
,
N. S.
Wingreen
, and
Y.
Meir
,
Phys. Rev. B
50
,
5528
(
1994
).
21.
W. R.
Frensley
,
J. Vac. Sci. Technol., B: Microelectron. Nanometer Struct.
3
,
1261
(
1985
).
22.
W. R.
Frensley
,
Rev. Mod. Phys.
62
,
745
(
1990
).
23.
I.
Knezevic
and
B.
Novakovic
,
J. Comput. Electron.
12
,
363
(
2013
).
24.
M. V.
Fischetti
,
J. Appl. Phys.
83
,
270
(
1998
).
25.
M. V.
Fischetti
,
Phys. Rev. B
59
,
4901
(
1999
).
26.
H.
Mizuta
and
C. J.
Goodings
,
J. Phys.: Condens. Matter
3
,
3739
(
1991
).
27.
I.
Knezevic
,
Phys. Rev. B
77
,
125301
(
2008
).
28.
G.
Schaller
and
T.
Brandes
,
Phys. Rev. A
78
,
022106
(
2008
).
29.
R.
Rosati
,
R. C.
Iotti
,
F.
Dolcini
, and
F.
Rossi
,
Phys. Rev. B
90
,
125140
(
2014
).
30.
R.
Rosati
,
D. E.
Reiter
, and
T.
Kuhn
,
Phys. Rev. B
95
,
165302
(
2017
).
31.
G.
Schaller
,
P.
Zedler
, and
T.
Brandes
,
Phys. Rev. A
79
,
032110
(
2009
).
32.
P.
Zedler
,
G.
Schaller
,
G.
Kiesslich
,
C.
Emary
, and
T.
Brandes
,
Phys. Rev. B
80
,
045309
(
2009
).
33.
E.
Vaz
and
J.
Kyriakidis
,
Phys. Rev. B
81
,
085315
(
2010
).
34.
W.
Kohn
and
J. M.
Luttinger
,
Phys. Rev.
108
,
590
(
1957
).
35.
D.
Gruss
,
K.
Velizhanin
, and
M.
Zwolak
,
Sci. Rep.
6
,
24514
(
2016
).
36.

These are not always true density/correlation matrices.

37.

It is convenient to have H¯lk instead of H¯kl in H to give the commutator in Eq. (1) rather than one with a transpose.

38.

The quantities ρ¯ and C give all the information required to reconstruct the full density matrix for noninteracting systems. If, though, one takes ρ¯ as trace one, then the total number of particles is also required to reconstruct ρ, but C always contains the number of particles.

39.
C. G.
Sánchez
,
M.
Stamenova
,
S.
Sanvito
,
D.
Bowler
,
A. P.
Horsfield
, and
T. N.
Todorov
,
J. Chem. Phys.
124
,
214708
(
2006
).
40.
J. E.
Subotnik
,
T.
Hansen
,
M. A.
Ratner
, and
A.
Nitzan
,
J. Chem. Phys.
130
,
144105
(
2009
).
41.
L.
Chen
,
T.
Hansen
, and
I.
Franco
,
J. Phys. Chem. C
118
,
20009
(
2014
).
42.
T.
Zelovich
,
L.
Kronik
, and
O.
Hod
,
J. Chem. Theory Comput.
10
,
2927
(
2014
).
43.
T.
Zelovich
,
L.
Kronik
, and
O.
Hod
,
J. Chem. Theory Comput.
11
,
4861
(
2015
).
44.
T.
Zelovich
,
L.
Kronik
, and
O.
Hod
,
J. Phys. Chem. C
120
,
15052
(
2016
).
45.
O.
Hod
,
C. A.
Rodríguez-Rosario
,
T.
Zelovich
, and
T.
Frauenheim
,
J. Phys. Chem. A
120
,
3278
(
2016
).
46.
U. N.
Morzan
,
F. F.
Ramírez
,
M. C.
González Lebrero
, and
D. A.
Scherlis
,
J. Chem. Phys.
146
,
044110
(
2017
).
47.
T.
Zelovich
,
T.
Hansen
,
Z.-F.
Liu
,
J. B.
Neaton
,
L.
Kronik
, and
O.
Hod
,
J. Chem. Phys.
146
,
092331
(
2017
).
48.
A. A.
Dzhioev
and
D. S.
Kosov
,
J. Chem. Phys.
134
,
044121
(
2011
).
49.
S.
Ajisaka
,
F.
Barra
,
C.
Mejía-Monasterio
, and
T.
Prosen
,
Phys. Rev. B
86
,
125111
(
2012
).
50.
S.
Ajisaka
and
F.
Barra
,
Phys. Rev. B
87
,
195114
(
2013
).
51.
D.
Dast
,
D.
Haag
,
H.
Cartarius
, and
G.
Wunner
,
Phys. Rev. A
90
,
052120
(
2014
).
52.
S.
Ajisaka
,
B.
Žunkovič
, and
Y.
Dubi
,
Sci. Rep.
5
,
8312
(
2015
).
53.
R.
Mahajan
,
C. D.
Freeman
,
S.
Mumford
,
N.
Tubman
, and
B.
Swingle
, e-print arXiv:1608.05074 (
2016
).
54.
C.
Zanoci
and
B. G.
Swingle
, e-print arXiv:1612.04840 (
2016
).
55.
G.
Lindblad
,
Commun. Math. Phys.
48
,
119
(
1976
).
56.
V.
Gorini
,
A.
Kossakowski
, and
E. C. G.
Sudarshan
,
J. Math. Phys.
17
,
821
(
1976
).
57.

The latter was shown in Ref. 45, where they call “leads” what we call “extended reservoirs.”

58.
H.-P.
Breuer
and
F.
Petruccione
,
The Theory of Open Quantum Systems
(
Oxford University Press
,
2002
).
59.

In most scenarios, including ours, the first-order term vanishes, trEk[H(t),ρ(0)]=0.

60.

If the integral over ω is taken first, we take that 0dtexp[ıωkt]δ(t)=1/2 [i.e., δ(t) is treated as an even function]. If the integral over t is performed first, we use 0dtexp[ı(ωkω)t]=πδ(ωkω)+ıPdω(ωkω)1, for which the principal value vanishes under Markovian conditions.

61.
S.
Gurvitz
and
Y. S.
Prager
,
Phys. Rev. B
53
,
15932
(
1996
).
62.
63.
H.
Sprekeler
,
G.
Kießlich
,
A.
Wacker
, and
E.
Schöll
,
Phys. Rev. B
69
,
125328
(
2004
).
64.
U.
Harbola
,
M.
Esposito
, and
S.
Mukamel
,
Phys. Rev. B
74
,
235309
(
2006
).
66.
K. A.
Velizhanin
,
S.
Sahu
,
C.-C.
Chien
,
Y.
Dubi
, and
M.
Zwolak
,
Sci. Rep.
5
,
17506
(
2015
).
67.
C.-C.
Chien
,
S.
Kouachi
,
K. A.
Velizhanin
,
Y.
Dubi
, and
M.
Zwolak
,
Phys. Rev. E
95
,
012137
(
2017
).
68.
C.-C.
Chien
,
K. A.
Velizhanin
,
Y.
Dubi
,
B. R.
Ilic
, and
M.
Zwolak
, e-print arXiv:1707.06669 (
2017
).
69.
D.
Gruss
,
A.
Smolyanitsky
, and
M.
Zwolak
,
J. Chem. Phys.
(to be published); e-print arXiv:1707.06650 (
2017
).
70.

When the reservoirs are symmetric, the forward and backward anomalous zero-bias currents cancel each other. However, this will still give anomalous correlations.