Three-dimensional Spectral Solution of Schriidinger Equation

We present a fast and robust method for the full-band solution of Schr6dinger’s equation on a grid, with the goal of achieving a more complete description of high energy states and realistic temperatures. Using Fast Fourier Transforms, Schr6dinger’s equation in the one band approximation can be expressed as an iterative eigenvalue problem for arbitrary shapes of the conduction band. The resulting eigenvalue problem can then be solved using Krylov subspace methods as Arnoldi iteration. We demonstrate the algorithm by presenting an example concerning non-parabolic effects in an ultra-small Metal-Oxide-Semiconductor quantum cavity at room-temperature. For this structure, we show that the non-parabolicity of the conduction band results in a significant lowering of high-energy electronic states.


INTRODUCTION
As electronic device dimensions approach nanometer scale, quantum effects are expected to dominate their electronic properties. Consequently there is considerable interest in the efficient numerical simulation of such structures, not only for exploring novel device architectures but also for maintaining reliability for reduced present-day devices. In this paper we present a fast and robust method for the numerical solution of the full-band Schrodinger equation, with the goal of achieving a more complete description of high energy states and realistic temperatures. We then demonstrate the use of the algorithm by examining nonparabolic effects on energy levels of an ultra-small Metal-Oxide-Semiconductor (MOS) quantum cavity at room-temperature.
In order to understand why full-band algorithms are needed, we have to consider that many quantum simulations assume that the kinetic energy of carriers has a quadratic momentum dependence E(p)=pZ/2m. However, this model is only correct for the lowest carrier energies. At higher energies, the effective mass m* (E) becomes energy dependent, or, with other words, the kinetic * Tel." (503)-613-7642, Fax: (503)-613-8950, e-mail: Alexandros.Trellakis@intel.com *Corresponding author. Tel.: (217)-244-5765, Fax: (217)-244-4333, e-mail: ravaioli@uiuc.edu 342 A. TRELLAKIS AND U. RAVAIOLI energy contains cubic and higher order terms in momentum p. If we increase the energy even further, even this approximation breaks down and a full-band description of the kinetic energy as a function of momentum E--E(p) becomes necessary.
Therefore, the accurate quantum simulation of high energy quantum states poses a challenge that has not been addressed yet in previous studies. Atomistic models might be a possible solution [1], however, the computational costs involved in these models are still prohibitively high for generalpurpose device simulation. Therefore, in order to simulate larger quantum structures it seems to be highly desirable to maintain a continuum description of the semiconductor that is corrected for high energy band structure effects. It is the goal of this paper to outline the implementation of such a method, and to illustrate its feasibility by simulating a silicon MOS quantum structure.

COMPUTATIONAL METHOD
In order to obtain this continuum full-band algorithm, we restrict carrier motion to one band (one band approximation) and ignore interband transitions. In this case, the kinetic energy of carriers in the band is described by a function E(p). This E(p) might not be available in analytical form but only as an interpolated table, as it is common for the output of band structure calculations. Nevertheless, independent from the origin of E(p) the usual quantization rules apply, and we may replace the momentum p by the momentum operator to obtain the Hamiltonian. Thus, the Schr6dinger equation for envelope wavefunctions describing carrier motion under the influence of an external potential V(x) can be written as if V(x) does not vary too much within a lattice constant a, and the dimensions of the device structures under consideration are much larger than a as well [2].
Finding a solution for this equation is not as straightforward as for parabolic energy bands. This fact becomes quite obvious, if we consider that now the kinetic energy operator E(O) may contain the momentum up to very high and possibly infinite orders. Since is proportional to the gradient vector in the position representation, this means that we have to solve an equally high order differential equation. Numerical solution strategies based on a finite differences or finite elements need to Taylor expand E(O) at least up to a order n+ in p in order to compute E(O) up to order n.
While it is possible to use such high order schemes [3,4], this approach can become inefficient for large n. Also, depending on the functional form of E(p), the convergence radius of its Taylor series might be too small for polynomial approximations. Finally, E(p) may only be available in tabulated form, making it impossible to derive high order polynomial approximations. While there exist solutions for these problems, they usually involve matrix representations of the kinetic energy operator that are quite complicated and possibly dense [5]. Therefore, considering these issues and the fact that the kinetic energy operator E(O) is diagonal in momentum space, a spectral solution approach using Fourier transforms seems to be a natural choice for this problem. If we choose to perform our calculation in position space, we may insert Fourier transforms (FT) around the kinetic energy operator, which then reduces to simple multiplication in momentum space. We then obtain for Schr6dinger's equation In fact, one could view space. Alternatively, we may also have expressed Schr6dinger's equation in momentum space by inserting Fourier transforms around the potential energy. This would give us a Hamiltonian similar to the ones employed in density functional calculations, except that there the kinetic energy is quadratic in the electron momentum due to the atomistic nature of these models [6]. However, for both position and momentum space formulation, the numerically difficult problem of applying a high order kinetic energy operator on a wavefunction has now been reduced to the more manageable one of calculating Fourier transforms.

IMPLEMENTATION
Even though electronic structure calculations are usually formulated in momentum space [6], it is not clear whether a position space or a momentum space formulation leads to a numerically superior algorithm for our continuum model. However, since all numerical routines used by the authors assume position space, it seems prudent to maintain compatibility by implementing the position space Eq. (2). In order to obtain an efficient numerical implementation it is mandatory to use Fast Fourier Transforms, which can be calculated using only O(N log N) computational steps on grid with N nodes, as compared to O(N2) for other methods.
Since the point-wise multiplication of two functions involves only O(N) steps, examining Eq. (2) shows that we can calculate the action of the Hamilton operator on a wavefunction with O(N log N) steps. Having this fast vector-matrix multiplication available, fast iterative eigenvalue solvers like the Arnoldi method become immediately the candidates of choice. In fact, since the Hamilton matrix is only implicitly available (unless one tries to calculate it), iterative eigenvalue solvers are the only practical choice. We use an Arnoldi solver with Chebyshev preconditioner to solve this eigenvalue problem [7].
The use of Fast Fourier Transforms imposes stringent restrictions in the choice of the computational grid. Most important, the grid lines have to be equidistant in each coordinate direction, and additionally, for many Fast Fourier Transforms commonly available, the number of grid lines has to be a power of 2 in each direction as well. These requirements can result in some unexpected difficulties. A too coarse grid samples only parts of momentum space and provides insufficient spatial resolution, which is expected. However, a too fine one results in an computational first Brillouin zone that is larger than the first Brillouin zone of the lattice. This poses the question how to deal with the high momentum computational states that are introduced this way. Mapping these states to physically equivalent low momentum states using the periodicity of the physical momentum space is not a viable option, since this approach would introduce spurious low energy modes into the solution. Instead, for the computation we assign very high energies to these modes, thereby eliminating their influence on lower energy bound states.
Finally, we would like to mention that the fullband approach as implemented here is not a good substitute for a finite difference based solver in the case of parabolic energy bands, since the requirement of equidistant grids is very restrictive. Also, the Fourier based matrix-vector multiplication used in the full-band solver scales as O(N log N) for N grid nodes, while the banded matrix-vector multiplication used in the conventional solver scales as O(N). Therefore, for large parabolic problems it is always faster to use a solver based on finite difference or finite element discretization. This is especially true in the one-dimensional case, where tridiagonal eigenvalue solvers are far superior to any iterative method.

EXAMPLE APPLICATION
As an example application of our full-band algorithm, we examine the influence of conduction 344 A. TRELLAKIS AND U. RAVAIOLI band valley non-parabolicity for the MOS quantum cavity shown in Figure 1. Unlike quantum dot floating gates for flash-memories [8,9], a quantum box of undoped silicon with dimensions 8 x 10 x 8 nm is here in direct contact with a 30nm thick epilayer of undoped silicon below, allowing for the exchange of electrons. The box is surrounded at the sides and the top with SiO2, with the oxide thickness in the normal direction being 2nm directly above the box and 10 nm everywhere else. The surface of the oxide is completely metallized by a 5 nm thick layer of A1. Vertical confinement towards the substrate is controlled by the electric field under the gate in conjunction with a highly p-doped deep substrate. Here, substrate doping is chosen as NA--1018cm -3 in order to improve 5  Vxc the exchange correlation potential in the local density approximation (see [10,11] for details concerning the numerical methods).
The non-parabolicity of the conduction band valleys in silicon is modeled using the following parameterizations described by Jacoboni and Reggiani [12], zml Here, Et, kt and mt describe the energy, wave number and effective mass in the transversal direction of a valley in the Silicon conduction band. Similarly, El, kl and ml describe the corresponding quantities in the longitudinal direction. a is a parameter with the dimension of an inverse energy that describes the non-parabolicity of the valley. Following the suggestion in [12], we choose a 0.7 eVfor our calculation.
No special treatment of the silicon-oxide interface was performed. Since at these interfaces the conduction band structure changes, a more precise model than the one described above might appear to be necessary in order to obtain correct results.
However, the bound states of the quantum dot barely penetrate into the oxide due to the more than 3eV high conduction band step at the interface. Therefore, for this application we may safely assume that the kinetic energy E(p) is material-independent without incurring major inaccuracies. An extension of the full-band method that also allows for position dependent E(p) will be discussed in a future publication [13].
The resulting energy spectra for a parabolic and a non-parabolic valley simulation are shown in Figures 2a and b. We find that the overall structure of the two spectra is very similar. The total occupation of the cavity is about five electrons in both cases, with about 1.5 electrons residing in the ground state, 2.5 electrons lying in a cluster of 5 states around 300 meV, and the states above 300meV being mostly unoccupied (the fractionate occupations seen here are to be understood in the context of thermodynamic equilibrium distributions as the Fermi distribution).
However, we find that due to the conduction band non-parabolicity all states except the ground state are shifted towards lower energies with respect to the bottom of the conduction band. This shift is not too surprising if we consider that in Eq. (6) the effective electron mass increases with increasing transversal momentum kt within each conduction band valley. Since the states in all three energy ladders are occupied, as exemplified by the cluster of states located at 300eV in Figure 2, a significant influence of non-parabolicity on occupied electron states can be expected. Despite this energy shift, the total occupation of the cavity does not change much, since the Fermi level has been shifted downwards as well, and the total number of electrons is mainly determined by the strength of the electric field under the gate in the cavity. Instead, the electrons are redistributed over the energy spectrum, with the ground state occupation being decreased and the occupation of the higher energy states being increased. Transforms to express the full-band Hamilton operator as an iterative eigenvalue problem for arbitrary shapes of the energy band. The resulting eigenvalue problem is then solved using a standard iterative solver as Arnoldi iteration. The applicability on the method is demonstrated on the example of a MOS quantum cavity at room-temperature. For this structure, we have shown that the nonparabolicity of the conduction band results in a significant lowering of high-energy electronic states.