Communication: Relaxation-limited electronic currents in extended reservoir simulations

“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 Wingreen^4–6 for interacting and non-interacting systems alike (or, in the case of non-interacting systems, the Landauer formula^7,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,

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 here^2,3). The density matrix, ρ, is for the $LSR$ system, and $ck†(ck)$ are the creation (annihilation) operators for state k. The Hamiltonian is

where $HS$ is for $S$ (including, potentially, many-body interactions), $HL(R)=∑k∈L(R)ℏωkck†ck$ are for the “extended reservoirs,” and $HI=∑k∈LR∑i∈S(ℏvkick†ci+ 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 states^2,3—if $HI$ is removed when $γk+≡γkfα(ωk)$ and $γk−≡γk[1−fα(ω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., sensing^9–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

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

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.² 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

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 behavior^1–3 $I1smallγ→I2→I3largeγ$⁠, 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

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

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

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

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

with $i,j∈LSR$⁠. 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.²⁷ 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 γ² when $ω̃i≠ω̃j$ and γ is small compared to the spacing $ω̃i−ω̃j$⁠. In this limit, the total current, I₁, is

where

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, I₂—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

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

That is, the relaxation must dominate all other energy scales. In practice, this means $γk≫W$⁠, 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 γ³⁰).

For just $S$⁠, the Green’s function is

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,

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, I₃, in the strong relaxation limit is then

where

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,³¹ 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),

where

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/γ². Similarly, Fig. 2(b) shows a total $LSR$ system that has no modes in the bias window, giving γ² dependence for small γ. The error inherent in numerical integration of a Lorentzian can also yield incorrect limits for small γ [Fig. 2(c)].

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

where $v0≈−3$ eV is the hopping energy between p_z 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.³³ 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×10−3$ fs⁻¹. Almost simultaneously, a γ² 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.

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 dynamics^39–41 (including of particle transport in cold-atom setups^42–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.