Wigner Function Methods in Modeling of Switching in Resonant Tunneling Devices

Issues associated with modeling of quantum devices within the framework of the transient Wigner equation are addressed. Of particular importance is the structure of the Wigner function, hysteresis, and device switching time, whose value is determined by the large signal device properties.

INTRODUCTION Numerical simulation of quantum devices is a noble pursuit, with a primary advantage in tran- sient studies.For while the operating principles of quantum devices are often ascertained with relatively simple approximations and analytical studies provide insight into the tunneling times, the details associated with non-uniform fields, scatter- ing, very short time scales, noise contributions, and large signal device-circuit interactions, require simulation.Given this, what is the best approach to studying transients?
The best approach should be ascertained prag- matically; constrained by how much can be done in an acceptable period of time, and the knowledge base.In the past the approach we have taken [1] has included: (i) studies of the quantum hydro- dynamic equations, (ii) the density matrix in the coordinate representation, and (iii) the Wigner distribution function [2].Of the three the weakest is the quantum hydrodynamics derived as quan- tum corrections [3].
The discussion of this document concerns the efficacy of Wigner equation simulations for the design of transient quantum devices; specifically, high frequency low phase-noise clocks, the stage set by recent experiments [4].The ingredients of this study require circuit (i.e., transmission line) equations, transient Wigner/Poisson equation, dissipation, and a noise model.We have simulated the large signal transient switching and recovery of a resonant tunneling device, the transient behavior of single barrier devices subject to a large signal sinusoidal input, the dc IV characteristics of double barrier and classical NIN devices, the role and presence of hysteresis, and have addressed numerical issues.We have required the algorithm to deliver current continuity.The issues raised by these studies are discussed below.

THE PHYSICAL MODEL AND THE DEVICE EQUATIONS
The basic quantum transport equation is the Wigner equation [5]: The potential energy in Eq. ( 1) consists of two contributions, the barrier/well configuration and the potential energy arising from Poisson's equation.In most cases the contributions from Poisson's equation were treated classically as a term 7xVpoIsso N Vkfw(k,x).As originally dis- cussed [5], the integration in Eq. ( 1) is over all space, and dissipation was not considered.But device issues include dissipation [6], finite size and open boundaries [7], as considered below.
In every case studied, the barriers were square with height V0 and width A. This permitted an analytical integration of the Wigner integral, with the following result used in all of our studies: In Eq. ( 2) x0 is the central position of the barrier, A the barrier width and V0 the barrier height.The center position, height and width of each barrier were defined separately.
The Wigner function is scattered by the barrier and the analytical form of this scattering as given by Eq. ( 2), indicates that for large values of A relative to L (none of which were considered here) and/or large values of Xo-X, the Wigner integral behaves as a delta function in kx-k'x.Local interactions are then enhanced.Further structure arises for a pair of barriers separated by a distance 2, leading to a term cos ((kx k'x)5).
In the studies presented below the device is broken into a classical and quantum region.The bounding reservoirs and cladding regions are treated classically, with the central region repre- senting the quantum mechanical region.Within this framework the Wigner integral is multiplied by a modulating function that is equal to unity within the "quantum region" and zero elsewhere.
The Poisson contribution is treated classically throughout the entire structure.For calculations discussed below in which the device length is 200nm, the quantum region is at least 120Am long.For each of the calculations discussed below the transition between the quantum and classical regions takes place within the heavily doped region.
What about dissipation?First we recognize that the Wigner function is not the Boltzmann dis- tribution, and the generic argument involving occupancy of the quantum states is less obvious.
The most detailed approach is that of [6].But the tack taken by many is to relax the Wigner func- tion.Here several approaches have been taken.One discussed in [8] and [9], as well as elsewhere, is to assume a relaxation time approximation: The relaxation time approximation in the above form leads to source and sink terms in the continuity equation.To avoid sink terms others [10] have multiplied the equilibrium distribution function by the ratio of the non-equilibrium carrier density to the equilibrium carrier density.
We have done both, but will confine ourselves to that represented by Eq. ( 3).
The question of interest is what is f0(k,x)?We must remember that unlike many early discussions in the theory of metals where it is often assumed that the electric field under equilibrium is near or approximately zero, in the case of electron devices where there are barriers, variations in background density, the local electric field can be in the range of mega-volts/cm.The form of the equilibrium distribution function is dependent upon the model used to connect current at the open boundary and must represent the spatially dependent distribu- tion associated with barriers, scattering and selfconsistency.The models used in this study have evolved.Initially we assumed a displaced Fermi Dirac distribution function on the boundaries.While, in principle, the displaced momentum should be zero, or at least zero within five or six significant figures, it wasn't.The boundary condi- tion was then changed to one in which the normal derivative (with respect to position) of the dis- tribution function was set to zero.This provided the requisite zero current conditions.Further, to enhance the possibilities of flat-band open bound- ary conditions the relaxation time in the vicinity of the boundaries was set at least an order of magnitude smaller than elsewhere in the devices.We note that initially (first time step) we assume a Fermi-Dirac distribution everywhere.Now, consider some details in differencing.The topic most workers tend to focus on is the term containing the spatial derivative of the distribution function.In the earliest studies specific attention was given to the importance of difference schemes representing particles coming in with negative/ positive values of momentum.Specifically, the spatial derivative for carriers entering with positive values of momentum is displaced by one spatial index, from particles coming in with negative values of momentum [11].When this condition was entered into our numerical schemes, the equilibrium distribution function was asymmetric with respect to momentum.Clearly, a non- physical result!To deal with this issue we introduced central differences, with the derivative defined halfway between grid points and expressed in terms of the differences between distribution functions at adjacent grid points.Under equilibrium conditions the distribution function is symmetric.
The next issue, with respect to differencing is concerned with continuity, a problem of particu- lar importance when dealing with device issues.Starting from the analytical expressions for the Wigner equation, most difference techniques will not provide numerically accurate expressions of continuity unless the number of grid points becomes excessive.The way around this is to write a Wigner equation in difference form that supports the analytical forms of the Wigner equation in its most common integro-differential form.We as well as others [7] workers have done this: and while the results appear promising; numerically based Brilllouin zone effects are introduced.
Another issue, thought by some to be based upon physical constraints, is the structure of the Wigner function.Is the 'numerical' noise that appears in most numerical simulations physical?
We cannot yet answer this question, but some of our results discussed below suggest that this 'numerical' noise is not physical.
The simulations discussed are for structures whose nominal length is 200 nm, although several studies were performed for structures whose length was 100 nm.The 100 nm structures were not regarded as representative of any device because conditions at the boundaries did not satisfy flat band constraints.Rather they are of interest only because of the insight they provided in the structure of the Wigner function.Indeed a careful examination of most device simulations with the Wigner equation for structures less than 200nm, do not show flat band conditions.First some numbers in a typical device simulation the number of spatial grid points, Ns, was 150, and the num- ber of momentum grid points, NM, was 150.The number of momentum grid points we need to reach a value of energy equal to the Fermi energy is readily calculated.(For parameters suitable for GaAs and for a Fermi energy equal to 40 meV the Fermi wave vector is reached with an index greater than "9".Thus the computation takes place for momentum values well in excess of the Fermi energy.)The matrix we dealt with was constructed from a square tridiagonal N x Ns with NM NM sub-matrices both along and straddling the diag- onal.The computational matrix includes matrices for Poisson's equation and the boundaries.In performing the simulations we discretize both the position and momentum indices.How we discre- tize them is arbitrary.We have used the following argument.If in equilibrium we assume periodic boundary conditions for a wave function then kx=n:r/L, where the device length is 2L.The differential increment in momentum space 2r/ Device Length has been the starting point of all of our numerical simulations with the smallest value being Ak :r/Device Length.
The key nature of numerical noise in the Wigner distribution function was displayed in a sequence of calculations of bias dependent distribution functions for a resonant tunneling structure that displayed the usual current drop back.This struc- ture was 200nm long, with a nominal doping of 101S/cm 3 everywhere, except in the center, where the background density was a lower 1015/cm 3.For this calculation the momentum grid spacing was half of the nominal value discussed above.The Wigner function at bias levels below the threshold for negative differential conductivity did not display any numerical noise.As the bias level increased the distribution functions retained their apparently smooth structure until the bias was below/above but near threshold.For very high values of bias the numerical noise in Wigner function was again absent.As reported earlier, this calculation did not show hysteresis in the dc characteristics [2].In a supplementary calculation with a device length of 100nm, and the same number of spatial and momentum grid points, the reduction in length resulted in an increased range for the momentum wave vector and a corresponding increase in the spatial resolution.While this calculation did not display flat-band conditions at the boundaries, the calculation, which consisted of two barriers and an interior quantum well, displayed negative conduc- tance, but did not sustain an obvious numerical noise in the Wignev distribution.Subsequent calculations with single barriers within a 200 nm structure, with increased momentum range and increased spatial resolution showed a corresponding decrease in Wigner distribution noise.

TRANSIENT DEVICE SIMULATIONS
In the design of high frequency clocks we are interested in the repetitive nature of the device response.That is how long does it take the device to relax prior to a re-interrogation.Is this number a function of bias and of switching range?If the answer is yes then simple descriptions of switching times are irrelevant.What we did was subject the RTD to a time dependent change in voltage.Figure shows the change in voltage.Remember- ing that the anode is negative with respect to the cathode and that negative values of voltage re- sult in movement across the region of negative differential conductivity, we start from a point just below that value of voltage that results in a drop in current, to a point well beyond the region of negative conductance.The voltage is both FIGURE Voltage pulse for switching.
decreased and increased at a controlled rate, as is the dwell time.
The response of the device to this change in voltage is shown in Figure 2, which displays the voltage change as well as the current response at one of the contracts.Because we have dissipation within the device current continuity does not imply current that is independent of position throughout the device.The current displayed in Figure 2 consists of two components, a particle current and a displacement current.The displacement current contributions arise from local changes in the electric field in the heavily doped regions at the boundaries.When the system relaxes toward equilibrium these displacement current contribu- tions vanish.You will notice that the dwell time for this calculation is approximately 200fs at which time the voltage was lowered to the value just below threshold for negative differential conductance.Several things should be noted.First, while we did not see dc hysteresis in the current versus voltage relation here the currents at the same values of voltage are different for increasing and decreasing voltages.This result is expected!Also note that the transient current, at least on this time scale appears to react almost instantaneously to the change in voltage.This is likely to be an artifact of the calculation.This result teaches that for this bias change the device does not relax within 200fs, the duration of the voltage dwell time.The relaxation time of this structure is more closely obtained from the results of Figure 3, where the transient voltage pulse differs from that of Figure 2 in the dwell time.It would appear from this calculation that the current relaxes in under 0.75 ps.These results speak to the switching time of the device, but not to whether this device will result in sustained oscillations at the shorter period excitation.Here as in the case of relaxation oscillations in space charge dominated classical devices, sustained oscillations are likely to depend upon the ability of the circuit to quench any incipient space charge effects.
The switching properties of the RTD are expected to depend not only on the dwell time but also on the voltage swings.For example in another calculation (Fig. 4), starting from the same initial voltage level but going down to a voltage level of approximately 400meV, the carriers appear to relax to steady state in a time close to 500 fs, which is shorter than that realized above.Also notice that the rise time for this calculation was shorter than for the Figure 3 calculation, resulting in a higher displacement   current contribution.This more rapid relaxation is likely to be a consequence of more rapidly removing the charge that built up in the quantum well.These calculations indicate that the measure of a quantum device switching time is also function of the device environment.
While the key issue is the interaction of the device with its external circuit we have begun to study the interaction of the device with a sinusoidal input.Others have done this as well [11].To date we have done this only for a single barrier device, and one for which there are an insufficient number of grid points.But it is because it highlights the advantages of performing simulations in which there isn't a contact preference.These results are displayed in Figure 5, which displays the Wigner function, because it highlights the advantages of performing simula- tions in which there isn't a contact preference.These results are displayed in Figure 5, which displays the Wigner function at two instants of time for a structure subject to a sinusoidal change in voltage with a period of approximately classically, with the quantum features determining the actual distribution of charge, along with quan- tum interference effects as represented by the momentum dependent structure in the Wigner function.

FIGURE 2
FIGURE 2 Voltage pulse and subsequent current response.FIGURE3As in Figure2, but with a longer dwell time.

FIGURE 3
FIGURE 2 Voltage pulse and subsequent current response.FIGURE3As in Figure2, but with a longer dwell time.

FIGURE 4
FIGURE 4 Transient response to a larger change in voltage.

FIGURE 5
FIGURE 5 Response of a single barrier structure to sinusoi-