The Methodology of Simulating Particle Trajectories Through Tunneling Structures Using a Wigner Distribution Approach

Understanding and describing steady-state and transient tunneling behavior are necessary for determining the fundamental transport properties and speed performance limits of electronic devices based on field emitter arrays. The Wigner distribution approach is used as a first step in developing a theoretical and calculational methodology applicable to a wide class of quantum tunneling problems. The tunneling behavior in nanostructure field emitters is fundamentally similar to the tunneling behavior in other tunneling structures, e.g., resonant tunneling diodes (RTD’s), where a large amount of experimental data exist. A modified relaxation time approximation is used for the scattering. Self-consistency is implemented, for the first time in fully time-dependent calculations, to investigate both the steady-state and transient behavior in RTD’s. Of particular interest is their ultra-fast switching times (approximately 0.1 ps), which is comparable to that promised in micrometer scale vacuum devices based on field emitter arrays. A methodology for constructing quantum “particle trajectories” is presented with the long-term goal of incorporating them into a particle Monte Carlo simulation. The trajectories for an RTD are presented and their behavior, as a function of scattering and selfconsistency, is shown to be consistent with the steady-state current-voltage/quantum well electron density characteristics of the RTD, and with the response of the RTD to a sudden bias switch. The trajectories also exhibit a conservation-of-energylike behavior. The trajectory formulation is thus shown to be potentially useful for incorporating into a multidimensional particle Monte Carlo simulation of quantum-based devices in which the tunneling region is small compared to the dimensions of the device, as is the case with heterojunctions and nanostructure field emitters.


I. INTRODUCTION PROBLEM of central concern to modeling in the field
A of vacuum microelectronics is the ability to characterize the space-and time-dependent tunneling process of electrons on a field emitter.Conventional treatments of tunneling, such as calculation of transmission coefficients and currents [ 11, fail to account for the space-and timedependent elementary events in a quantum tunneling process which will ultimately limit the speed performance of vacuum microelectronic devices.In most problems of practical and important technological applications, tunneling may occur between any two different media, and general techniques for dealing with this problem obviously have broad applications.In principle, then, if a means for simulating tunneling in a metal-vacuum interface may be found, then it could be used to treat tunneling in, for example, a semiconductor heterostructure such as the Resonant Tunneling Diode (RTD), or vice versa.Indeed, the two cases are not that different from a simulation point of view: As we shall see, the two cases differ only in the use of the boundary conditions, but not the method of solution for the Wigner function.From the point of view of particle trajectories, therefore, the differences between a metal-vacuum interface and an RTD are immaterial, as a phase space description is applicable to both.In this paper, we introduce trajectory methodology to describe elementary space-and time-dependent events in a tunneling process in the RTD, due to the large body of experimental data that exists which allows us to validate our model.Work is being done on extending this approach to the metal-vacuum interface.
Investigating transport in solids typically uses the Boltzmann equation, which allows for the inclusion of scattering phenomena [2], [3] and has a natural particle interpretation.By assuming slow spatial and temporal variations, the quantum Boltzmann equation may be formulated, but the slow variation assumptions preclude full quantum tunneling [4].The Wigner distribution function [ 5 ] , which acts like a quantum probability distribution function, has emerged as one of the most promising candidates for modeling quantum transport [6]-[ 151, and the inherent phase-space description in the Wigner distribution function (WDF) approach allows for a natural particle interpretation and takes into account scattering and self-consistency [ 101, [ 1 11, [ 141.We have emphasized the particle trajectory representation of the quantum transport processes for use in developing a full-fledged ensemble particle Monte Carlo device analysis because of the following: One of the most important and useful features of the Wigner function is that it can rigorously model quantum transient phenomena [6], [12], such as when the bias across an RTD is suddenly changed.However, present computer resources to simulate these phenomena limit the size of the simulation U.S. Government work not protected by U.S. copyright region.Further, assumptions must be made about the velocity distribution of the electrons entering the simulation region [ 1 3 ] .Ensemble particle Monte Carlo (EPMC) is recognized as a powerful numerical technique for including all charge-carrier scatterings, many-body effects, selfconsistency, and multi-dimensionality, and it is able to model devices which are very "large" by comparison to the simulation regions used in quantum transport.However, EPMC cannot model space-and time-dependent quantum tunneling, and the "particle trajectory" interpretation we have been advancing is designed to bridge the gap between the two methods by extending a classical trajectory concept (used by Monte Carlo) into the "quantum trajectory regime."

THE WIGNER DISTRIBUTION FUNCTION: CONTINUUM EQUATIONS
The Wigner function f (x, k; t ) may be defined as the Weyl transform of the one-particle density matrix operator p(x, x'), representing either pure, mixed states, or a pair-correlation function of a complicated many-body scattering system [ 151 For present purposes, it may be thought of as a "quantum probability distribution function" in momentum and position, though one should note that it can take on negative values, which distinguishes it from a true probability distribution.Its time evolution is given by where L is the "Liouville operator," C is the collision (scattering) operator [ 1 4 ] , and V(x) is the energy bandedge profile of the RTD.The presence and form of the scattering operator may be justified rigorously in the context of a phase-space nonequilibrium Green's function approach to quantum transport [ 1 5 ] , and it is analogous to the equivalent term in Boltzmann's equation in the zerothorder gradient expansion.The time evolution equation may be recast in the following term, which is convenient for a numerical evaluation: ( 3 ) The particle trajectories [9] may be introduced by drawing an analogy to the Liouville's theorem on the conservation of particle density in phase space, namely, s, Comparing (4) with ( 2 ) , the velocity of the fictitious particles are given by ( h k / m * ) .The "effective quantum force" is given by By virtue of V(x k y ) , the quantum force is manifestly nonlocal and therefore accounts for tunneling and also affects the behavior of electrons far from the barrier region.In effect, we can treat the particles as though they have well-defined momenta and positions, i.e., we treat them "classically." The quantum-mechanical tunneling enters through the modified "quantum force" these particles experience.Such a treatment is close in spirit to the Bohm formulation of quantum mechanics in terms of "quantum potentials" and trajectories If scattering is to be included inside the simulation region, then a scattering operator must be appended to the time evolution equation.The general collision operator, which is used to model scattering within the simulation region, is given by ~2 4 1 .
( 7 ) The right-hand side (RHS) may be approximated by a modified relaxation time approach, where by "modified" we mean that the collision operator is not simply given by -(ffo)/7; the operator used must satisfy the "condition of detailed balance" [ l o ] .If we define the density of particles as p(x; t ) = (2?r)-' j F f ( x , k; t ) dk, and require a constant relaxation time 7, then the form of W is uniquely determined, and the modified relaxation time approximation to the collision operator becomes We must then determine what the equilibrium distribution fo is.A natural choice is to let it be given by the unbiased steady-state distribution, i.e., L fo = 0, so as to account for all space charges and the effect of particle depletion near the barrier regions.
We have found 7 by numerically integrating over the momentum relaxation times 7i(x), where x = P A k 2 / 2 m * , " i " refers to the various processes (acoustic, zero-order optical, polar optical phonons, and piezoelectric and ionized impurity scattering), and 0 = l / k T (9) where po(x) is the doping density, c0 is the permittivity of free space, and K, is the static dielectric constant.The band-edge discontinuity A E&) interface conditions are then added to +(x).In the self-consistent approach, fix, k ) must be known to evaluate p(x), which gives the self-consistent potential needed to evaluate f (x, k ) .We have chosen to calculate the initial self-consistent potential using our transient simulation program starting with a neutral charge distribution and then allowing the system to time-evolve self-consistently to steady state (this approach is substantially different from that of either [7] or [ 1 3 ] ) .The program for the Wigner function then has the virtue of being completely self-contained and applicable to any arbitrary potential profile.The virtue of this approach is that it dynamically finds the self-consistent steady-state solution, if one exists.As it turns out, the current for constant bias in the negative differential resistance region (NDR) of the RTD oscillates at a very high frequency, and therefore only averaged quantities are uniquely defined.These oscillations have been simulated for the first time using the above approach, and address a long-standing controversy about the nature of bistability and current oscillations in RTD's (see [25] and references therein).The detailed discussion of this important result will not be treated here.
Finally, we must evaluate the trajectories.Using (4), we see that the Wigner distribution-function equation can be identified with the Liouville's theorem on the conservation of particle density in phase-space.Thusf(x + 6 x , k + 6k; t + St) = f ( x , k; t ) .The particle trajectories consequently lie on contour lines of the Wigner function.However, merely plotting contour lines is not enough.A more physically meaningful approach to getting trajectories is to estimate the next point in phase space by solving Newton's equations, and then making a correction by pushing the point onto a contour line.Consequently, our program for evaluating trajectories is to make an initial estimate for the position x' and momentum k' of the particle using ( 1 1 ) and then moving along the line of steepest descent to the contour.This approach is not as susceptible to error as a pure dynamics approach is because the contour line is guaranteed to correspond to the trajectory, and therefore, longer time steps may be taken; it is also better than a pure contour approach which has ambiguities due to saddle points, surface convolutions, and numerical "noise."

DISCRETIZATION OF THE CONTINUUM EQUATIONS
There are four tasks before us: 1 ) we must discretize and solve (2); 2 ) we must evaluate the relaxation time; 3) we must develop a procedure to incorporate self-consistency in a transient simulation, and 4) we must find a numerical method for finding the contour lines of a timeevolving Wigner function.

A . Discretization
We now turn to the discretization of the "Liouville operator."We are interested in simulating a region between x = 0 and x = L for both positive and negative momenta, since particles will flow in both directions at the boundary.If we let Nn and N be the number of x and k points on a grid in phase space, then we will evaluate the Wigner function on the points x ( i ) and k ( j ) such that The value of 6k may be determined by noting that V(x, k k ' ) represents a Fourier transform; to insure that the discrete Fourier transform is invertible, we thus require that 6k = a / N S x .Consequently, the discretized Liouville operator becomes where A* is the finite difference approximation to the derivative.A simple approximation to this is the two point (UDS, or Euler) scheme [6], given by ( where UDS stands for "up/downwind differencing scheme."The scheme is accurate to order O(E).This scheme is advantageous in terms of programming ease and computer memory, but it has difficulty modeling the RTD near resonant voltages (i.e., the voltage at which the current peaks in the negative differential resistance (NDR) region).We have thus introduced a three-point secondorder differencing scheme (SDS) accurate to order O(E*) 1151: if a functionf(x) is given on three (equally spaced) points, then a unique quadratic polynomial may be fitted through these points.f (x) is then approximated by where x = x ( i ) + p 6 x and -1 I p I 1. Taking the derivative with respect to x yields Letting p = 0 gives the familiar central difference schemes (CDS).This scheme, as we shall see below, is only valid for simulating unbiased RTD's [9].To simulate a biased RTD, we need to mimic the up/downwind behavior of the Euler approach; thus the Second-order Differencing Scheme (SDS) is given by letting p = rfr: The choice of ( f ) depends intimately on the boundary conditions.In our simulations, we make the assumption that the distribution of particles entering the device is known a priori; particles enter at the boundary given by (x = 0, k > 0) and (x = L, k < 0).For the other two boundaries (x = L, k < 0) and (x = L, k > 0), the distribution of exiting particles is determined by the behavior of the device.These "ohmic" boundary conditions are well-known.The idea is widely used in plasma physics and Monte Carlo device simulations, and is discussed at length in [13], and also in [7].For i = 1, A-requires values off@, k) for x < 0, and for i = Nx, A + requires values forx > L. Consequently, we shall use A-when k < 0 and A+ for k > 0 (Fig. 1).The distribution of incoming particles shall be assumed to be the three-dimensional Fermi-Dirac distribution with the transverse momentum components integrated out.

B. Scattering discretized density of electrons is given by
The scattering operator may be found by noting that the Consequently The calculation of 7 for GaAs parameters has been given in detail before [lo]; for present purposes, we note that for a temperature of 77 and 300 K, r = 525.2and 97.4 fs, respectively.
The calculation of the current across the device is a bit more complicated than the simple formula for density.At steady state, the current across the device is constant, i.e., aJ/ax = 0.By invoking the continuity equation, we then find that in steady state a a ( 2 0 ) Evaluating 13,f in discrete form follows the rules previously set up for A'.Further, we approximate a,J by [J(i + 1 / 2 ) -J(i -1 /2)] / 2 6 x .Discretizing the right-hand side of (18), regrouping terms, and demanding that a,J vanish in steady state, we find [15] hsk 21) wherej' = j + N/2, i.e., k ( j ) corresponds to k < 0 and k ( j ' ) corresponds to k > 0.

C. Selfconsistency of the following matrix equation:
The inclusion of self-consistency requires the solution where +(O) = 0 and +(Nx + 1) = bias difference across the device (compare to (10)).A fast and stable algorithm which makes full use of the simplicity of Poisson's equation in one dimension is given below: let the right-hand side of (22) be given by z(i).On input, then, z ( i ) is proportional to the density; on output, z ( i ) = + ( i ) : The potential V(x) is then the sum of +(x) and the band edge discontinuity AE,(x).
To evaluate the time evolution of the Wigner distribution function, we use the Cayley form of an exponential operator, where, for an operator As the matrix corresponding to the "Liouville operator" is very large, we rewrite (3) so as to minimize matrix operations

D. Trajectories
The method of evaluating trajectories is simply expressed in (4).However, the method for evaluating the derivatives off@, k) given the pointsf(i, j ) needs elaboration.In one dimension, there is a unique nth-order polynomial that passes through IZ fixed points; an analogous result does not obtain in more than one dimension [26].The method we have used in the current work requires that f ( x , k) be relatively smooth so that the interpolating polynomials do not oscillate dramatically.This requires that a finer grid mesh be used if f ( x , k) becomes more oscillatory, for example, near resonance.
The latter is used in the evaluation of the "quantum force," both are used in calculating the correction to (x', k ' ) .If we let x ( t + St) = x ' + e,, and k(t + St) = k' + ek, where e, and ek represent the small terms which will correct the estimates x' and k', respectively.Then the e's may be found from and the fact that the line along which the correction is made is perpendicular to the contour lines (direction of steepest descent).

I v . APPLICATION TO DOUBLE-BARRIER STRUCTURES
In this section, we shall apply the formalism of the previous sections to a problem of current interest, namely, the Resonant Tunneling Diode (RTD).If the drain end of the RTD were replaced by a vacuum, then the problem would become analogous to problems currently of interest in vacuum microelectronics; indeed, the formalism is ideally suited to deal with any one-dimensional potential structure with differing "materials" at the source or drain end.We shall investigate in what ways the characteristics of the device may be understood by the behavior of the trajectories, and we shall examine how these characteristics are affected by scattering and self-consistency .

A. Properties of the RTD
Our model of the RTD will be a one-dimensional qua!tum well structure in which a layer of GaAs around 50 fi is sandwiched between layers of AlGaAs around 30 A, which in turn is sandwiched by GaAs again.A simple understanding of the steady-state behavior of the quantum well under bias may be obtained by investigating the transmission coefficient [ 13, [ 191.Note that the treatment below is meant to be illustrative (more accurate formulations exist [20]-[22], and predict resonance currents greater by roughly a factor of 2 , making them comparable to the currents calculated in the Wigner approach).In general, however, we have found that calculations of currents, densities, and self-consistent potentials from the wave-function approach are time-consuming compared to the Wigner function approach.Moreover, using the wave function approach is inefficient if not impractical by comparison to the Wigner function approach for time-dependent problems.
For a given potential profile V(x), Schrodinger's equation is For incident plane waves from the left given by &(c) = exp (ikx), far to the right where

m * E / h 2 ] , and k' = J[2m*(E -V h ) / h 2 ] . The quantity T(k) = I t(k)(' is defined as the transmission coefficient [ 2 7 ] .
For incident kinetic energies below the barrier height, solutions to Schrodinger's equation are exponentially decaying, implying that T(k) is typically small.However, at particular momenta (depending on the bias), resonance occurs, and T(k) becomes larger; for the unbiased case, T(k) can approach unity even though the incident kinetic energy is below the barr!er height (for our parameters, this occurs at k = 0.0394 A-').As a result, an incoming distribution of particles can still generate a current flow through the device even though the number of particles with a kinetic energy greater than the barrier height is insignificant.The difference between net particles flowing to the left and flowing to the right gives rise to an overall current across the device, which may be crudely calculated by (where Ek = h 2 k 2 / 2 m * and where f,(E(k)) is the incoming distribution of particles)

(hk/m*) dk [fb(Ek) -&(Ek + eV,,)l T(k). (28)
Analyzing the overlap between the transmission coeficient and the incoming particles as a function of bias clearly indicates that one may expect the current to initially increase as a function of bias, encounter a region where the current decreases as the bias increases (Negative Differential Resistance, NDR), and finally, resume increasing as the bias increases (Figs. 2 and 3 ) .Self-consistency has been included in the transmission coefficient calculations; when it is, the qualitative effects are quite similar to what we observe using the Wigner distribution function [ 2 11, [ 2 2 ] .A number of features are worth pointing out: a) the current peak and valley for both 77 and 300 K lie near 0.12 and 0.22-0.24eV; b) the magnitudes of the current peaks are roughly the same for both temperatures; c) the NDR region has a typically smooth and welldefined current profile.The inclusion of scattering and self-consistency can be expected to alter these observa- tions; the notion of trajectories provides us with a convenient language for discussing the changes.

B. Simulation Parameters for the WDF Approach
In the WDF approach, we have endeavored to simulate an RTD with properties similar to or identical to those of experimental quantum-well devices.The parameters used, unles!otherwise noted, are: a computational b9x size of 550 A; well and barrier sizes of 50 and 30 A, respectively; for the linear bias drop, or "hot injected electron" [ l o ] , approximatio! to the bias profile, the bias linearly decreases from 30 A before the first barrier to 30 A after the second, and is constant in the rest of the computational box; the self-consistent calylations assumed that the doping w: s constant up to 30 A from the first barrier and after 30 A from the second barrier with a magnitude equal to the electron density in the boundaries (2 X 10I8 ~m -~) and the quantum well region is undoped (such spacer layers are experimentally advantageous for several reasons: they reduce the migration of impurities from the contact regions into the quantum well and they improve the peak to valley ratio of the I-V curves [ 2 2 ] , [ 2 3 ] ) ; the barriers were assumed to be 0.3 eV, and were merely added on to the bias potential profile.The effective mass of the electrons was assumed to be constant throughout the device, with a value of m* = 0.0667mo.The transient time step 6t was taken to be 1 fs.N x and N were 86 and 72, respectively.

C. Behavior of the Steady-State Trajectories
The behavior of the trajectories for an unbiased nonself-consistent RTD at 300 K shows that near the resonance momentum of 0.04 A -' (i.e., where T(k) approaches unity for the unbiased case), trajectories begin to crqss the device (note that the barriers here are thicker (50 A ) than in subsequent calculations) (Fig. 4).Their momentum increases in the quantum well and prior to the barriers, but in the barriers, the momentum decreases.Some trapped particles exist in the quantum-well region, indicated by closed loops, but, interestingly enough, other trapped particles appear to exist in the barrier regions.The latter do so in the form of "complimentary trajectories," in which the particle jumps from positive to negative momentum instantaneously, as though the particle were trapped in between impenetrable walls (this calls to mind the Gaussian wavepacket simulations, in which a portion of the wave packet was trapped inside the potential barrier and oscillated for a while before escaping).
We have calculated trajectories with and without scattering for low bias (0.04 eV) conditions for self-consistent Wigner functions.For the case without scattering at 77 K (Fig. 5(a)-(c)), the behavior of the tunneling trajectories in the well is markedly different from the unbiased case; for example, the particles appear to dwell for a while in the well before tunneling out and continuing their motion (we shall refer to this as "resonance behavior" in which the particle spends a comparatively long time in the well).This is also similar to the Gaussian wavepacket simulations, in which a portion of the wavepacket trapped in the quantum well continued to leak out slowly, a process taking on the order of 100 fs [ 181, Trajectories with momenta higher than the resonance momentum exhibit a behavior more like the unbiased case.Parenthetically, we may note that the statistics of these particle trajectories may allow us to roughly estimate the current through the device by adding up all of the individual currents for each tunneling trajectory.Performing such a calculation gives a result surprisingly close to the actual current [9] calculated by (18).Further, at higher biases, we note that some of the tunneling trajectories increase their momentum at the drain end; the difference in the incoming and outgoing kinetic energy of these trajectories roughly corresponds to the bias, so that the trajectories exhibit a "conservation of energy" behavior [ 151.
Including scattering in the previous calculation causes numerous differences: first, a number of trajectories are prevented from tunneling through, which will result in a reduction in current through the device.Second, the momentum of these tunneling trajectories is greater in the well than for the case of no scattering.This will cause the density of particles in the well to decrease, as the particles tunnel out more readily, and will thus affect the switching  time (i.e., how long the RTD takes to return to equilibrium after a sudden bias change) of the RTD, which depends on the amount of time the well takes to deplete.
We now increase the temperature to 300 K with no scattering inside the device (Fig. 6(a)-(c)).Note that nonresonant trajectories appear, that is, the tunneling appears to be "sequential, " whereas at resonance, the tunneling particles appear to experience the entire potential region.Putting in scattering obliterates the resonant trajectories, and the remaining trajectories overall resemble the unbiased one.Note that, as in the 77 K case, a number of  trajectories are prevented from tunneling.Thus we may expect that the resonant current is degraded as scattering becomes more pronounced, i.e., as the temperature increases.Further, the density in the well may also be affected as a function of temperature.

D. Effects of Scattering on Particle Density and Current Density
In light of the previous comments, the electron density and current density in the device is examined.The effects of scattering on the I-V characteristics erodes the peakto-valley ratio significantly as a function of temperature (Fig. 7).Compared to the case where scattering is absent, the current for 77 K has decreased by a factor of 2 at the resonant voltage.In the higher bias regions, where no resonant peaks in the transmission coefficient occur, the current is not as greatly affected.When self-consistency is added (Fig. 8), the current peaks are moved to a higher bias, which may have been expected in light of the selfconsistent potential profile (Fig. 9(a), (b)), in which the well potential is not as low as for the linear bias model.Consequently, a larger bias is needed to line up the Fermi level with the resonant state.Further, note that the voltage separation between the peak and valley has decreased, which duplicates experimental results.This effect of selfconsistency on I-V characteristics is similar to the experimentally well-known effect of a resistor connected in series with an Esaki tunnel diode.In the 77 K case, note that the value of the current is not uniquely given in the NDR region: the value of the current depends on whether the bias is increased through the NDR or decreased.The plateau-like structure that is apparent on the increasing bias sweep in the NDR region is the time average of an oscillating current [25] and replicates experimental observations [28]; this is the first time this effect has been seen in a numerical simulation of RTD's, and is indicative of the power of our approach to including self-consistency .Fig. 8. Current versus applied bias for 77 and 300 K with scattering and self-consistency.Note that the resonant bias has been moved to a higher voltage, and the difference between peak and valley voltages has decreased compared to the non-self-consistent calculations. (+ ) means that the current for this bias point was calculated starting with a distribution at a lower bias point; ( -) means that the current was calculated starting with a distribution at a higher bias point.For the ( + ) case, the plateau in the NDR represents a time-averaged current; all other points are steady-state currents.not be entirely due to the external circuitry, but may also be induced by the inherent characteristics of the RTD.The density of electrons in the well region first increases with bias until the resonant bias is achieved, and then decreases through the NDR region (when self-con- sistency is not included, however, the unbiased density is largest in the case of 300 K).The effect of scattering on the trajectories may help explain why the density at resonance of 300 K is smaller than for 77 K (Figs. 9 and 10): the dwell time of the "particles" in the well is reduced, i.e., the particles are scattered into momentum states whereby they may more easily tunnel out, an explanation which is inherently "trajectory" based.

E. Efects of Scattering on Transient Response
Since switching times are a measure of how long the well takes to deplete, the switching time will be more greatly affected by scattering since the resonant trajectories are greatly affected.and 12) are much larger than for non-self-consistent calculations [ 101, presumably due to the transient electric fields in the source and drain regions.For 300 K, the time evolution of the current is smooth and settles down after some 300 fs.For 77 K, the oscillations in the well region are present and long-lived.These oscillations indicate that the density in the well is oscillating, and mimics an RC-like behavior in which the density repeatedly increases and decreases.The larger mag- nitude of the initial current for 77 over 300 K is a reflection of the larger well density for the former temperature.Scattering thus seems to have the unusual effect of making the switching time of an RTD faster, primarily by dampening out the long-time oscillations in the well region.The more abrupt and smaller voltage swing in the selfconsistent I-V peak-to-valley switching also renders the self-consistent switching time to be faster, if not approximately the same, as the non-self-consistent switching times.These effects are explained by a combination of a V. COUPLING TO MONTE CARLO SIMULATIONS In the present approach to trajectories, the quantum region of the device is determined by the extento of the quantum influences, and is typically around 600 A .In the rest of the device, the quantum influences are'negligible or absent.The whole quantum region is characterized by a Wigner distribution function at each time step, and is acted on by a "quantum force" calculated from (5).The success of this simulation depends on being able to simulate the actual time-dependent boundary condition for the quantum region at each time step.
The algorithm would essentially consist of conventional Monte Carlo moves for the particles outside the quantum region, and Wigner Distribution Function (WDF) moves inside.For the WDF moves, the Wigner function must be solved at each time step with, say, a "drifted Fermi-Dirac distribution function," with the drift velocity equal to the product of the mobility and the updated electric field at the boundary of the quantum region.The forces on the particle are then calculated for the momentum and position of the particle, and its trajectory is then updated, analogous to the Monte Carlo moves outside the quantum region.This work is in the incipient stages, and we hope to report on the results in the future.

VI. CONCLUSION
We have presented a methodology to simulating onedimensional barrier problems in terms of "quantum trajectories" and have shown in what ways the trajectory concept was useful in understanding the behavior of a particular device of interest, namely, the RTD.The trajectory concept and the behavior of trajectories are useful in explaining and anticipating the behavior of the device under question, and, significantly, the methodology has application to other potential profiles, in particular, to those of interest in vacuum microelectronics, in which the drain end of the RTD device is replaced by a vacuum; work in this direction has begun.Potentially the most useful feature of the trajectory analysis is that it may provide a means to couple these quantum transport calculations with Monte Carlo calculations, which require a trajectory approach.Such a fruitful union would allow for the simulation of large, multidimensional, and time-dependent high-speed systems in which quantum tunneling processes in a small region dominate the behavior of high-performance devices.
[ 161, [ 1 7 ] ; indeed, we have shown the correspondence between Bohm trajectories and Wigner trajectories for a traveling Gaussian wave packet previously [ 1 8 ] , namely, that the Bohm trajectories are a weighted average over Wigner trajectories, that is, the velocity vB of a Bohm trajectory is given by the average velocity of the Wigner trajectories, as indicated in the following equation: l m A where hk is the momentum of the particle; U g ( t ) = h d , S / m * ; and where S is the real position-dependent phase of the wave function $(x) = R(x) exp (is@)).Wave packet simulations, which involve the numerical solution of the Schrodinger's equation, model the interaction of single particles with potential barriers, and have been studied before in the context of vacuum microelectronics Two sources of dissipation are included in the simulation.Electrons which enter the simulation region are thermalized due to scatterings outside the boundaries [ 6 ] , [ 131.
7. = ~ j" dx x3/27;(x) e-'.I r(5/2) All scattering processes are combined via the Matthiessen's Rule: 7-I = E 1 / ~~.The parameters used in the calculation of T ~, and the methods used for their tedious calculation have been given elsewhere [ 121.Self-consistency may be incorporated by simultaneously solving Poisson's equation at each time step in a transient simulation given the electron density p(x)

Fig. 1 .
Fig. 1 .Phase-space grid for the simulation of the Wigner distribution function.A priori known boundary values are given by ( 0 ), initially unknown distribution points are given by ( 0 ) ; i .e ., only the particles entering the device have a known distribution.
23)where r = 2h/St = 1.316 eV for a time step of 1 fs.Finally, we mention the matrix solution of the Wigner function.Let b be a column vector with Nx -N components, such that b((i -l ) N + j ) = f ( i , j ) .L becomes a multiplying matrix plus a boundary value vector:L f ( x , k) * A * bb,.The matrix part A is block-pentadiagonal in the SDS formulation (for CDS [9] and UDS [6],[13], the matrix is block tridiagonal), there being Nx diagonal blocks, each block being N x N in size.The diagonal blocks correspond to the action of the integral operator and part of the differential operator, and the off- diagonal blocks correspond solely to the differential operator.Standard banded matrix solving packages, such as those in IMSL, may then be used.

.Fig. 5 .
Fig. 5. (a) Unbiased self-consistent trajectories for 77 K.The number following the "k" is a trajectory label for use in (b) and (c).(b) Self-consistent trajectories at 77 K for a bias of 0.04 V without scattering.(c) Same as (b) but with scattering.The well momenta increase, and some of the lower momentum trajectories are prevented from tunneling.

Fig. 6 .
Fig. 6 .(a) Same as Fig. 5(a), but at 300 K. (b) Same as Fig. 5(b), but at 300 K. Two types of trajectories exist: "Sequential," in which the particle momentum decreases only in the barrier region.and "Resonant." in which the particle momentum is decreased throughout the quantum well region.(c) Same as Fig. 5(c), but at 300 K.The trajectories almost completely resemble the equilibrium case.All traces of "resonant" tunneling behavior have vanished.

Fig. 7 .
Fig. 7 .Current versus applied bias for various temperatures.Scattering is included, but not self-consistency.The dashed line represents the 77 K calculation without scattering.Note that the current peak-to-valley ratio has decreased significantly in the presence of scattering (T = 77 versus T = 77 (no scattering)).
If a resistor is connected with the RTD, this plateau can be made to disappear [23].These experimental results may

Fig. 9 .
Fig. 9. (a) Self-consistent potential profiles at various biases corresponding to a low bias, and the biases for the peaks and valleys of Fig. 8 for 77 K.The dashed line refers to the linear bias model for a bias drop of 0.32 V. (b) Same as (a), but for 300 K.The linear bias is for a bias drop of 0.30 V.
Fig. IO.(a) Self-consistent density versus position for 77 K.The dashed lines represents the doping profile and the band-edge discontinuity.The densities are calculated for peak and valley biases (see Fig.8), and for a low bias of 0.02 V. (b) Same as (a) but for 300 K.Note that the well density at resonance is much lower than for the 77 K case; likewise, the well density for the valley bias is much larger.

Fig.Fig. 12 .
Fig. I I .Self-consistently calculated current as a function of position and time for a sudden bias switch from peak to valley at 77 K. Scattering is included.Note the oscillations in the well region, indicating that the particles are sloshing back and forth, and that the density in the well is oscillating.