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 sequencing^{1–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 approach^{21–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,

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 $\rho \xaf0$ ($C0$)^{36} at some rate *γ*. We have written this expression in terms of both the single-particle density matrix, $\rho \xaf$, and the correlation matrix, $Ckl\u2009=\u2009tr\u2009[ck\u2020cl\u2009\rho ]$, since many recent studies have focused on noninteracting systems. The single-particle Hamiltonian $H\xaf$ is defined through $H\u2009=\u2009\u2211k,lH\xaflkck\u2020cl$, where $ck$ ($ck\u2020$) 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 $\rho \xaf0$ is “polarized” by different chemical potentials. This approach was applied to mean field electrons in Refs. 39 and 40 by casting $\rho \xaf0$ in a specific form, while other studies employ an alternative $\rho \xaf0$ 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

where $HS$ is the Hamiltonian for $S$, potentially including many-body interactions, $HL(R)\u2009=\u2009\u2211k\u03f5L(R)\u210f\omega kck\u2020ck$ are the “extended reservoir” Hamiltonians, and $HI\u2009=\u2009\u2211k\u03f5LR\u2211i\u03f5S(\u210fvkick\u2020ci+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.

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

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$, $\gamma k+\u2261\gamma f\alpha (\omega k)$ and $\gamma k\u2212\u2261\gamma [1\u2212f\alpha (\omega k)]$, where *f*^{α}(*ω*_{k}) is the Fermi-Dirac distribution in the $\alpha \u03f5{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 $trck\u2020cl[H,\rho ]=[H\xaf,C]kl$, yielding

The quantity $R[C]$ is the relaxation,

where $\gamma i\xb1\u2009=\u20090$ when $i\u2009\u03f5\u2009S$, $\delta k\u03f5\alpha $ is 1 when $k\u03f5\alpha $ (and zero otherwise), and *δ*_{kl} is the typical Kronecker delta.

Taking the block form

where $C\alpha ,\alpha \u2032$ are for a subset of states, i.e., in the regions $\alpha ,\alpha \u2032\u03f5{L,S,R}$, the relaxation component becomes

The “relaxed” distributions are $C0\alpha =diag[f\alpha (\omega 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 approach^{58} 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\u03f5L(R)$]. The Hamiltonian for this dressed, extended reservoir state [Fig. 1(b)] is

with $H\u2032\u2009=\u2009ck\u2020E+E\u2020ck$ and $E\u2009=\u2009\u2211l\u03f5Ek\u210f\nu lcl$. In the following, we work in the interaction picture, $E(t)\u2009=\u2009UEk\u2020EUEk$ with $UEk=exp(\u2212\u0131HEkt/\u210f)$, so that $E(t)\u2009=\u2009\u2211l\u210f\nu lclexp(\u2212\u0131\omega lt)$ [similarly, $ck(t)\u2009=\u2009ckexp(\u2212\u0131\omega kt)$].

We begin with the Born–Markov master equation^{59}

where $\rho 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

where the Hamiltonian component of the evolution will be recovered after returning from the interaction picture. The notation $[\u22c5](t)$ indicates that all operators within the brackets are in the interaction picture. The relaxation $\gamma k\xb1=(2/\u210f2)\u222b0\u221edt\u2032\u2009J\xb1(t\u2032)exp(\u2213\u0131\omega kt\u2032)$ is given in terms of the correlation functions $J+(t\u2032)=trEkE\u2020(t\u2032)E\rho Ek$ and $J\u2212(t\u2032)=trEkE(t\u2032)E\u2020\rho Ek$. An ideal Markovian environment exhibits only time-local correlations,

where a similar expression holds for $J\u2212(t\u2032)$ [but with $f\u2009\u2192\u2009(1\u2009\u2212\u2009f)$ 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(\omega )\u2261\u210f2\gamma k+/2\pi $, and $f(\omega )\u22611$ 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 $\xb1\u221e$, as adopted in other efforts.^{32}

Considering $J\u2212(t\u2032)$, we find that $J(\omega )\u2261\u210f2\gamma k\u2212/2\pi $ and $1\u2212f(\omega )\u22611$. 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

where **G**^{r(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 $\Gamma ijL(R)(\omega )=\u0131\u2211l\u03f5L(R)vjlvli[glr(\omega )\u2212gla(\omega )]$, defined in terms of the “unperturbed”—but dressed—extended reservoir state Green’s functions, $gkr(a)(\omega )=[\omega \u2212\omega k\xb1\u0131\gamma /2]\u22121$. One may also obtain the lesser Green’s function, $gk<(\omega )=\u2212fL(R)(\omega )[gkr(\omega )\u2212gka(\omega )]$. 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 formula^{19,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 solution^{65} and holds for thermal^{66–68} as well as electronic transport.^{35,69}

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

In the Markovian case, the extended reservoir state is being relaxed to the occupation $fL(R)(\omega k)$ and then broadened by *γ*. For the non-Markovian case, the *γ*-dressed state has the proper occupation $fL(R)(\omega )$. 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 $\gamma \u2009\u226a\u2009kBT/\u210f$, with temperature *T* and *k*_{B} as the Boltzmann constant (or, in terms of time scales, $\gamma \u22121\u2009\u226b\u200925$ 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, *N*_{r}. So long as the extended reservoirs are sizable, one can take $\gamma \u2009\u226a\u2009kBT/\u210f$ 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 $\gamma \u2009\u2248\u2009W/Nr$, where *W* is the bandwidth of the reservoirs, i.e., when *γ* is on the order of the mode spacing in the extended reservoirs^{35,69} (of course, inhomogeneity in the mode spacing can change this^{69}). Such behavior was recognized in earlier studies that employed $\gamma \u2009>\u2009W/Nr$.^{39,40,48}

Putting these conditions together gives $Nr\u2009\u2248\u2009\u210fW/kBT$ (hence, *N*_{r} 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 $\gamma \u2009\u226a\u2009kBT/\u210f$ with *N*_{r} 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.

## REFERENCES

These are not always true density/correlation matrices.

It is convenient to have $H\xaflk$ instead of $H\xafkl$ in *H* to give the commutator in Eq. (1) rather than one with a transpose.

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

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

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

If the integral over *ω* is taken first, we take that $\u222b0\u221edt\u2032\u2009exp[\u0131\omega kt\u2032]\delta (t\u2032)=1/2$ [i.e., $\delta (t\u2032)$ is treated as an even function]. If the integral over $t\u2032$ is performed first, we use $\u222b0\u221edt\u2032\u2009exp[\u0131(\omega k\u2212\omega )t\u2032]=\pi \delta (\omega k\u2212\omega )+\u0131P\u222b\u2212\u221e\u221ed\omega (\omega k\u2212\omega )\u22121$, for which the principal value vanishes under Markovian conditions.

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