A Comparison of Resonant Tunneling Based on Schrödinger ’ s Equation and Quantum Hydrodynamics

Smooth quantum hydrodynamic (QHD) model simulations of the current–voltage curve of a resonant tunneling diode at 300 K are compared with that predicted by the mixed-state Schrödinger equation approach. Although the resonant peak for the QHD simulation occurs at 0.15 V instead of the Schrödinger equation value of 0.2 V, there is good qualitative agreement between the current–voltage curves for the two models, including the predicted peak current values.


INTRODUCTION
Quantum semiconductor devices like resonant tunneling diodes and transistors and superlattice devices are playing an increasingly important role in state-of-the-art VLSI chips. These devices rely on quantum tunneling of charge carriers through potential barriers for their operation. Advanced microelectronic applications include multiplestate logic and memory devices and high frequency oscillators and sensors.
A fundamental approach to modeling these tunneling devices is through the mixed-state Schrödinger equation model coupled to Poisson's equation for the electrostatic potential. In three spatial dimensions, this approach involves a total of six dimensions (for position x and momentum p) plus time, so a hydrodynamic approximation, where the density, velocity, and temperature of a charge carrier are functions of x and t, offers enormous computational speedups in simulating 2D and 3D devices.
However, the (quantum) fluid dynamical approximation will only be valid if the active device length scales are not too small. As a benchmark, we will compare quantum hydrodynamic (QHD) model simulations of the currentvoltage curve of a resonant tunneling diode at 300 K with a quantum region of 200 Å with that predicted by the mixedstate Schrödinger equation approach.

THE SMOOTH QHD MODEL
Quantum transport effects including electron or hole tunneling through potential barriers and charge buildup in quantum wells can be incorporated into the hydrodynamic description of charge propagation in semiconductor devices. References [1,2] present a quantum extension of the classical hydrodynamic model to handle, in a mathematically rigorous way, the discontinuities in the classical potential energy which occur at heterojunction barriers in quantum semiconductor devices. This "smooth" QHD model is valid to all orders of É 2 =ðmT 0 l 2 Þ (where m is the effective mass of electrons or holes, T 0 is the ambient temperature, and l is a typical length scale for the problem) and to first order in the classical potential energy.
The smooth QHD equations have the same form as the classical hydrodynamic equations: where n is the electron density, u i is the velocity, m is the effective electron mass, P ij is the stress tensor, V ¼ 2ef is the potential energy, f is the electrostatic potential, e . 0 is the electronic charge, W is the energy density, and q i is the heat flux. Boltzmann's constant k B is set equal to 1. Indices i, j equal 1, 2, 3, and repeated indices are summed over. Electron scattering is modeled by the standard relaxation time approximation, with momentum and energy relaxation times t p and t w .
The stress tensor and energy density are where T is the temperature of the electron gas and the "quantum" potential V is given by ðb ¼ 1=T 0 Þ Vðb; xÞ ¼ The heat flux incorporates both classical and quantum effects and is derived in Ref. [3] by a Chapman -Enskog expansion. The quantum contribution is necessary for internal consistency of the QHD model. We will model the relaxation times t p and t w by where v s is the saturation velocity, and the coefficient k by where m n0 is the low-field electron mobility and k 0 . 0 is a phenomenological constant.
The transport equations (1) -(3) are coupled to Poisson's equation for the electrostatic potential energy where e is the dielectric constant and N is the density of donors. The total potential energy V consists of two parts, one from Poisson's equation V P and the other from the potential barriers V B : V B has a step function discontinuity at potential barriers.
To derive the stress tensor and energy density, we constructed an effective density matrix as an O(bV ) solution to the Bloch equation. Then using the momentum-shifted effective density matrix, we took moments of the quantum Liouville equation to derive the smooth QHD equations [1].
There are two contributions to the quantum potential V: the double barrier potential and the "self-consistent" electrostatic potential from Poisson's equation. Note that second derivatives of V appear in the stress tensor and energy density, which then are differenced in the smooth QHD transport equations. Thus we compute B is just computed once since it only depends on the barriers and not on the applied voltage or state variables (n, u, T, V P ). In computing V 00 P , we first use Poisson's equation to obtain V 00 P ðb; xÞ ¼ The convolution (13) can be computed efficiently using discrete Fourier transform algorithms.
In Fig. 1, we present the smooth QHD simulation of a GaAs resonant tunneling diode with Al 0.3 Ga 0.7 As double barriers at 300 K. Note the experimental signal of quantum resonance-negative differential resistance, a region of the current -voltage (I -V) curve where the current decreases as the applied voltage is increased. The barrier height is equal to 280 meV. The diode consists of n þ source (at the left) and drain (at the right) regions with the doping density N ¼ 10 18 cm 23 ; and an n channel with N ¼ 5 £ 10 15 cm 23 : The channel is 200 Å long, the barriers are 25 Å wide, and the quantum well between the barriers is 50 Å wide. Note that the device has 50 Å spacers between the barriers and the contacts.

THE SCHRÖ DINGER EQUATION MODEL
In this section, we consider another approach based on the Schrödinger picture. Here we assume that electron transport is ballistic in the spacer and double barrier zone, while highly doped regions are assumed to be at equilibrium. Therefore the highly doped regions are considered as electron reservoirs injecting electrons into the active regions with given momentum profile (equilibrium profile). More precisely, let c p be the stationary wave function corresponding to an electron beam arriving at the left boundary a with a momentum p $ 0: The energy of such a beam is where V a is the applied bias at x ¼ a: The wave function c p is then a solution of the stationary Schrödinger equation The boundary conditions satisfied by c p are transparent ones (see Refs. [7,8] for reviews on such boundary conditions and Refs. [6,5] for their mathematical analysis). They are deduced from the fact that where p b ( p ) is defined by and r p and t p are reflection and transmission amplitudes to be computed. Note that by choice, the incoming wave has amplitude one. The boundary conditions are obtained by eliminating the unknowns r p and t p . This leads to the following For electrons injected from the right boundary b with a momentum p , 0; the same type of computation leads to where The electron density and current density are given by a superposition of the densities corresponding to the individual wave functions where g( p ) is the statistics of the left reservoir for p . 0 and of the right reservoir for p , 0: In our case, since the reservoirs are assumed to be at thermal equilibrium and since the doping in the source and drain has the same value, the reservoir profile g is given by the formula where E F is the Fermi level defined implicitly by the neutrality condition for the doping density N þ in the source and drain regions: gð pÞ dp: The current density can be computed thanks to the following formula j ¼ ð þ1 21 pTð pÞgð pÞ dp ð26Þ where the transmission coefficients are given by The electrostatic potential f ¼ f e þ f s is the sum of a given external potential f e , which includes the double barriers and the applied bias, and a self-consistent potential solving In order to solve the problem numerically, we use an iterative procedure. Starting from an initial guess for the electrostatic potential, we compute the wave functions c p , deduce the electron density n, and then compute a new value of the electrostatic potential. The iterative procedure goes on until the difference between two successive potentials is below a certain threshold. Two problems arise in the implementation of this idea. The first one concerns the computation of the wave function. It is due to the existence of peaked resonances which requires a refined mesh in the p variable. These resonances depend on the electrostatic potential and we have developed a procedure in order to detect them numerically and automatically adapt the mesh size in the p variable. More precisely, we evaluate the logarithmic derivative of the transmission coefficient, which reads (for the positive momentum case) When jd log T=dpj is less than a given threshold, the mesh size is wide, while if jd log T=dpj is above the threshold, we use a small mesh step. A logarithmic derivative is necessary in order to resolve the resonance. The adaptive meshing in the p variable allows for a reduction of a factor of four in the computation time.
The second problem to be solved concerns the Schrödinger -Poisson coupling procedure. Indeed, the naive way to perform the coupling is to compute the new value of the potential f new s by solving the Poisson equation with a density given by n ¼ n½f old s where f old s is the electrostatic potential obtained at the end of the previous iteration. This idea fails because of the exponential dependence of the density on the potential, and a relaxation procedure requires a tiny relaxation coefficient which results in a very slow convergence. In order to accelerate the convergence, we use the Gummel procedure [9], which takes into account this exponential behavior. Namely, we write the Poisson equation in the equivalent form We noticed that the linear version of the Gummel method works well for applied small applied voltages for which the resonance levels are situated between the two barriers. When the applied voltage is high enough, a resonant level builds up in a quantum well which forms at the interface between the source spacer and the first barrier. In such situation, the nonlinear version of the Gummel iteration gives much better results. Other details of this method as well as transient simulations, based on the discretization of open transient boundary conditions [4], can be found in Ref. [10]. The Schrödinger equation model can be generalized to multidimensional devices for which the access zones can be considered as waveguides. Transparent boundary conditions can be derived in such situations [11] and lead to mathematically well-posed problems [12]. Numerical simulations can be found in Ref. [13].
In Fig. 2, the current -voltage curves for the resonant tunneling diode are displayed using the mixed-state Schrödinger equation model.

CONCLUSION
It is perhaps surprising that quantum hydrodynamics can produce negative differential resistance at all. It is perhaps even more surprising that, although the resonant peak for the QHD simulation occurs at 0.15 V instead of the Schrödinger equation value of 0.2 V, there is good qualitative agreement between the current -voltage curves for the two models, including the predicted peak current values.
The smooth QHD model discussed here is linearized in the potential V. Important nonlinear effects in V were dropped in the derivation of the model. As a consequence, the smooth QHD model is too "smooth": the charge density predicted by the model in the well of the resonant tunneling diode is too small as compared to kinetic theory/Schrödinger equation approaches. An open question is how to extend quantum hydrodynamics to incorporate higher order terms in the potential. We expect that better agreement between the QHD and Schrödinger equation model simulations will be obtained for quantum hydrodynamics which incorporates nonlinear terms in the potential.
Naoufel Ben Abdallah is professor of mathematics at the Université Paul Sabatier, Toulouse. His current research interests concern the modeling of charge transport in semiconductor nanostructures, including asymptotic analysis, coupling of models, and numerical simulations of quantum ballistic transport.
Olivier Pinaud is preparing a Ph.D. thesis in applied mathematics. His main research interests involve the simulation of the stationary and time-dependent Schrödinger equations with open boundary conditions and the mathematical analysis of 2D electron gas models.
Carl L. Gardner is Professor of Mathematics at Arizona State University. His current research interests lie in classical and quantum semiconductor device simulation, computational fluid dynamics, astrophysical flows, and the modelling and simulation of ion transport in the channels of cellular membranes.
Christian Ringhofer is Professor of Mathematics at Arizona State University. His current research interests include classical and quantum transport equations and moment models for semiconductor device modelling.