Theory of Electron Transport in Small Semiconductor Devices Using the Pauli Master Equation

It is argued that the Pauli master equation can be used to simulate electron transport in very small electronic devices under steady-state conditions. Written in a basis of suitable wavefunctions and with the appropriate open boundary conditions, this equation removes some of the approximations which render the Boltzmann equation unsatisfactory at small length-scales. The main problems consist in describing the interaction of the system with the reservoirs and in assessing the range of validity of the equation: Only devices smaller than the size of the electron wavepackets injected from the contacts can be handled. Two one-dimensional examples solved by a simple Monte Carlo technique are presented.


INTRODUCTION
As noted by Frensley discussing quantum trans- port in open systems [1], the Pauli master equation [2], (PME) could constitute an intuitive description of electron transport in semiconductor devices of size comparable to the electron wavelength: One could capture the wavelike nature of transport, lost in the Boltzmann-transport-equation (BTE) picture, so bypassing the weak-field, long devices limitations.In retaining the weak-scattering and completed-collision limits one would not reach a full-fledged quantum description, as required to study very fast time transients [3,4].Yet, as devices approach the sub-50 nm length-scale, the advantage of avoiding the concept of point-like electrons appears overwhelming.
Although arguments have been raised against the correctness of the PME [1], it is applicable to very small devices and to steady-state phenomena: Small devices, because the PME rests on the absence of off-diagonal elements of the density matrix p.In turn, their absence implies that the contacts inject 'plane waves', which corresponds physically to the injection of spatially highly delocalized wavepackets.Hence, the device must be smaller than the 'dephasing length' in the contacts [5].For heavily doped Si, this is of the order of 50 nm, which is of technological interest.Steady state, since any nontrivial time transient (e.g., relaxation of photo-excited carriers, fast switching of a device) would excite those off- diagonal terms of p which the PME cannot handle.Indeed, as also note by Frensley [1], the PME is inconsistent with current continuity in this case.

MASTER EQUATION AND CONTACTS
Let's consider a semiconductor region unbounded in the (x, y)-plane, homogeneous in this plane, in contact with two reservoirs at z 0 and z L. Let's write the single-particle Hamiltonian as: e V(z) + Hint -[-Hres H0 t_ Hint -[-gres, within the envelope/effective-mass approximation, m* being the effective mass.Hint describes the electron-phonon interaction or interactions with impurities, etc. nre represents the interaction between the device and the external reservoirs.
The time evolution of the system is given by the 'transport' equation for the 'reduced' (electronic) The last term on the right-hand-side describes the effect of the reservoirs.Assuming Hint a H', we can derive a transport equation from Eq. ( 2), by considering only terms to the leading order in a, and following the 'standard' derivation of Kohn and Luttinger [7], but now using as basis functions not plane waves, but the eigenstates IK) of H0, such that HoI#K)= [E, + h2K2/(2m*)]l#K)= E,(K)I#K).This simplifyies Eq. ( 2), since the commutator now contains only the interaction term Hint.The wave-functions are simply of the form u(z) exp (i K. R)/(27r), where K and R are the wavevector and position on the (x, y)-plane, and the functions (u(z) obey the one-dimensional Schr/Sdinger equation for the potential V(z), subject to the boundary conditions V(0)= VL and V (L) VR.Taking -e Vr > -e Vn, for any given energy Eu>-eVr, two independent solutions represent waves incident from the left ( + ) or right (-), and partially reflected and transmitted, with reflection and transmission amplitudes r + and +.
For -e V > E u > -e V, the nondegenerate solu- tions represent non-current-carrying states.In the following the index a 4-will label leftand right-traveling waves in the energy range Eu> e VL, and (zl#a) (z).
To the leading order in a, the diagonal elements, pux, (symbol used in place of puxux,) of the density matrix obey the Eqs. 7-9]: uK'a' ##Ka ( Op,x + \ Ot Je" here Wx,,;ux represents the probability per unit time of a transition from the state [#Ka) to the state [u K'a ) using Fermi golden rule.While the derivations provided in [7-9] can be followed very easily, now one must also assume that the reservoirs do not alter the ratio a between off- diagonal and diagonal matrix elements.Van Hove [8] has shown the validity of Eq. (3) in describing approach to equilibrium, in the case of closed systems (for which the term (Op/Ot)res is absent), provided the initial state is chosen as 'quasi- diagonal': At t=0 we must have puxvX,o=O for Eu(K) Ev(K') > 6, where 6/h is of the order of the relaxation rate, 1/r.Similarly, Eq. ( 3) is valid as long as the contacts do not 'inject' off-diagonal elements, that is (Opuxx,,/Ot)res 0 for Eu(K) E(K) > 6.This is the case for very delocalized wavepackets entering the contact, plane waves being the ideal limiting case.In this limit, the effect of the contacts can be modeled phenomenologically by assuming that the left reservoir attempts to restore charge neutrality by injecting a flux of particles of any K into the right-bound component of the states associated to waves incident from the left, +, fixed by the Fermi level at the cathode, E.
Accounting also for the flux out of the device via either contact, one can write" Ot res where [2m*(E.+eVR/L)] 1/2, and the A's are normali- zation constants.Similar expressions also hold for states traveling from the right.The occupation of the states at thermal equilibrium, as determined by the 'left' reservoir, has been expressed in terms of the equilibrium Fermi function (integrated over the 'transverse' wavevectors K), f(th,L) (Eu) pos- sibly shifted in k-space, in order to account for the open nature of the contacts [1,6].By analogy, a similar expression (with the terms p,+-Tp,_ replaced by Pu) is also used to express the linear relaxation of the of states with energy -e Vn < Eu < -e VL.Finally, the net current flowing through the device is given by: //kS + + 2 J -e p.+,m--7 It.A,   where obviously only the doubly-degenerate states with E, >-eVr contribute to the sum.
Criticisms have been raised against Eq.(3) based on the violation of continuity [1].Yet, it should be noted that at steady state the reservoirs and the reverse scattering-induced transitions provide the 'missing' current and current continuity is trivially satisfied.However, continuity is indeed violated when using the PME to describe 'fast' time-dep- endent phenomena: For example, any time depen- dent perturbation V (t) will cause an electron in I#Kr) at 0 to have nonzero amplitude also in a state IK'r').At long times the requirement that this amplitude be negligible translates into the condition leO6V/Otl << h/7", expressing the fact that the applied bias cannot change appreciably over a relaxation time, so that the PME is unable to handle nontrivial time-dependent phenomena.

MONTE CARLO SOLUTION OF ONE DIMENSIONAL EXAMPLES
In order to implement Eq. (3) in one-dimensional situations, starting from a given potential V(z), such as a Thomas-Fermi solution, the Schr6dinger equation is solved with vanishing Neumann boundary conditions.The eigenvalues {E} will be spaced densely enough (to accurately reproduce the bulk density of states in the reservoirs) if the length L of the device is selected appropriately.The extended states can be finally decomposed into left-and right-traveling waves following the 'quantum transmitting boundary method' by Lent and Kirkner [10].The levels I#) are populated following P6tz [6] and the Poisson equation can be solved again, using well-known potential-damping schemes to improve convergence.Iterating this Poisson/Schr6dinger scheme yields a solution in the ballistic case.Scattering is introduced by calculating all transition rates Wx,,;ux between the levels, and adding the PME as a third equation over which one must iterate.A Monte Carlo algorithm is used here.Scattering and injection/ exit from contacts are treated stochastically in order to determine the new populations Pux.
Typically, 300,000 'electrons' are employed, using a time step t 10 -16 S in the Monte Carlo algorithm.After about 100-500 Monte Carlo steps, the Poisson and the Schr6dinger equations are solved again until convergence is obtained.Finally, the current flowing through the device is obtained from Eq. ( 5), while information about carrier density, kinetic energy, and velocity as a function of position inside the device can be obtained from the expectation values of the corresponding operators as traces over their products with the density matrix.
As a first example, let's consider a simple nin resistor at 300K consisting of a 300 nm-long region between z 0 and z L with n-type regions 100 nm-long n-type doped to 1017 cm-3, and an intrinsic region also 100 nm long.A uniform effective mass of m* 0.32 m0 (where m0 is the free electron mass) and a dielectric constant 11.7 eva (where Cvao is the permittivity of vacuum) are used.Nonpolar scattering with optical phonons is included using a photon energy 60 meV, an optical deformation potential of 5 108 eV/cm, a crystal density Px 2.33 g/cm3.These parameters imply a 'dephasing length' , in reservoirs of the same material of the order of 785 nm, larger than the active region of the device, about 100 nm long.
Figure shows the potential, carrier density, kinetic energy, and velocity at the end of the iteration.The dashed lines are the corresponding Potential and charge density (top frame), kinetic energy and drift-velocity (bottom frame) as a function of position z inside an nin device calculated using the master (solid lines) and Boltzmann (dashed lines) equation.The charge accumulation at the device boundaries is an artifact caused by boundary conditions employed to solve the Schr6dinger equation.
quantities obtained from a conventional Monte Carlo solution of the BTE using identical para- meters.The current density is about 6.78 x 104 A/ cm2, about 10% smaller than the value obtained using the BTE, 7.58 x 104 A/cm 2. This difference originates mainly from a slightly higher built-in barrier at the cathode/intrinsic-region junction seen in the 'quantum' results, caused by the penetration of wavefunctions into the barrier itself.The penetration of charge into classically forbidden region is also the origin of the redis- tribution of charge seen in the top frame of Figure 1" The charge density obtained using the PME is larger in the intrinsic region, as wavefunc- tions penetrate into the gap beyond the classical turning point.This charge redistribution is the major factor responsible for the differences seen in average kinetic energy and drift velocity.
The second example is a simple GaAs/A10.3_Ga0.TAs resonant tunneling structure.Two 150 nm thick layers of GaAs-like semiconductor with n- type doping (2 x 1017 cm-3) are separated by two A1GaAs tunnel barriers, each 2.8 nm thick, enclosing a 4.5 nm wide semiconductor well.Two 3 nm-thick undoped 'spacer' layers, in turn, enclose the double-barrier region.'Conventional'   parameters for GaAs and A10.3Ga0.7Asare used [12], also including nonpolar scattering with acoustic phonons and polar scattering with long- itudinal optical phonons.A lattice temperature of 300 K is used.The contact mobility used to displace the distributions injected by the reservoirs is taken to be el 5000 cm2/Vs, corresponding to a dephasing length in the reservoirs, 130 nm, larger than the active region of the device.Since even at room temperature scattering is not sufficiently strong to fully populate the states in the accumulation layer, as discussed by Frensley [1], the curves in Figure 2 show that a large fraction of the applied bias drops over the cathode region to the left of the double-barrier region.Although similar results have been obtained by Frensley [1] and Kluksdahl et al. [11], they appear poorly credible, since the potential profile clearly depends on the size of the simulated region.The curves are parametrized by the applied bias from 0.0 to 0.6 V in steps of 0.05 V.Note the 'gaps' at the left of the barriers around 0.2 and 0.45 V. Convergence could not be obtained close to these biases, because of an instability caused by scattering, as discussed in the text.The dotted solutions for biases of 0.2 and 0.45 V are only representative of the average potential during the iteration.
Figure 3 shows the current-voltage character- istics: As the resonance is approached (at about 0.2 V), the associated build-up of charge in the well screens the cathode field, increasing the field scattering.Note the 'snap back' of the current at biases of 0.2 and 0.45 V, corresponding to the instability regions of Figure 2.
The error bars are used to illustrate the range of the oscillation of the current.The inset shows the Calculated current density vs. simulation time at two bias points, slightly away from (0.175 V) and close to (0.200 V) resonance.Note the instability of the current at the bias point near resonance.
reached at which the rate at which scattering feeds charge into the states localized in the well becomes smaller than the rate at which the well-charge leaks to the right.Charge is removed from the well, and the potential profile now 'snaps', a larger voltage now dropping at the left of the barriers.The population of the current-carrying scattering states, and so the current itself, now drops.As the applied bias is increased further, another 'snapping point' is reached, at about 0.45 V.These critical bias points can be seen as 'gaps' in the cathode potential-profile in Figure 2.
As loosely indicated by the error bars in Fig- ure 3, for biases close to and slightly beyond resonance the situation is unstable: As 'resonant wavefunctions' leak out to the anode, scattering to non-current-carrying states in the collector region is enhanced, thanks to a larger overlap (form) factor.The resulting redistribution of charge alters the potential profile, driving the system away from resonance.As this happens, the charge in the anode side of the barriers drops once more, and the system is now driven back towards a resonant condition, and the cycle repeats once more.The inset of Figure 3 illustrates some features of this instability.The good convergence properties of the computed current density for a bias slightly away from resonance (0.175 V) are compared to the unstable behavior for a bias close to the first resonance (0.200 V).In the latter case, with some imagination one may detect an oscillatory behav- ior between two states [11,13].

SUMMARY AND CONCLUSIONS
The main conclusion of this paper is that there is an interesting class of problems for which the Pauli master provides a solution of the electron trans- port problem overcoming some of the strongest limitations of the semiclassical Boltzmann equation: For devices with active regions smaller than a few hundreds of nm, the master equation makes it possible to account for the wave nature of the electrons, although only in the weak scattering limit and in steady state or slow (adiabatic) time transients.The major price one pays is the inability to handle fast time-transients and strong scattering situations, collisional broadening and coherent transport in the femtosecond time scale being the most notable examples.The advantages consist in the ability to treat inelastic scattering processes, arbitrarily strong electric fields (including intra- collisional field effects) tunneling phenomena, elastic and not, and, in general, all those effects which depend on the wavelike nature of the electrons.

FIGURE 2
FIGURE 2 Potential energy vs. position in a resonant tunnel diode calculated including electron-phonon scattering.The curves are parametrized by the applied bias from 0.0 to 0.6 V in steps of 0.05 V.Note the 'gaps' at the left of the barriers around 0.2 and 0.45 V. Convergence could not be obtained close to these biases, because of an instability caused by scattering, as discussed in the text.The dotted solutions for biases of 0.2 and 0.45 V are only representative of the average

FIGURE 3
FIGURE 3 Current density in the GaAs/A1GaAs resonant tunneling structure with the inclusion of electron-phonon