Full text of DOI:10.1051/epjap:2001103

Eur. Phys. J. AP 15, 117–121 (2001)
THE EUROPEAN
PHYSICAL JOURNAL
APPLIED PHYSICS
c
ꢀ EDP Sciences 2001
Simulation of 2D quantum effects in ultra-short channel
MOSFETs by a finite element method?
A. Ponceta, C. Faugeras, and M. Mouis
Laboratoire de Physique de la Mati`ereb, INSA-Lyon, 20 avenue Albert Einstein, 69621 Villeurbanne, France
Received: 20 November 2000 / Revised: 18 April 2001 / Accepted: 27 April 2001
Abstract. This paper presents a flexible numerical technique which is especially suited to analyze lateral
modulation of quantum effects in short channel MOS transistors. We discuss boundary conditions for the
Schro¨dinger equation and the impact of the finite element meshing. We show how channel length shortening
alters the sub-band structure, thus giving an evaluation of the limits of a 1D quantum approach.
PACS. 85.30.De Semiconductor-device characterization, design, and modeling – 02.70.Dh Finite-element
and Galerkin methods
1 Introduction
two dimensions [3,4] by the use of finite elements which
are especially well suited to analyze devices of arbitrary
shape, such as self- organized quantum boxes or MOS-
FETs with non uniform gate oxide.
It is now recognized that quantization must be accounted
for to analyze correctly the operation of field-effects tran-
sistors using inversion channels along an insulating or wide
bandgap semi-conducting layer (MOSFETs and MOD-
FETs). It has been shown for instance that the gate
insulator thickness of a MOSFET, as deduced from ca-
pacitance measurements, was significantly over-estimated
if quantization was neglected. Although quantization in
MOS channel was present from the origin (and usually ne-
glected), the corrections which it introduces to standard
extraction procedures are no longer negligible, due to the
reduction of gate oxide thickness in advanced technologies.
Different simulation approaches have been used in the
literature. Due to carrier confinement against the gate in-
sulator, quantum effects are usually simulated using self
consistent numerical solutions of Poisson and Schro¨dinger
equations in 1D across the channel [1,2]. This approach
may not remain valid in very short gate length MOSFETs,
where significant confining field variations are observed
between source and drain, on very short distance scales.
Besides, a strict 1D approach make it impossible to ac-
count for effects such as lateral variations of oxide thick-
ness or dopant profiles.
2 Numerical methods
2.1 Galerkin formulation of Schro¨dinger equation
Let us consider a domain Ω (of any dimensionality)
in which quantum effects need to be simulated. The
Green formula can be used to integrate Schro¨dinger equa-
tion over Ω with a test function w in a suitable func-
tional space [7]. The resulting Galerkin formulation of
Schro¨dinger equation for any valley in the silicon lattice
reads then:
!
Z
n
X
η
1 ∂ψ ∂w
m?xi ∂xi ∂xi
−
+ (E − V ) · ψ.w · ∂Ω
2
Ω
i=1
I
η
∂ψ
+
· w · ∂Γ = 0 (1)
2m?η ∂η
Γ
where Γ denotes boundaries and interfaces of Ω, and η is
the outer normal direction to this boundary. We accounted
for the ellipsoidal symmetry of the conduction band struc-
ture: m?xi is the electron effective mass in direction xi for
the valley under consideration. V is the potential energy
deduced from Poisson’s equation, and E is the eigen en-
ergy associated with wave function ψ.
However, only few attempts have been made to ac-
count for the lateral variation of quantum confinement
along the channel [3].
This paper presents results of a full-2D simulation of
quantum effects and of the influence of a gate length re-
duction down to a few ten nanometer. This work differs
from previous attempts to solve Schro¨dinger equations in
We used classical boundary conditions for Schr¨odinger
equation:
?
This paper has been presented at 3es Journ´ees Nationales
“H´et´erostructures a` semiconducteurs IV-IV”, Orsay, July 2000.
– ψ = 0 for an infinite barrier;
(2)
(3)
∂ψ
1
– continuity of
∂ψ
across interfaces;
m?η ∂η
a
e-mail: poncet@lpm.insa-lyon.fr
UMR-CNRS 5511
b
–
= 0 for reflecting boundaries.
(4)
∂η

118
The European Physical Journal Applied Physics
Correction on Electrical Potential during
Poi on rodinger teration
1E18cm-3 doping - 14nm Tox - 0->4V gate bias
All three conditions lead the boundary integrals in (1) to
vanish, making implementation in a finite element simu-
lator straightforward.
Condition (4) is especially well suited to account for
symmetries. It can be used to save CPU time by reducing
the size of the simulation domain to only half a device (in
the case of symmetric devices such as double-gate MOS-
FETs or Permeable Base Transistors, ...) or by shortening
ohmic areas.
If needed, periodic boundary conditions might also be
introduced to simulate Bloch waves. However, with such
boundary conditions, the matrices loose their band struc-
ture. This makes implementation more complex and more
CPU time consuming.
c
ca e
1
1
1
1
1
1
E-01
E-04
E-0
E-10
E-13
E-1
0
0
40
0
80
100
Iteration Number
Fig. 1. Convergence of the Schro¨dinger-Poisson decoupling
scheme for 1D simulation of a NMOS capacitor Q-V charac-
teristic.
2.2 Finite element simulation
The finite element discretization of equations (1) to (4)
is straightforward: a basis of piecewise polynomial shape
functions wi is constructed. These functions are piecewise
linear on a triangular or tetrahedral mesh, and piecewise
bi-linear on a quadrangular or hexahedral mesh. Wave
functions are projected on this basis to get a discrete
eigensystem. However, attention must be paid to the map-
ping from eigenvectors to wave functions. In this linear al-
gebraic problem, the eigenvectors do not provide directly
element method is in fact minor in comparison with its
advantages.
In another respect, a finite difference method could be
even simpler and slightly faster (as explained in the pre-
vious paragraph). We did not choose this approach how-
ever, in order to keep the latitude of simulating arbitrary
shaped devices.
the wave function values at mesh points. If numerical inte-
R
2.4 Iterative decoupling scheme
gration uses a trapezoidal rule, the matrix of
wiwj∂Ω
terms is diagonal, but is not the identity matrix.ΩThe com-
Several methods have been proposed to solve the coupled
system of Poisson and Schro¨dinger equations iteratively.
The most straightforward method consists in re-injecting
the carrier concentrations calculated from Schro¨dinger
equation in the right hand side member of Poisson equa-
tion. However, this is known to result in instabilities
which can be overcome by using under-relaxation meth-
ods. Convergence is then rather slow. A much more effi-
cient method has been proposed in [5]. It consists in devel-
oping a Newton-like method in which the Jacobian matrix
is replaced by a first order approximation. A huge amount
of mathematical work is needed to end up with approxi-
mations that are both valid and computable.
Here, we propose a new method which is altogether
simple and fast. Once electron concentration has been
computed by summing up all the wave functions contri-
butions, a “Fermi-like” potential is extracted using the
relationship which relates carrier concentration to the
electrostatic potential when a volume Boltzman statistics
holds. This quantity (which is physically meaningless in
the present situation, and must only be seen as a variable
mapping) is injected in Poisson equation, as usually done
in classical device simulators. This make Poisson equa-
tion highly non-linear but very easy to solve with a pure
Newton-Raphson procedure. The procedure is iterated un-
til convergence is reached.
ponent number i of any eigenvector is the product of the
corresponding wave function at node number i by a factor
R
Ω wi2∂Ω. On the other hand, if the integrals are computed
exactly, the above mentioned matrix is even not diagonal,
and must be first diagonalized. This is the only CPU time
overhead in comparison with finite difference schemes.
2.3 Discussion of alternative methods
An alternative way to discretize Schr¨odinger equation has
been proposed in the literature [3,4]. There, the piecewise
polynomial shape functions wi are replaced by sine func-
tions, with space coordinates splitting, in order to decor-
relate the minimum wave function wavelength from the
choice of the mesh size. However, although attractive at
first sight, this method gives full matrices while finite el-
ement and finite difference schemes lead to sparse ma-
trices. Moreover, the only boundary which can be han-
dled is condition (2), and it seems very difficult to extend
this method to arbitrarily shaped structures. Finally, the
smallest period which can be introduced with this method
is in fact the same as with our more classical discretiza-
tion scheme. This can be illustrated in the 1D case. If L
is the structure length, the smallest wavelength which can
be simulated is 2L/N in both cases, N representing the
number of terms in the sine expansion in one case and
Convergence is extremely fast, as shown in Figures 1
N +1 representing the number of mesh points in our case. and 2. Figure 1 shows the number of iterations needed to
It means that both methods actually lead to similar ma- bring the relative correction below 10−12 during the sim-
trix sizes, so that the only potential drawback of the finite ulation of the Q-V characteristic of an NMOS capacitor

A. Poncet et al.: Simulation of 2D quantum effects in ultra-short channel MOSFETs
119
Sub
b
Correction on Electrical Potential during
Poi on rodinger teration
20nm 1E18cm-3 channel - 14nm Tox - 0->4V gate bias
c
ca e
1
1
1
1
1
1
1
E-01
E-04
E-0
1
a
L=0. 2m
L=0.02 m
E-10
E-13
E-1
=20
=2
0
10
20
30
40
0
0
teration
u
er
1
1
Fig. 2. Convergence of the Schro¨dinger-Poisson decoupling
scheme for the 2D simulation of a short channel NMOS tran-
sistor gate characteristic.
Fig. 4. Eigen energy values as a function of the sub-band
index (a, b) and contribution to the total electron charge of the
corresponding wave function (c, d) for long and short channel
length (0.92 µm for a-b and 20 nm for b-d).
Here, we restrict ourselves to equilibrium conditions
(VDS = 0 V). Our method can be used then to calculate
the 2D distribution of electrons and holes, together with
the corresponding 2D potential and electric field distri-
butions, as a function of gate voltage. This gives for in-
stance indications about gate control of the channel carrier
charge. Electrical parameters such as the sub-threshold
slope and threshold voltage at low VDS can be extracted
and short-channel effects at low VDS can be evaluated:
indeed, 2D effects, resulting in a parasitic control of the
source/channel potential barrier are accounted for. More-
over, the existence of non-zero wave- functions in regions
of the device where the sub- band energy is in the band-
gap, can bring an evidence of the presence of quantum
coupling between some regions of the devices and of pos-
sible tunneling effects. However, full current-voltage char-
acteristics will be available only once this model has been
coupled to a quantum transport model. The following dis-
cussion focuses on the numerical conditions to get a phys-
ically significant solution under thermodynamic equilib-
rium conditions, with given computer resources (limited
number of mesh points).
Fig. 3. Boundary conditions for Poisson equation. For
Schro¨dinger equation, we used Neuman boundary conditions
on the outer boundaries.
(with 21 successive voltage steps). Dimensions and dop-
ing levels were taken in [6] so that we could check our
results. The oxide layer was 14 nm thick, so that high
voltages (up to 4 V) needed to be applied to the gate. In
the weak inversion regime, only a couple of iterations are
required. The error exhibits an exponential decay even at
high gate voltages. We obtained a very similar behavior
in 2D, during the simulation of the gate characteristics a
short channel NMOS transistor, with 11 successive volt-
age steps (Fig. 2). This proves that this method is efficient
in both the 1D and 2D cases.
In order to investigate how the calculated eigen en-
ergy spectrum depends on channel length, we performed
2D simulations with electrical channel lengths (L) rang-
ing from 20 nm to 1 µm, with a fixed number of mesh
points in the direction of the channel (∼50). Channel dop-
ing was 1017 cm−3, and the gate oxide was 2.5 nm thick.
The NMOS transistor was biased with 1 V on the gate
(assuming zero work function and no interface charges).
This bias voltage corresponds to inversion conditions.
3 Simulation results
When narrow channel effects are neglected, it becomes
obviously useless to solve the problem in 3D. The 2D
Schro¨dinger equation is known to be appropriate to study
quantum wire structures, i.e. carrier confinement in 2 di-
rections. The use of this 2D approach to investigate quan-
tum effects in a short channel MOSFET is less straight-
forward, because carriers are not confined laterally in
the channel. Therefore, boundary conditions and meshing
must be defined carefully.
3.1 Long channel
The boundary conditions used for Poisson and
Schro¨dinger solutions are given in Figure 3. These bound- The 2D solution provides a series of sub-bands with a
ary conditions will be discussed below with the support of minimum energy εi (corresponding each to a given wave
simulation results.
function Ψi in the (x,y) plane) and a continuum for wave

120
The European Physical Journal Applied Physics
a
b
Fig. 5. Example of wave functions obtained (a) in the S/D ohmic regions (region 1 of the energy spectrum in Fig. 4), (b)-(e)
in the first sub-band subset (wave functions with one lobe perpendicular to the gate, region 2 of the energy spectrum) and (f)
in the second sub-band sub-set (wave functions with two lobes perpendicular to the gate, region 3 of the energy spectrum).
vector kz in the z direction (along the device width). Fig-
+
N+
N
ure 4a shows the εi values obtained as a function of the
sub-band index in the long channel case (1 µm). The en-
ergy spectrum is clearly stepwise.
For moderate gate voltages, the lowest energies (re-
gion 1 of the energy spectrum) correspond to wave func-
tions which vanish outside the source and drain areas
(such as the one shown in Fig. 5a). Next come energies
which are associated to wave functions with only one lobe
Fig. 6. An example of a 2D wave function calculated in the
channel of a MOSFET (with an effective channel length of
in the direction normal to Si/SiO2 interface (such as those 20 nm).
shown in Figs. 5b to 5e). Then come energies correspond-
ing to wave functions that show two transverse lobes (such
as in Fig. 5f), and so on...
drain outer boundary (instead of periodic boundary con-
ditions) is therefore of little importance: these boundary
conditions do not disturb results in the channel (unless the
ohmic region size is really too small). However, the larger
these ohmic regions, the more eigen values will be asso-
ciated to them. It must therefore be kept in mind that if
they are taken too large, it will be necessary to calculate a
huge number of eigen energy before a sufficient number of
channel eigen energies has been obtained to get physically
significant results in the channel.
The other important result is that the eigen energy
stepwise variation is progressively smoothed (Fig. 4b) and
that the energy spectra of the different sub-sets overlap:
as shown in Figure 4d, the sub-set corresponding to wave
functions with two lobes perpendicular to the oxide (re-
gion 3) starts before the sub-set corresponding to wave
functions with only one lobe (region 2) has ended.
Therefore, the band structure degenerates into a se-
ries of sub-band sub-sets. Each sub-set corresponds to the
sub-band which would be computed using a 1D approach.
However, each sub-set contains a finite number of sub-
bands, which is equal to the number of mesh points in the
channel direction. Moreover, each wave function in this
sub-set has a similar weight in the total carrier concen-
tration (Fig. 4c). In the case shown, the 2D numerical
solution is physically meaningless, due to mesh-induced
truncation at the top of each sub-set. Indeed, as explained
below, the criterion to get physically significant results is
the mesh size. Therefore, simulation of long channel MOS-
FETs (with local non-uniformities for instance) is very de-
manding in computer resources. With 100 nodes, strong
errors are still observed at 120 nm channel length.
As explained before, this 2D approach cannot be used
for very long gates, due to the mesh-induced truncation of
the energy spectrum of each sub-set. Let us now discuss in
more details the conditions which must be fulfilled to get
physically significant results in the channel. For this pur-
pose, we consider the respective contributions to electron
3.2 Short channel
The picture is rather different for short channel MOS-
FETs. The first visible results is that the wave functions
become strongly two- dimensional, as shown in Figure 6.
However, it should be noted that, as in the long chan- concentration of the first and last computed sub-bands
nel case, the wave functions associated to S/D ohmic re- in sub-set Si, which we name here αlow(Si) and αhigh(Si).
gions vanish in the channel while channel wave function The ratio αhigh(Si)/αlow(Si) reflects the error due to mesh
vanish in the source and drain. The fact that we chose zero induced truncation in Si. Figure 7 shows the variation of
normal derivative boundary conditions at the source and αhigh(S1)/αlow(S1) with mesh size for different numbers

A. Poncet et al.: Simulation of 2D quantum effects in ultra-short channel MOSFETs
121
1,E
1,E
0
0
In addition, note that if several sub-set are occupied,
the condition αhigh(S1)/αlow(S2) ꢀ 1 indicates that all
significant sub-bands have actually been computed, and
gives a validity criterium for the whole computation.
int r
t
L=
50 n
m
1,E 00
1,E-0
L= 2 nm
00 n
intra
t
1,E-0
L=20 nm 50 n
1,E-0
4 Conclusions
L=2 nm
00 n
1,E-0
It can be concluded that 2D numerical solution of the
coupled Poisson and Schr¨odinger equations may be ob-
tained with reasonable computing resources (here, simu-
lations of the NMOS structure were done on a desktop
computer). Provided a sufficiently small mesh size has
been chosen, physically significant solutions can be ob-
1,E-10
0,1
1
10
100
Mesh Size (nm)
top / bottom, 1st sub-set (50 mesh nodes laterally)
top / bottom, 1st sub-set (100 mesh nodes laterally)
tained for nanoscale gate lengths, in the range where a
top 1st sub-set / bottom nd sub-set (50 mesh nodes
)
1D approach is actually inadequate.
Fig. 7. Evaluation of the mesh truncation error to the calcula-
tion of the electron charge as a function of mesh size, channel
length and number of nodes.
References
of mesh nodes (50 and 100) in the lateral direction. This
ratio depends on mesh size rather than on gate length or
mesh number. This was expected from first order eval-
uation. Indeed, the minimum simulated wavelength is di-
rectly correlated to mesh size and the highest energy wave
function in a given subset is expected to have the shortest
“wavelength” in the lateral direction. It is therefore nat-
ural to get a correlation between the residual error in a
given subset and the lateral mesh size. However, speak-
ing of wavelengths here is an image as the 2D functions
are not really periodical and it is interesting to note that
this trend remains valid in strongly 2D situations. Given
a number of nodes, the diagram of Figure 7 provides the
maximum channel length which can be simulated with a
given accuracy on electron distribution.
1. M. Moglestue, J. Appl. Phys. 59, 3175 (1986).
2. C. Raynaud, J.-L. Autran, P. Masson, M. Bidaud, A.
Poncet, Mat. Res. Soc. Symp. Proc., Vol. 592, Boston, 1999.
3. A. Abramo, A. Cardin, L. Selmi, E. Sangiorgi, IEEE Trans.
Electron Devices 47, 1858 (2000).
4. A. Abou-Elnour, K. Schuenemann, J. Appl. Phys. 74, 3273
(1993).
5. A. Trellakis, A.T. Galick, A. Pacelli, U. Ravaioli. J. Appl.
Phys. 81, 7880 (1997).
6. M.J. Van Dort, P.H. Woerlee, A.J. Walker, Solid-State Elec-
tron. 37, 411 (1994).
7. P. Ciarlet, The finite element method for elliptic problems,
Studies in Mathematics and its applications (North Holland,
1978), Vol. 4.
To access this journal online:
www.edpsciences.org