Non-Markovian electron dynamics in nanostructures coupled to dissipative contacts
aa r X i v : . [ c ond - m a t . m e s - h a ll ] M a y Non-Markovian electron dynamics in nanostructures coupled to dissipative contacts
B. Novakovic ∗ and I. Knezevic † Department of Electrical and Computer Engineering, University of Wisconsin – Madison, Madison, WI 53706, USA
In quasiballistic semiconductor nanostructures, carrier exchange between the active region and dissipativecontacts is the mechanism that governs relaxation. In this paper, we present a theoretical treatment of transientquantum transport in quasiballistic semiconductor nanostructures, which is based on the open system theoryand valid on timescales much longer than the characteristic relaxation time in the contacts. The approach relieson a model interaction between the current-limiting active region and the contacts, given in the scattering-state basis. We derive a non-Markovian master equation for the irreversible evolution of the active region’smany-body statistical operator by coarse-graining the exact dynamical map over the contact relaxation time. Inorder to obtain the response quantities of a nanostructure under bias, such as the potential and the charge andcurrent densities, the non-Markovian master equation must be solved numerically together with the Schr¨odinger,Poisson, and continuity equations. We discuss how to numerically solve this coupled system of equations andillustrate the approach on the example of a silicon nin diode.
I. INTRODUCTION
In nanoscale, quasiballistic electronic systems under bias,the process of relaxation towards a nonequilibrium steadystate cannot be attributed to scattering, because these struc-tures are small compared to the carrier mean free path [1,2]. Rather, the active region of a nanostructure is an openquantum-mechanical system that exchanges particles and in-formation with the dissipative reservoirs of charge, usually re-ferred to as contacts [3, 4]. While qualitatively clear, a quanti-tative description of the irreversible evolution of the electronicsystem in this regime, where dissipation in the contacts cou-pled with the carrier exchange between the active region andcontacts is the mechanism governing relaxation, is very chal-lenging [5–7].In this paper, we present a theoretical treatment of thetransient-regime evolution of the electronic system in a two-terminal ballistic nanostructure coupled to dissipative contactsand illustrate it on the example of a semiconductor nin diode.The approach is rooted in the open system theory [8, 9]. Westart from the closed-system, Hamiltonian dynamics of themany-body statistical operator for the ballistic active regionand the dissipative contacts together, with a model interactiondescribing the injection of electrons into the active region.The model interaction Hamiltonian differs from those typi-cally employed [10, 11]: it is specifically constructed to con-serve current during the process of carrier injection from/intothe contacts, and its matrix elements are readily calculatedfrom the single-particle transmission problem for structureswith and without resonances alike (Sec. II). As is commonlydone, we trace out the contact degrees of freedom and obtainthe exact non-Markovian dynamical map that describes theevolution of the active region’s statistical operator. However,while exact, this map is not useful in practical calculations. Inorder to obtain a tractable theoretical approach, we employ thefact that relaxation in the contacts of a nanostructure typicallyoccurs on the shortest timescales in the whole system. We as- ∗ [email protected] † [email protected] sume that the contacts are highly doped, so the fastest scatter-ing mechanism is electron-electron scattering [12–14]. Withinthe momentum relaxation time, the contacts adjust themselvesto the new level of current flowing through the structure. Themomentum relaxation time is virtually instantaneous from thestandpoint of the nanostructure as a whole; if we are not tolook into the microscopic details of relaxation in the contacts,but want to include their effect on the overall evolution ofthe nanostructure, the momentum relaxation time can be con-sidered the shortest meaningfully resolvable time. Therefore,we coarse-grain the evolution over the contact momentum re-laxation time and obtain a dynamical map that is piecewiseMarkovian but globally a non-Markovian, completely posi-tive map (Sec. III). We present a numerical algorithm for thecalculation of relevant response quantities such as the chargedensity, potential, and current density based on the presentedmodel, and illustrate the approach with a calculation of theresponse of a realistic semiconductor nin diode in Sec. IV. II. INTERACTION BETWEEN THE ACTIVE REGIONAND THE CONTACTS
It has been well-established that the active region of ananostructure is an open quantum-mechanical system [15,16]. Usually, the effect of openness is treated through openboundary conditions; examples of such treatment are the ex-plicit source terms in the density matrix [17] or Wigner func-tion formalisms [15, 18]. Alternatively, a dynamical quantityis ascribed to the coupling between the active region and thecontact: in the popular tight-binding variant of the nonequi-librium Green’s function formalism, pioneered by Datta [19],the active region-contact coupling is described through a spe-cial self-energy term. In the Meir-Wingreen [10, 11] approachand its derivatives, one employs a coupling Hamiltonian be-tween the contacts and the active region, but no general recipeis available for the derivation of the Hamiltonian matrix el-ements. Also, this approach has so far been applied onlywhen the active region supports a small number of discretestates, so the model has little practical value for structureswith no resonances, such as an nin diode, or to account forthe continuum states in structures with mixed spectrum, such qV active regionleft contact right contact W ikx-k e x L x R r - ikx e k Y xikk et '' E F ∆ k c k ,L d † k FIG. 1. Schematic of the coupling between the active regionof a generic two-terminal nanostructure and the contacts. In thecase of ballistic injection through the open boundaries, a forward-propagating state Ψ k is coupled with the state exp( ikx ) in the leftcontact via a hopping model interaction (2). After [20]. as a double-barrier tunneling structure (also known as theresonant-tunneling diode). We present an alternative interac-tion Hamiltonian that does nor require that a structure a prioripossesses resonances, and whose matrix elements are straight-forwardly derived from the single particle transmission prob-lem.Consider a generic two-terminal nanostructure under bias(Fig. 1). For every energy E k above the bottom of the leftcontact, the active region’s single particle Hamiltonian hastwo eigenfunctions, a forward ( Ψ k ) and a backward ( Ψ − k )propagating state, that can be found by (in general numeri-cally) solving the single-particle Schr¨odinger equation for agiven potential profile in the active region. Associated with Ψ k ( Ψ − k ) in the active region are the creation and destruc-tion operators d † k and d k ( d †− k and d − k ), so the active regionmany-body Hamiltonian is H S = X k> ω k ( d † k d k + d †− k d − k ) . (1)Spin is disregarded, and ω k = E k / ~ . In the case of ballisticinjection through the open boundaries, each state Ψ k is natu-rally coupled with the injected states exp( ikx ) from the leftcontact. For Ψ − k , the coupling is with exp( − ik ′ x ) from theright contact ( k ′ − meV / ~ = k = 2 m E k / ~ , where V isthe applied bias). To model this coupling via a hopping-typeinteraction, we can write quite generally (see Fig. 1) H int = X k> ∆ k d † k c k,L + ∆ − k d †− k c − k ′ ,L + h.c. (2) c † k,L ( c k,L ) and c †− k ′ ,R ( c − k ′ ,R ) create (destroy) an electronwith a wavevector k in the left and − k ′ in the right contact,respectively. The hopping coefficients ∆ k and ∆ − k are pro-portional to the current carried by each mode, i.e. ∆ k = I k e T k , (3) where T k is the transmission coefficient of mode k . ∆ k canbe written in terms of the scattering-state injection amplitude[21]. III. THE TRANSPORT MASTER EQUATION
In general, the dynamics of a nanostructure’s active regionis non-unitary and non-Markovian (i.e., memory effects areimportant, meaning that the system remembers how it got to acertain state and its future direction of evolution depends notonly on the state it is currently in, but also on how it got to thatstate to begin with). A non-Markovian, non-unitary map thatwould describe the active region in a ballistic nanostructurein the presence of contacts can be derived by tracing our thecontact degrees of freedom from the unitary evolution of theclosed ”active region and contacts” system. A general formof the non-Markovian evolution of the active region statisticaloperator ρ is given by ρ ( t ) = W ( t, ρ (0) , where map W isof the form W ( t,
0) = T c exp (cid:18)Z t K ( t ′ ) dt ′ (cid:19) . (4)Here, K ( t ) is the generator of the map W ( t, . In general,it is impossible to obtain W ( t, exactly. If one is interestedin retaining the non-Markovian nature of (4), typically an ex-pansion up to the second or fourth order in the interaction isundertaken [9]. On the other hand, a Markovian approxima-tion to the exact dynamics can quite generally be obtained inthe weak-coupling limit. This limit has been used previouslyby several authors [22, 23] to derive Markovian rate equationsfor tunneling structures in the resonant-level model, althoughthe weak-coupling approximation is not generally applicableto nanostructures [22].However, here we point out the Markovian approximationto the long-time evolution of nanostructures can be justifiedmore broadly, by employing the approximation of a memory-less environment for the contacts. Basically, electron-electronscattering in the highly doped contacts of semiconductor de-vices ensures that the carrier distribution snaps into a driftedFermi-Dirac distribution [12] within the energy-relaxationtime τ ≈ − femtoseconds [13, 14] (the actual valuedepends on the doping density and temperature). This time isvery short with respect to the typical response times of thesedevices, which is on the timescales of τ AR ≈ − ps (”AR”stands for the active region), so on these timescales contactscan be considered memoryless . For low-dimensional nanos-tructures, fabricated on a high-mobility two-dimensional elec-tron gas (2DEG) and operating at low temperatures, phononsare frozen so the energy relaxation in the contacts is also gov-erned by the inelastic electron-electron scattering [24, 25].The ratio τ /τ AR is not as small as in devices, but is still lessthan unity.To practically obtain the Markovian approximation due toan environment that loses memory after a time τ , we use thecoarse-graining procedure: we can partition the time axis intointervals of length τ , t n = nτ , so the environment interactswith the system in exactly the same way during each interval [ t n , t n +1 ] [26], dρ S dt ≈ ρ S,n +1 − ρ S,n τ = K τ ρ S,n , (5)where K τ = R τ K ( t ′ ) dt ′ τ = R tn +1 tn K ( t ′ ) dt ′ τ is the averaged valueof the map’s generator over any interval [ t n , t n +1 ] ( K is resetat each t n ). If the coarse-graining time τ is short enough,then the short-time expansion of K can be used to performthe coarse-graining [20], so we finally arrive at the desiredMarkovian kinetic equation dρ S ( t ) dt = ( − i L eff − Λ τ ) ρ S ( t ) . (6)where L eff = [ H S + hH int i , . . . ] = L S + [ hH int i , . . . ] isan effective system Liouvillian, containing the noninteracting-system Liouvillian L S and a correction due to the interaction[ h . . . i = Tr E [ ρ E (0) . . . ] denotes the partial average with re-spect to the initial environmental state ρ E (0) ]. The matrixelements of superoperator Λ , in a basis αβ in the system’s Li-ouville space (Liouville space is basically a tensor square ofthe Hilbert space), are determined from the matrix elementsof the interaction Hamiltonian: Λ αβα ′ β ′ = 12 n(cid:10) H (cid:11) αα ′ δ β ′ β + (cid:10) H (cid:11) β ′ β ′ δ αα ′ (7) − X j,j ′ ( H int ) j ′ αjα ′ ρ jE ( H int ) jβ ′ j ′ β − (cid:0) hH int i (cid:1) αα ′ δ β ′ β + 2 hH int i αα ′ hH int i β ′ β − (cid:0) hH int i (cid:1) β ′ β δ αα ′ o , where ρ jE are the eigenvalues of the initial environment sta-tistical operator ρ E (0) . Λ contains essential information onthe directions of coherence loss. Strictly speaking, the abovecoarse-graining procedure holds if k Λ k τ ≪ min { , kL eff k τ } . (8)Since the interaction Hamiltonian is linear in the contact cre-ation and destruction operators, and we can approximate thateach contact snaps back to a ”drifted” grand-canonical sta-tistical operator, we have hH L/R int i = 0 . This means that L eff = L S , and also leaves us with only the first three termsin Eq. (7) for Λ to calculate. It can be shown [20] that eachterm in Λ is a sum of independent contributions over individ-ual modes [ Λ = P k Λ k ] that attack only single-particle stateswith a given k . The same holds for L S . Consequently, in re-ality we have a multitude of two-level problems, one for eachstate Ψ k , where the two levels are a particle being in Ψ k (”+”)and a particle being absent from Ψ k (”-”). Each such 2-levelproblem is cast on its own 4-dimensional Liouville space, with ρ k = (cid:0) ρ ++ k , ρ + − k , ρ − + k , ρ −− k (cid:1) T being the reduced statisticaloperator that describes the occupation of Ψ k . According to(6), dρ k dt = [ − i L S,k − Λ k τ ] ρ k , (9) where L S,k = ω k − ω k
00 0 0 0 , Λ k = A k − B k C k C k − A k B k , (10)and A k = ∆ k (1 − f Lk ) , B k = ∆ k f Lk , C k = ( A k + B k ) / k / . (11)The rows/columns are ordered as | + i h + | , | + i h−| , |−i h + | , |−i h−| . Clearly, off-diagonal ele-ments ρ + − k and ρ − + k decay as exp ( ∓ i ω k − τ C k ) t and reachzero in the steady state. The two equations for ρ ++ k = f k ( t ) and ρ −− k = 1 − f k ( t ) are actually one and the same, and eitherone yields df k dt = − τ ( A k + B k ) f k + τ B k = − τ ∆ k f k + τ ∆ k f Lk , (12a)where f k is the distribution function for the active region. Ananalogous relationship holds for the backward-propagatingstates: df − k dt = − τ ∆ − k f − k + τ ∆ − k f R − k ′ (12b)Equations (12) may at first glance appear to be Markovianin form, but they are generally not, as we will discuss in thenext section. However, if we are in the low-bias regime andassume that: (1) the potential and thus the scattering states,transmission coefficients, and the coupling strengths ∆ ± k arevirtually constant throughout the transient, and (2) the currentdensity is low, so any changes to the contact distribution func-tions that result from a current flow can also be neglected, thenevolution (12) will indeed be Markovian [20]. In fact, we cansolve the above equations analytically in the limit of low biasand low current densities. In that case, the contact distributionfunctions are nearly constant, and the solution to Eqs. (16)can be found as f k ( t ) = (cid:0) f k (0) − f Lk (cid:1) e − τ ∆ k t + f Lk , (13) f − k ( t ) = (cid:0) f − k (0) − f R − k ′ (cid:1) e − τ ∆ − k t + f R − k ′ . As expected, the steady-state values of the distribution func-tions are the contact distribution functions f k ( ∞ ) = f Lk , f − k ( ∞ ) = f R − k ′ . (14)A detailed discussion of the relationship of the model with theLandauer-B¨uttiker formalism can be found in [20]. IV. EXAMPLE: TRANSIENT IN AN
NIN
DIODE
As the current starts to flow through the structure, the con-tact distribution functions quickly adjust to accommodate thecurrent flow. A good approximation for the distribution func-tion in bulklike contacts in which electron-electron scatteringis the most efficient mechanism is the drifted Fermi-Dirac dis-tribution function f Lk ( k d ) = 1exp n ~ [( k − k d ) − k F ]2 m || k B T o + 1 , (15)where k d , the drift wave vector, depends on the total currentdensity J flowing through the structure as k d = m || J/e ~ n . m || is the effective mass in the direction of current flow, and n is the contact carrier density. k d changes during the transientand brings about non-Markovian character to Eqs. (16): df k dt = − τ ∆ k f k + τ ∆ k f Lk ( k d ) , (16) df − k dt = − τ ∆ − k f − k + τ ∆ − k f R − k ′ ( k d ) . where it should be understood that k d changes with time.As the transient progresses, the current and the charge den-sity in the structure change, which in turn changes the po-tential profile, the scattering states available to electrons, thetransmission coefficients, and, to a small degree, the interac-tion matrix elements ∆ ± k , as well as the aforementioned con-tact distribution functions. Therefore, all these quantities haveto be carefully updated during the simulation.Equations (16) must in general be solved numerically. Aflowchart of the numerical algorithm is presented in Fig. 2.Upon the application of bias to the contact, but before the cur-rent starts to flow ( t = 0 + ), we solve the Schr¨odinger andPoisson equations with equilibrium initial distribution func-tions f ± k (0) . At this point current continuity between thecontacts and device is not necessary and k d = 0 ( J = 0 ).We then proceed to the next time step, with non-zero current,and solve first the Schr¨odinger equation using the potentialfrom the previous time step. Using the previous value for k d ,we update the distribution functions at the new time step andcalculate the current and charge densities, then find the currentdensity due to the change in the device charge density, iteratefor a new k d until the current density in the contacts is equalto the sum of the current density in the device and the cur-rent density due to the change in the device charge density. Ineach new iteration, we use a k d that is formed as a weightedsum of k d ’s from the current and previous iterations. Whenthe current and k d are self-consistently obtained, we use thenewly obtained device charge density in the Poisson equationto obtain a new guess for the potential, and repeat until thepotential converges (using the globally convergent Newton’smethod [27] with a semiclassical Jacobian [28, 29]). We re-peat for all time steps until a steady state is reached.There are several nontrivial numerical considerations. Oneis the ability to achieve a high enough density of scatteringstates to properly represent physical quantities such as thecharge and current densities or the potential profile. We aretrying to capture a continuum of scattering states, which atfirst glance might seem doable by indiscriminately increasingthe density of k ’s by choosing larger and larger simulationdomains. Unfortunately, this brute-numerical-force approachdoes not work; what does work instead is generating a ”smart”discrete set of scattering states by first solving the Schr¨odinger Thomas-Fermi initial guessSchrödingerequationUpdate distributionfunctions at time i*dtCalculate Kd fromcurrent continuityVoltage stepappliedPoissonequation Potentialconverged?Kd converged?Steadystate?I(t), V(t), n(t)No Yes NoYesYesNo i=i+1Self-consistent solution withequilibrium distribution functionsi=1
FIG. 2. Flowchart of the numerical algorithm for the calculation ofthe electronic response of a biased two-terminal nanostructure duringa transient. equation in the simulation domain with the condition that thefirst derivative be zero at the boundaries, and then projectingthese states onto the forward and backward moving scatter-ing states. Details of this discretization of the scattering statecontinuum can be found in [29].A related question is how to populate a small set of boundstates that can emerge in a biased nanostructure (e.g. notethe potential pocket on the left-hand-side of Fig. 3a at timesbelow 100 ps). Those states are in reality filled by electron-electron and electron-phonon scattering, essentially the samemechanisms as in the contacts. Since we are not treating scat-tering explicitly in this approach, we populate the bound statesaccording to the Fermi level in the nearest contact. More detailon the finer points of the numerical simulation can be foundin [21].Figures 3 and 4 depict the potential, charge density, andcurrent density for a single ellipsoidal valley in an nin silicondiode at room temperature. The left and right contacts aredoped to cm − , whereas the middle region is intrinsic(undoped). The momentum relaxation time in the contacts istaken to be τ =
120 fs, based on the textbook mobility val- x [nm] po t en t i a l [ m e V ] t = 0 pst = 3 pst = 10 pst = 100 pst = 1000 ps0 50 100 150 200 25000.511.52 x 10 x [nm] e l e c t r on den s i t y [ m − ] t = 0 pst = 3 pst = 10 pst = 100 pst = 1000 ps FIG. 3. (a) Potential and (b) charge density in the nin diode as afunction of time upon the application of -25 mV to the left contact.The n -type regions are doped to cm − . The contact momentumrelaxation time is τ =120 fs, as calculated from the textbook mobilityvalue corresponding to the contact doping density. ues for the above doping density. Note that the characteristicresponse time of the current is of order hundreds of picosec-onds, so three orders of magnitude greater than the contactrelaxation time. The transient duration is long because of therelatively weak coupling between the active region and thecontacts; the transient duration can be thought of as the in-verse of a typical ∆ k τ among the k ’s participating in the cur-rent flow. V. CONCLUSION
We presented a theoretical treatment of the transient-regimeevolution of an electronic system in a two-terminal ballisticnanostructure coupled to dissipative contacts. The approachis rooted in the open system theory and is based on two keyingredients: (1) A model interaction Hamiltonian between the active region and the contacts, constructed specifically to con-serve current during the process of carrier injection from/intothe contacts, whose matrix elements are readily calculatedfrom the single-particle transmission problem for structureswith and without resonances alike. (2) In the absence of scat-tering in the active region, it is the rapid energy relaxationin the contacts (due to electron-phonon or, in good, highly-doped contacts, due to electron-electron scattering) that is the t [ps] J [ A / m ] FIG. 4. Current density versus time for the nin diode from Fig. 3upon the application of -25 mV to the left contact. The n -type regionsare doped to cm − and τ =120 fs. indirect source of irreversibility in the evolution of the current-limiting active region, owing to the contact-active region cou-pling. We account for the influence of the rapid relaxation inthe contacts by coarse graining the exact active region evo-lution over the contact momentum relaxation time. The re-sulting equations of motion for the distribution functions ofthe forward and backward propagating states in the active re-gion, Eqs. (16), have non-Markovian character as they incor-porate the time-varying contact distribution functions throughthe time-dependent drift-wavevector that depends on the in-stantaneous current flowing. In order to obtain the responsequantities of a nanostructure under bias, such as the potentialand the charge and current densities, the non-Markovian mas-ter equations must be solved numerically together with theSchr¨odinger, Poisson, and continuity equations. We presentedan algorithm for the numerical solution of this coupled systemof equations and illustrated the approach on the example of asilicon nin diode. VI. ACKNOWLEDGMENT
This work has been supported by the NSF, award ECCS-0547415. [1] M. Lundstrom, Fundamentals of Carrier Transport (CambridgeUniversity Press, Cambridge, 2000).[2] D. K. Ferry and S. M. Goodnick, Transport in Nanostructures(Cambridge University Press, Cambridge, UK, 1997). [3] M. V. Fischetti, J. Appl. Phys. , 270–291 (1998).[4] M. V. Fischetti, Phys. Rev. B , 4901–4917 (1999).[5] D. K. Ferry, R. Akis, J. P. Bird, M. Elhassan, I. Knezevic,C. Prasad, and A. Shailos, J. Vac. Sci. Technol. B , 1891– , S254–S256 (2004).[7] R. Gebauer and R. Car, Phys. Rev. Lett. , 160404 (2004).[8] R. Alicki and K. Lendi, Quantum Dynamical Semigroups andApplications, Lecture Notes in Physics, Vol. 286 (Springer-Verlag, Berlin, 1987).[9] H. P. Breuer and F. Petruccione, The Theory of Open QuantumSystems (Oxford University Press, Oxford, 2002).[10] Y. Meir and N. S. Wingreen, Phys. Rev. Lett. , 2512 (1992).[11] A. P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B , 5528(1994).[12] P. Lugli and D. K. Ferry, IEEE Trans. Electron Devices ,2431–2437 (1985).[13] M. A. Osman and D. K. Ferry, Phys. Rev. B , 6018 (1987).[14] A. M. Kriman, M. J. Kann, D. K. Ferry, and R. Joshi, Phys. Rev.Lett. , 1619–1622 (1990).[15] W. R. Frensley, Rev. Mod. Phys. , 745 (1990).[16] W. P¨otz, J. Appl. Phys. , 2458 (1989).[17] R. Brunetti, C. Jacoboni, and F. Rossi, Phys. Rev. B , 10781(1989).[18] M. Nedjalkov, H. Kosina, S. Selberherr, C. Ringhofer, and D. K. Ferry, Phys. Rev. B , 115319 (2004).[19] S. Datta and M. P. Anantram, Phys. Rev. B , 13761 (1992).[20] I. Knezevic, Phys. Rev. B , 125301 (2008).[21] B. Novakovic and I. Knezevic (2012), in preparation.[22] X. Q. Li, J. Y. Luo, Y. G. Yang, P. Cui, and Y. J. Yan, Phys. Rev.B , 205304 (2005).[23] J. N. Pedersen and A. Wacker, Phys. Rev. B , 195330 (2005).[24] B. L. Altshuler and A. G. Aronov, JETP Lett. , 514 (1979).[25] B. L. Altshuler and A. G. Aronov, Solid State Commun. , 11(1981).[26] D. A. Lidar, Z. Bihary, and K. B. Whaley, Chem. Phys. , 35(2001).[27] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flan-nery, Fortran Numerical Recipes, Vol. 1, Numerical Recipes inFortran 77: The Art of Scientific Computing, 2 edition (Cam-bridge University Press, New York, 1992).[28] R. Lake, G. Klimeck, R. C. Bowen, and D. Jovanovic, J. Appl.Phys. , 7845 (1997).[29] S. E. Laux, A. Kumar, and M. V. Fischetti, J. Appl. Phys.95