New Quantum Mechanical Simulator in ATLAS


With the advent of bandgap engineering and the high electron mobility transistor (HEMT) it became possible to take advantage of very high electron mobilities in FET like devices by separating charge carriers in the device channel from their associated dopants in adjacent buffer layers. With these devices, design tools such as S-Pisces and Blaze, have been introduced to account for variations of band-gap in drift diffusion and energy balance type device simulation[7]. These modifications give users good insight into the operation of such devices. However, these earlier simulators failed to properly account for certain quantum mechanical (QM) effects associated with the confinement of carriers in potential wells associated with heterojunctions.

Heterojunctions introduced by band-gap engineering tend to confine electrons in potential wells. These wells, typically on the order of tens of nanometers wide, cause the allowed energy states for carriers to become quantized in the direction perpendicular to the heterojunction. This quantization directly effects the distribution of carriers in the channel and buffer layers and also affects the predicted terminal currents. In addition proper accounting for these quantum effects should enable correct prediction of quantum mechanical tunneling into and out of the channel at the source and drain.

Quantum mechanical effects can also be observed under the gate of ultra-thin gate oxide MOS transistors. The effect of QM confinement in the MOS channel is manifest in a shift of the peak carrier concentration away from the silicon-oxide interface. This effect is observed as shifts in the threshold voltage, gate capacitance[6] and also changes carrier mobility. With the trend toward smaller devices with thinner gate oxides these effects will become increasingly important in the near future.

In order to properly predict QM effects as described above, Silvaco, in its continuing pioneering efforts, has introduced a QM solver into the ATLAS device simulator. This solver allows solution of Schrodinger's equation (Eq. 1) along a series of 1-D slices (in the Y direction).

Schrodinger's equation is an eigenvalue equation and as such the solutions to Schrodinger's equation appear as a sequence of eigenenergies and eigenfunctions. The eigenenergies represent the discrete energies where electrons may reside and the eigenfunctions represent the probability distribution of electrons as a function of position along each slice. To introduce self consistency into the solution of Schrodinger's and Poisson's equations the eigenfunctions are used to predict the electron distributions in space which are then put into the right hand side of Poisson's equation where the classical electron densities would normally appear. Iteration between the Schrodinger's equation and the Poisson's equation are made until convergence is obtained.

1D Schrodinger Solver

The 1D Schrodinger equation forms a tri-diagonal matrix, however since the mesh spacing and effective mass may vary throughout the solution domain the matrix is not symmetric. This complicates the eigenvalue solution since the majority of published solvers concentrate upon symmetric matrices. Further, since only the lowest dozen or so eigenvalues are required it is inefficient to solve for all the eigenvalues, an operation performed in all full matrix eigensolvers. Hence the "Sturm" sequence approach is adopted which permits non-symmetric matrix solutions and solves for each selected eigenvalue in turn. The solution of Schrodinger's equation is therefore performed in two sections, the calculation of the desired eigenvalues followed by their corresponding eigenvectors.

Eigenvalues form the roots of the characteristic polynomial of a matrix and hence the method chosen simply isolates each root in turn using a bisection algorithm. The Sturm sequence[1,2] evaluates the characteristic polynomial for a given trial eigenvalue and returns the number of roots throughout the sequence. (A root is recognized by the accompanying change of sign in the characteristic polynomial).

This, executed within the bisection algorithm, is repeated until the error tolerance on the eigenvalue is met, typically 10-12V, requiring approximately 40 iterations. Although the number of iterations seems excessive, this procedure is not expensive in computer resources since the scheme is linear with the array size and hence quick.

The initial boundaries for the bisection algorithm are provided from Gershgorin's theorem. However these limits are dynamically re-located during the solution process. After each eigenvalue has been determined the lower, minimum limit is moved to its value since all subsequent eigenvalues will have a value greater than the current eigenvalue. Also the upper limit is reduced if any trial eigenvalue produces more than the total number of desired roots. The question of degeneracy presents itself; however practical experience of the nature of the equations shows that degenerate states rarely occur in real systems.

Once an eigenvalue has been obtained it may be substituted into the Schrodinger equation and the corresponding eigenvector calculated. Even though this forms a set of linear equations the eigenvector cannot be calculated directly since direct methods return the null vector. Note: the null vector is an acceptable solution to Schrodinger's equation, however it is not the solution we are looking for. Consequently a Newton scheme is used to solve for the root, adding a correction to a trial function. This only requires 1 iteration since the equations are linear. Further, since the equations are linear the solution is not critical upon the choice of the trial vector and hence this is set to the null vector for convenience.

Once calculated, the solutions are checked to ensure that all the eigenvectors are mutually orthogonal.

Alternative Schrodinger Solver

An alternative Schrodinger solver was also implemented based upon the full matrix evaluation of the eigenvalues. This is more efficient if more than 75% of the eigenvalues are wanted. This is a general approach and although currently coded for a 1D solution could easily be modified for 2D equations. It has been used to validate the values produced by the Sturm sequence approach and has shown corellation down to 10-14, the precision for a double variable.

Carrier Densities

The carrier densities are calculated from the boundstate energies, their associated wavefunctions and a reference potential, Ef, the Fermi-level. If no current is present the Fermi-level has zero gradient, flat across the y-direction. The evaluation of this Fermi-level causes most problems associated with device simulations.

An algorithm has been written to calculate the Fermi-level by matching the sheet electron density calculated by the Schrodinger equation with that provided by the classical transport approach. This is based upon Newton's method and requires the calculation of the derivative of electron density with respect to potential.

Expanding the electron density equation it can be seen that the derivative of electron density, dn/dV, depends upon the derivative of both the eigenvalue and the wavefunction with respect to potential, dl/dV and dw/dV. First order perturbation theory shows these parameters to be -1.0 and 0.0 respectively leading to a relatively straight-forward evaluation. These parameters are also used in the self-consistent Schrodinger-Poisson solver.

Matrix Solver

The eigenvector solution requires a tri-diagonal matrix solver. This is performed by LU decomposition, an exceedingly fast procedure in 1D[3] These routines are written for general tri-diagonal matrix inversion.

Physical Issues

It is worth noting that the kinetic energy Hamiltonian is not strictly defined in a system with spacially varying effective mass. However it is reasonable to assume the electron wavefunction and its first derivative are continuous across a change in effective mass and hence one accepted modification is to move the effective mass from its normal position outside the div grad operator to within the first derivative. The Hamiltonian therefore reads: "div ( 1/mass grad ) ..." as opposed to "1/mass div(grad ..." [4].

The solver assumes a-priori that the electrons are bound. The solution domain is placed into an infinite potential well, inferred from the solution of the wavefunctions on the boundaries, all set to zero.

The solution of Schrodinger's equation within this region only depends upon the conduction band energy and the effective electron mass. Additional terms can be added, the exchange-corellation energy to account for the change from single to multi particle statistics. This cannot be included precisely since the multi-particle Schrodinger equation would have to be solved and hence an empirical form is available[5]. This adds a little complication to the solution since it is dependent upon the electron density, a parameter derived from the Schrodinger solution. Hence if the exact Schrodinger equation is desired, one would have to perform a iterative loop on Schrodinger's equation with exchange-corellation potential. In practice this is not really necessary since the exchange-corellation potential is relatively small when compared with the electron potential energy.

One further implication from the equation set is that confined electrons have a fixed momentum vector in the y direction, set by Schrodinger's equation. The wavefunctions in the x- and z-directions are assumed to be Bloch states. The Bloch state energy levels are close enough for the additional energy gained or lost from the environmental electric field to be accommodate. The electron can skip from one energy state to another, and hence electrons are allowed to move in these directions. However since the energy gained from the electric-field is not sufficient to promote an electron between the bound state energy levels in the y-direction the electron is not permitted to move. Motion in this direction would promote the electron to a state that is not present.

However, one can argue that once the inter-bound-state energy level drops below kT, the electrons do have enough energy to skip between states and motion can occur. Consequently a parameter has been included in the simulation that denotes the ratio of electrons in 2D states to the total electron density. This parameter is included in the transport equations (only in the y-direction) to account for the carrier confinement that occurs within the quantum well.

Application Examples

To demonstrate the new self-consistent Schrodinger-Poisson solver we simulated a planar HEMT device. The device is composed of an undoped GaAs substrate, capped by a 35 nanometer layer of AlGaAs with a composition fraction of 0.3 doped with a doping concentration of 1.5 X 1018 donors per cubic centimeter. In this example, a half micron long gate is separated from the source and drain by 0.5 microns on either side. The device structure and mesh around the channel are shown in Figure 1. The 1-D Schrodinger's equation is solved at periodic vertical locations in the structure. The electron concentrations calculated from the resultant eigenfunctions are then put into the 2D Poisson solver.


Figure 1. Example HEMT structure.

Figure 2 shows the simulated electron eigenenergies as a function of depth below the middle of the gate. The gate was biased at 0.5 volts. We see the first 15 eigenenergies in the potential well formed at the interface between the AlGaAs and GaAs layers. For increasing eigenenergy the spacing between subsequent eigenenergies decreases. Electrons in the lower eigenstates can be considered as bound states and cannot travel in the y direction. Electrons in the higher states are closer in energy and may contribute to current in the y direction by hopping states as discussed above.

Figure 2. Quantum states at the heterojunction.

Figure 3 shows the normalized eigenfunctions associated with the first 5 eigenenergies in Figure 2. These eigenfunctions are weighted and summed to arrive at a quantum electron density based on the self consistent solution of Schrodinger's and Poisson's equations. A comparison of the electron density as a function of depth for the classical (Boltzman's statistical) and quantum statistical predictions of the electron density is given in Figure 4. Here we can see significant differences between the classical and quantum predictions.

Figure 3. Electron distribution functions for each Quantum state.

Figure 4. Comparision of classical vs. QM electron distribution.

Conclusions and Future Work

In this article we have described a new self consistent (1-D) Schrodinger,(2D) Poisson solver implemented in the ATLAS device simulator. This solver can be used to analyze effects of quantum mechanical confinement on carrier distributions in heterojunction devices.

In this article the effects of quantum statistics are only introduced into Poisson's equation. These effects will also be introduced self consistently into the transport equations. In addition current flow in the y direction will be weighted to account for electrons residing in the bound states. These capabilities should be extendible to both drift-diffusion and energy-balance transport models.


  1. J.H. Wilkinson, "The Algebraic Eigenvalue Problem", Springer Verlag, p.300ff, 1971
  2. J. Stoer & R. Burlish, "Introduction to numerical analysis", Springer Verlag, p.385ff, 1980
  3. W.H. Press. B.P. Flannery, S.A. Teusolsky & W.T Vetterling, "Numerical Recipes in C, the Art of Scientific Computing", Cambridge University Press, p.48ff, 1986/1989
  4. R.A.Morrow & K.R.Brownstein, "Model effective-mass Hamiltonians for abrupt heterojunctions and the associated wave-function-matching conditions", Physical Review B, Vol 30, No 2, p.678, July 1984
  5. F.Stern & S.Das Sarma, "Electron energy levels in AlAs-AlGaAs heterojunctions", Physical Review B, Vol 30, No 2, p.840, July 1984
  6. F.Stern & S.Das Sarma, "Electron energy levels in AlAs-AlGaAs heterojunctions", Physical Review B, Vol 30, No 2, p.840, July 1984
  7. R. Rios and N. Arora "Determination of Ultra-Thin Gate Oxide Thickness for CMOS Structures Using Quantum Effects", International Electron Devices Meeting 1994 Technical Digest, Dec. 1994, pp. 613 - 616.
  8. ATLAS User's Manual, Version 4.0, June 1995, pp. 5-1 - 5-6.