- Email: [email protected]
Compact subthreshold current and capacitance modeling of short-channel double-gate MOSFETs
Mathematical and Computer Modelling 51 (2010) 901–907
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Compact subthreshold current and capacitance modeling of short-channel double-gate MOSFETs U. Monga ∗ , H. Børli, T.A. Fjeldly UniK – University Graduate Center and Norwegian University of Science and Technology, N-2027 Kjeller, Norway
article
info
Keywords: Device modeling Nanoscale MOSFET Current modeling Subthreshold slope Capacitance modeling Conformal mapping
abstract A precise two-dimensional subthreshold current and capacitance modeling of shortchannel, nanoscale double-gate MOSFETs is presented. The model covers a wide range of geometries and material combinations. The subthreshold model is based on conformal mapping techniques. The results are in excellent agreement with numerical simulations. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction In the near future, double-gate (DG) MOSFETs will become one of the potential candidates to take over the classical single-gate MOSFET. [1]. Due to their near-ideal subthreshold slope and negligible junction capacitance, DG MOSFETs are also being considered as very promising candidates for subthreshold applications [2]. Compact models are necessary for future CAD applications, and thus we here present Spice compatible compact models for drain current, subthreshold slope and capacitance in the subthreshold regime. The drain current is modeled using the drift-diffusion theory, although in the nanoscale regime, the drain current will have the characteristics of both driftdiffusion and ballistic transport. It has been shown that the drift-diffusion model with constant mobility gives rise to similar characteristics as those of ballistic models [3], and thus the drift-diffusion model can still be used in the ballistic/quasiballistic regime for compact modeling. In the subthreshold regime, the maximum charge density is at the center of the device and goes to minimum at the gate contacts, thus hampering the electronic confinement; also for the body thickness of 12 nm and above, the structural confinement is quite small and there is a small difference between the charge sheet densities calculated classically and quantum mechanically [3,4]. Thus we have neglected the quantization effects for the present modeling scheme. We have assumed the following device specifications for the DG MOSFET (see Fig. 1): silicon substrate thickness tsi = 12 nm, insulator thickness tox = 1.6 nm, relative dielectric permittivity of the insulator εox = 7.0, and p-type silicon body doping Na = 1015 cm−3 . To simplify the calculation, we replaced the insulator thickness tox by the electrostatically equivalent 0 = tox εsi /εox where εox = 11.8 is the dielectric permittivity of silicon. Hence the effective thickness of the thickness tox 0 extended body is H = tsi + 2tox . Idealized, equipotential drain and source Schottky contacts with work function 4.17 eV are used in order to facilitate the model verification against the numerical simulations (Silvaco Atlas), thus avoiding depletion regions in the contacts. As gate material, we selected a near mid-gap metal with the work function 4.53 eV (corresponding to that of molybdenum). 2. Drain current modeling In the subthreshold regime, the electrostatics in the device body is dominated by the inter-electrode capacitive coupling between the electrodes. This can be described by a two-dimensional Laplace equation using the potentials of the
∗
Corresponding author. E-mail address: [email protected] (U. Monga).
0895-7177/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2009.08.043
902
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 1. Schematic view of the DG MOSFET.
electrodes as boundary conditions. An analytical solution is obtained by performing a conformal mapping of the device cross-section from the normal (x, y)-plane to the upper half-plane of the complex (u, iv)-plane defined by the following Schwartz–Christoffel transformation [5–7], z = x + iy =
L F (k, w) 2 K (k)
.
(1)
Here F (k, w) is the complex elliptic integral of the first kind, K (k) = F (k, 1) is the corresponding complete elliptic integral, and w = u + iv . The modulus k is a constant between 0 and 1, which is determined by the geometric ratio L/H [5]. The drain current based on the drift-diffusion theory can be expressed as: Id = −qµn W
ZZ
n (x, y)
dVF (x) dx
dydx.
(2)
The double integral runs over the length and height of the silicon body. W is the device width, µn is the electron mobility, n(x, y) is the electron distribution, and VF (x) is the quasi-Fermi potential, which is assumed to be constant over any given cross-section perpendicular to the x-axis. Invoking current continuity along the channel and separating (2) into a coordinatedependent and a VF -dependent part, we obtain [4]:
h
Id =
Wqµn Vth 1 − exp − Vds V
i
th
R L/2
−L/2 R t +t 0 si ox 0 tox
.
dx
(3)
h i n2 i exp ϕ(x,y) dy Na Vth
Here Vth is the thermal voltage, ni is the intrinsic carrier concentration in silicon, Vds is the drain–source voltage, and ϕ(x, y) is the potential distribution determined from the inter-electrode coupling [3,4]. The integrals in (3) do not have simple analytical solutions because of the functional form of ϕ(x, y). From (3) we observe that the current is determined from the inverse of the electron sheet density along the channel. Thus, to obtain a compact model, we look for a functional form of the inverse charge sheet density along the channel. Fig. 2 shows the variation of the numerically computed inverse charge sheet density along the channel; it is evident that the inverse charge sheet density can be approximated by a normal distribution with a mean value at the top of the barrier. This occurs due the fact that the top of the barrier carries the least amount of charge and thus the inverse charge sheet density has its maximum at that position. The distribution is given as: 1 ns (x)
=
1 nsm
(x − xm )2 exp − σ2
(4)
where nsm is the charge sheet density at the top of the barrier, and σ defines the position where the charge has increased by a ratio of e compared to that of the top of the barrier. To find nsm , we need a potential distribution along the y-axis. A parabolic potential distribution in the y-direction has been shown to be a good approximation in subthreshold [6], i.e.,
"
ϕ(x, y) ≈ ϕc (x) 1 − 1 −
2y H
2 #
+ Vgs − VFB .
(5)
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
903
Fig. 2. Numerically computed inverse charge sheet density, compared with modeled results.
where ϕc = ϕ(x, H /2) − Vgs + VFB is the difference between the potential at position x on the D-S symmetry line and that of the gate-silicon interface of the extended body, Vgs is the gate-source potential, and VFB is the flat-band voltage of the gate. As the current is determined by the inverse of charge sheet density, which has its maximum at top of the barrier, any error at the source/drain contacts can be neglected for all practical cases. [4], also there is some underestimation of inverse charge density at the top of the barrier, when using (4), this leads to an overestimation of the charge and the current, but all of these errors are within acceptable limits. Using (5) together with Boltzmann statistics for the electron distribution, nsm can be expressed as nsm =
n2i N
e
ϕm Vth
"
0 +t tox si
Z
exp − 0 tox
ϕc (x)
Vth
1−
2y
2 # dy
H
(6)
where the integral has an analytical solution in terms of error functions, the value of above integral is given as:
nsm =
n2i Na
ϕm
H
e Vth
o n √ √ ϕc√ (xm )tsi π Vth Erf H Vth . √ 2 ϕc (xm )
(7)
The inter-electrode potential distribution is given by the following analytical expression [4],
1 − ku 1 + ku 1−u 1+u π − tan−1 − tan−1 + tan−1 + tan−1 1 kv kv v v ϕ (u, v) = π −1 1 − ku −1 1 − u −1 1 + u −1 1 + ku + V tan − tan + V + V − tan tan ( ) bi bi ds kv v kv v Vgs − VFB
. (8)
Here Vbi is the built-in voltage of the source and drain electrodes. Note that in (8) we have neglected the effect of the finite insulator thickness in the four corners of the extended device body. However, in the calculations of the drain current, √a simple correction for this effect has been implemented [6]. The S-D symmetry line maps into a semi-circle of radius 1 / k √ around the origin in the (u, iv )-plane. By setting v = k−1 − u2 in (8) we have the potential profile along the S-D symmetry line. The barrier maximum along the S-D symmetry line is then found by differentiating (8) with respect to u, um =
(1 + k)Vds , 2k(2Vbi + Vds − 2Vgs + 2VFB )
r vm =
1 k
− u2m .
(9)
The negative sign is added to um , so as to make its direction consistent with direction of source. Using the mapping function of (1), we obtain the position xm . The corresponding potential is ϕm = ϕ (um ). The above analysis slightly overestimates the potential in the channel and thus it also overestimates the current. This occurs due to the fact that near the source and drain contacts, the potential is relatively large, allowing a significant amount of electrons to accumulate. Thus it is necessary, to consider their electrostatic effect even under subthreshold conditions. In the subthreshold regime this effect is included simply by adjusting the boundary conditions at source and drain as follows [4]: Vcs = Vbi −
Vth 2
EO ES
2 Vcd = Vbi + Vds −
Vth 2
EO
2 (10)
ES
√
where Vcs and Vcd are new contact potentials at source and drain boundary, Eo is given as E0 = 2qNc Vth /εsi , Nc is total electron concentration at each contact, ES and ED represent the total surface electric field at the source and drain contacts.
904
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 3. Comparison of modeled (symbols) and numerically simulated (solid curves) Ids −Vgs characteristics. The upper curve for each L represents Vds = 0.5 V and the lower represents Vds = 0.1 V.
When applying (4) and (7) in the integrals of (3), and making necessary correction as described earlier, we obtain the following expression for the drain current:
h
Wqµn Vth 1 − exp − Vds
Id = √
V
i
th
π
2
σ Erf
L−2xm 2σ
nsm
L+2xm 2σ
+ Erf
.
(11)
Finally we need the value of σ , as discussed before. At x = xm + σ the charge increases by a factor of e compared to that at top of the barrier. Again resorting to parabolic functions, we use an equation similar to (5) for the S-D symmetry line. σ is then given as:
p σ =
0 +t )−ϕ ϕ(xm + σ , tox si m (L/2 − xm ) (Vbi − ϕm )
(12)
0 where the potential ϕ(xm + σ , tox + tsi /2) is obtained from the following expression:
√ n2i Na
0 +t /2) ϕ(xm +σ ,tox si
H
Vth
e
√ 0 +t /2)t ϕ(xm +σ ,tox si si √ π Vth Erf H Vth p = ensm . 0 + t /2) 2 ϕ(xm + σ , tox si
(13)
Fig. 3 shows the modeled current compared with numerical simulations for DG MOSFETS of lengths L = 25 nm and 50 nm. 3. Subthreshold slope modeling Using (2) and (7) at x = 0, the drain current in the subthreshold regime can also be expressed as:
Ids = q
n2i Na
µn W
dVF (x) dx
e
− VV(x) tf
ϕSD (0)
e
x =0
Vth
H
n √ o √ c (0)tsi π Vth Erf Hϕ√ Vth . √ 2 ϕc (0)
(14)
Since the current is constant along the x-direction, it can be evaluated at any point along x and the results will be same. Thus to evaluate the subthreshold slope, we have chosen x = 0. Since this point is at the gate-to-gate symmetry line, we obtain the model for subthreshold slope in the following comprehensible, compact form: S=
ln(10)Ids
(15)
dIds /dVgs .
Differentiating (14) with respect to Vgs results in:
∂ Ids mIds = ∂ Vgs Vth
1+
Vth ∂ F /∂ Vgs mF
(16)
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
905
Fig. 4. Comparison of modeled (symbols) and numerically simulated (solid curves) subthreshold slope with channel at Vds = 0.1 V and Vgs = VFB .
where H F =
o n √ √ c (0)tsi π Vth Erf Hϕ√ Vth . √ 2 ϕc (0)
(17)
and m=1+
2
π
−1
√
−1
k − tan
tan
1
√
k
.
(18)
From the above we get an expression of subthreshold slope as: S=
Vth ln(10)
m 1+
Vth ∂ F /∂ Vgs mF
.
(19)
This model encompasses the effects of the device dimensions, bias voltages and source–drain junction doping on the subthreshold slope. Fig. 4 shows the variation of subthreshold slope with channel length and body thickness at Vds = 0.1 V and Vgs = VFB . It is evident from Fig. 4 that subthreshold slope decreases with increase in channel length, and hence the device shows ideal characteristics around 50 nm. The improvement in subthreshold characteristics is the result of decreased short-channel effects (SCEs), as we move to longer devices. 4. Corner correction For realistic gate insulator thicknesses, (8) is a good approximation for most of the body interior, but tends to fail close to the four corners of the body. The capacitance modeling requires a more accurate description of the potential in the corner areas. Therefore a more accurate description of the potential profile across the insulator gaps of the four corners must be included in the boundary conditions of the Laplace equation. The potential distribution across the insulator gaps is analyzed considering a single corner with infinitely long gate contact and infinitely deep source or drain contacts [8]. The potential distribution in this region can be modeled by performing a new conformal mapping defined by the following Schwartz–Christoffel transformation [7]
∂ z1c = ∂w1c
√ w1c − 1 , w1c
or z1c = x1c + iy1c =
0 nhp 2tox
π
w1c − 1 − tan−1
p
w1c − 1
i
o +1 .
(20)
Here, the subscript 1c indicates that this is a one-corner transform. An analytical solution to the potential distribution in the (u1c , v1c )-plane can then be found [7],
ϕL1c (u1c , v1c ) =
1 2
Vgs − VFB + VS /D −
1
π
Vgs − VFB − VS /D tan−1
u1c
v1c
.
(21)
Mapping this potential back to the (x1c , y1c )-plane we obtain an estimate of the potential profile across the four oxide gaps in the four corners of the DG device. With some simplifications these potential profiles are applied as boundary conditions of the four-corner transform to obtain an improved potential distribution in the corner regions [2].
906
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 5. Cross-section of DG and GAA device with equivalent circuit of charge conserving intrinsic capacitances, CXY . Extrinsic overlap capacitances, Cextr , and boundaries between the intrinsic and extrinsic capacitances are indicated.
5. Intrinsic capacitance modeling The intrinsic capacitances are calculated from charges determined by the vertical electric displacement field on the electrode surfaces. This displacement field can be derived from the inter-electrode electrostatics. Charge conservation requires that the charges associated with the inter-electrode coupling add up to zero charge when summed over all electrodes. For symmetric gate biasing, the DG MOSFET can be considered a three-terminal device. Imposing charge conservation [9], such a device will have nine capacitances CXY as indicated in the equivalent circuit in Fig. 5 [10]. Of these, four are independent. The self- (X = Y ) and trans-capacitances (X 6= Y ) reflect the change of the charge on electrode X with a small variation in the voltage on electrode Y . Using the approximate potential distribution in (8), we derive the following analytical expression for the electrode charges associated with the inter-electrode coupling [8],
∂ϕL du ∂v v→0 zmin umin umax iεsi (u − 1)(ku + 1) u+1 u−1 . = VG ln + VS ln + VD ln π (u + 1)(ku − 1) ku + 1 ku − 1 umin
QX = εsi
Z
zmax
E⊥ dz = iεsi
Z
umax
(22)
Here, VG = Vgs − VFB , VS = Vbi , VD = Vds + Vbi . The limits of integration are the appropriate dimensions of the electrodes 0 over which the intrinsic charges are distributed. For the drain and source electrodes, the integration runs from y = tox to 0 tsi + tox for x = ±L/2, respectively, or between the corresponding coordinates on the u-axis in the four-corner (u, v)-plane. The latter are obtained from the transformation in (1). An analysis of the electric field distributions in the corner regions [8] indicates that the charges located on the gate electrode inside the boundary of the extrinsic overlapping capacitance as indicated in Fig. 5, correspond to field lines that terminate on the sides of the source/drain electrodes. Therefore, in order to preserve intrinsic total charge neutrality, these charges should be excluded and assigned to the extrinsic capacitances. From the one-corner analysis, we find that the integration for the intrinsic gate charge should run between x = −L/2 + x0 and L/2 − x0 for the two gates, where 0 x0 = 0.339tox [8], or between the corresponding coordinates along the u-axis in the four-corner (u, v)-plane. The capacitances are found by taking the derivatives of the electrode charges with respect to the various electrode potentials, i.e., CXY = ±dQX /dV Y , where the plus sign is used for X = Y (self-capacitances) and the minus sign is used for X 6= Y (trans-capacitances). For the inter-electrode coupling, QX is a linear function in VY , which means that all the capacitances are bias independent. Owing to symmetry and charge conservation, we find that CGS = CGD = CSG = CDG = −CGG /2, CDS = CSD , and CSS = CDD . The modeled intrinsic inter-electrode capacitances at different device lengths are compared to numerical simulations (Silvaco Atlas) in Fig. 6. We observe that the source and drain self-capacitances,. CSS and CDD , are not noticeably affected by the change in gate length. As expected we also observe that the capacitances associated with the gate increase for increasing gate lengths until they flatten out at L ≈ 30–35 nm. This concurs with CDS and CSD approaching zero, and can be considered the long-channel limit above which short-channel effects are no longer significant. 6. Extrinsic overlap capacitance The extrinsic overlap capacitances which are indicated in Fig. 5 can be modeled by the one-corner conformal mapping procedure, which is described above. Using ϕL1c of (21), the charge associated with the gate contact, QG1c , is derived following
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
907
Fig. 6. Modeled intrinsic inter-electrode capacitances (symbols) as a function of DG device length. Oxide and substrate thicknesses are held constant, tsi = 12 nm and tox = 1.6 nm. Model is bias independent. Numerical simulations were carried out at Vds = 0 V and Vgs = −0.2 V and are plotted with solid lines.
the same procedure as in (22) QG1c = iεsi
Z
umax umin
dϕL1c
VG − VS / D umax = i ε ln . si dv v→0 π umin
(23)
Here umin and umax correspond to the boundary between the intrinsic and extrinsic capacitance as indicated in Fig. 5, and the overlap length between the gate and the source/drain contact, respectively. The extrinsic overlapping capacitance associated with one of the corners is then found by differentiating QG1c in (23) with respect to VG . 7. Conclusion We have developed a precise, compact two-dimensional subthreshold current and capacitance model for nanoscale DG MOSFETs. The two-dimensional modeling is based on a precise conformal mapping of the inter-electrode capacitive coupling, including a special treatment of the corner regions. The current and capacitances calculated by the present method show excellent agreement with numerical simulations (Silvaco Atlas). The scaling of the modeled subthreshold slope with channel length, also shows good agreement with simulated results. Acknowledgements This work was supported by the European Commission under contract no. 506844 (SINANO) and the Norwegian Research Council under contract No. 159559/130 (SMIDA). We acknowledge the donation of TCAD tools from Silvaco. References [1] The International Technology Roadmap for Semiconductors, 2003. Available: http://public.itrs.net. [Online]. [2] Jae-Joon Kim, Kaushik Roy, Double Gate-MOSFET subthreshold circuit for ultralow power applications, IEEE Tran. Electron Devices 51 (9) (2004). [3] T.A. Fjeldly, S. Kolberg, B. Iniguez, Precise 2D compact modeling of nanoscale DG MOSFETs based on conformal mapping techniques, in: Proc. NSTINanotech 2006, Boston, MA, May 7-11, 2006, vol. 3, pp. 668–673. [4] S. Kolberg, Modeling of electrostatics and drain current in nanoscale double-gate MOSFETs, Ph.D. Thesis, Norwegian University of Science and Technology, 2007. [5] S. Kolberg, T.A. Fjeldly, 2D modeling of nanoscale double gate silicon-on-insulator MOSFETs using conformal mapping, Phys. Scripta T126 (2006) 54–60. [6] H. Børli, S. Kolberg, T.A. Fjeldly, B. Iniguez, Precise modeling framework for short-channel double-gate and gate-all-around MOSFETs, IEEE Trans. Electron Devices 55 (10) (2008). [7] E. Weber, Electromagnetic Fields. 1, Wiley, New York, 1950. [8] H. Børli, K. Vinkenes, T.A. Fjeldly, Physics based capacitance modeling of short-channel double-gate MOSFETs, Phys. Status Solidi 1–4 (2008). [9] D.E. Ward, R.W. Dutton, A charge-oriented model for MOS transistor capacitances, IEEE J Solid-State Circuits 13 (1978) 703–708. [10] M. Nawaz, T.A. Fjeldly, A new charge conserving capacitance model for GaAs MESFETs, IEEE Trans. Electron Devices 44 (1997) 1813–1.
Contents lists available at ScienceDirect
Mathematical and Computer Modelling journal homepage: www.elsevier.com/locate/mcm
Compact subthreshold current and capacitance modeling of short-channel double-gate MOSFETs U. Monga ∗ , H. Børli, T.A. Fjeldly UniK – University Graduate Center and Norwegian University of Science and Technology, N-2027 Kjeller, Norway
article
info
Keywords: Device modeling Nanoscale MOSFET Current modeling Subthreshold slope Capacitance modeling Conformal mapping
abstract A precise two-dimensional subthreshold current and capacitance modeling of shortchannel, nanoscale double-gate MOSFETs is presented. The model covers a wide range of geometries and material combinations. The subthreshold model is based on conformal mapping techniques. The results are in excellent agreement with numerical simulations. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction In the near future, double-gate (DG) MOSFETs will become one of the potential candidates to take over the classical single-gate MOSFET. [1]. Due to their near-ideal subthreshold slope and negligible junction capacitance, DG MOSFETs are also being considered as very promising candidates for subthreshold applications [2]. Compact models are necessary for future CAD applications, and thus we here present Spice compatible compact models for drain current, subthreshold slope and capacitance in the subthreshold regime. The drain current is modeled using the drift-diffusion theory, although in the nanoscale regime, the drain current will have the characteristics of both driftdiffusion and ballistic transport. It has been shown that the drift-diffusion model with constant mobility gives rise to similar characteristics as those of ballistic models [3], and thus the drift-diffusion model can still be used in the ballistic/quasiballistic regime for compact modeling. In the subthreshold regime, the maximum charge density is at the center of the device and goes to minimum at the gate contacts, thus hampering the electronic confinement; also for the body thickness of 12 nm and above, the structural confinement is quite small and there is a small difference between the charge sheet densities calculated classically and quantum mechanically [3,4]. Thus we have neglected the quantization effects for the present modeling scheme. We have assumed the following device specifications for the DG MOSFET (see Fig. 1): silicon substrate thickness tsi = 12 nm, insulator thickness tox = 1.6 nm, relative dielectric permittivity of the insulator εox = 7.0, and p-type silicon body doping Na = 1015 cm−3 . To simplify the calculation, we replaced the insulator thickness tox by the electrostatically equivalent 0 = tox εsi /εox where εox = 11.8 is the dielectric permittivity of silicon. Hence the effective thickness of the thickness tox 0 extended body is H = tsi + 2tox . Idealized, equipotential drain and source Schottky contacts with work function 4.17 eV are used in order to facilitate the model verification against the numerical simulations (Silvaco Atlas), thus avoiding depletion regions in the contacts. As gate material, we selected a near mid-gap metal with the work function 4.53 eV (corresponding to that of molybdenum). 2. Drain current modeling In the subthreshold regime, the electrostatics in the device body is dominated by the inter-electrode capacitive coupling between the electrodes. This can be described by a two-dimensional Laplace equation using the potentials of the
∗
Corresponding author. E-mail address: [email protected] (U. Monga).
0895-7177/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.mcm.2009.08.043
902
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 1. Schematic view of the DG MOSFET.
electrodes as boundary conditions. An analytical solution is obtained by performing a conformal mapping of the device cross-section from the normal (x, y)-plane to the upper half-plane of the complex (u, iv)-plane defined by the following Schwartz–Christoffel transformation [5–7], z = x + iy =
L F (k, w) 2 K (k)
.
(1)
Here F (k, w) is the complex elliptic integral of the first kind, K (k) = F (k, 1) is the corresponding complete elliptic integral, and w = u + iv . The modulus k is a constant between 0 and 1, which is determined by the geometric ratio L/H [5]. The drain current based on the drift-diffusion theory can be expressed as: Id = −qµn W
ZZ
n (x, y)
dVF (x) dx
dydx.
(2)
The double integral runs over the length and height of the silicon body. W is the device width, µn is the electron mobility, n(x, y) is the electron distribution, and VF (x) is the quasi-Fermi potential, which is assumed to be constant over any given cross-section perpendicular to the x-axis. Invoking current continuity along the channel and separating (2) into a coordinatedependent and a VF -dependent part, we obtain [4]:
h
Id =
Wqµn Vth 1 − exp − Vds V
i
th
R L/2
−L/2 R t +t 0 si ox 0 tox
.
dx
(3)
h i n2 i exp ϕ(x,y) dy Na Vth
Here Vth is the thermal voltage, ni is the intrinsic carrier concentration in silicon, Vds is the drain–source voltage, and ϕ(x, y) is the potential distribution determined from the inter-electrode coupling [3,4]. The integrals in (3) do not have simple analytical solutions because of the functional form of ϕ(x, y). From (3) we observe that the current is determined from the inverse of the electron sheet density along the channel. Thus, to obtain a compact model, we look for a functional form of the inverse charge sheet density along the channel. Fig. 2 shows the variation of the numerically computed inverse charge sheet density along the channel; it is evident that the inverse charge sheet density can be approximated by a normal distribution with a mean value at the top of the barrier. This occurs due the fact that the top of the barrier carries the least amount of charge and thus the inverse charge sheet density has its maximum at that position. The distribution is given as: 1 ns (x)
=
1 nsm
(x − xm )2 exp − σ2
(4)
where nsm is the charge sheet density at the top of the barrier, and σ defines the position where the charge has increased by a ratio of e compared to that of the top of the barrier. To find nsm , we need a potential distribution along the y-axis. A parabolic potential distribution in the y-direction has been shown to be a good approximation in subthreshold [6], i.e.,
"
ϕ(x, y) ≈ ϕc (x) 1 − 1 −
2y H
2 #
+ Vgs − VFB .
(5)
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
903
Fig. 2. Numerically computed inverse charge sheet density, compared with modeled results.
where ϕc = ϕ(x, H /2) − Vgs + VFB is the difference between the potential at position x on the D-S symmetry line and that of the gate-silicon interface of the extended body, Vgs is the gate-source potential, and VFB is the flat-band voltage of the gate. As the current is determined by the inverse of charge sheet density, which has its maximum at top of the barrier, any error at the source/drain contacts can be neglected for all practical cases. [4], also there is some underestimation of inverse charge density at the top of the barrier, when using (4), this leads to an overestimation of the charge and the current, but all of these errors are within acceptable limits. Using (5) together with Boltzmann statistics for the electron distribution, nsm can be expressed as nsm =
n2i N
e
ϕm Vth
"
0 +t tox si
Z
exp − 0 tox
ϕc (x)
Vth
1−
2y
2 # dy
H
(6)
where the integral has an analytical solution in terms of error functions, the value of above integral is given as:
nsm =
n2i Na
ϕm
H
e Vth
o n √ √ ϕc√ (xm )tsi π Vth Erf H Vth . √ 2 ϕc (xm )
(7)
The inter-electrode potential distribution is given by the following analytical expression [4],
1 − ku 1 + ku 1−u 1+u π − tan−1 − tan−1 + tan−1 + tan−1 1 kv kv v v ϕ (u, v) = π −1 1 − ku −1 1 − u −1 1 + u −1 1 + ku + V tan − tan + V + V − tan tan ( ) bi bi ds kv v kv v Vgs − VFB
. (8)
Here Vbi is the built-in voltage of the source and drain electrodes. Note that in (8) we have neglected the effect of the finite insulator thickness in the four corners of the extended device body. However, in the calculations of the drain current, √a simple correction for this effect has been implemented [6]. The S-D symmetry line maps into a semi-circle of radius 1 / k √ around the origin in the (u, iv )-plane. By setting v = k−1 − u2 in (8) we have the potential profile along the S-D symmetry line. The barrier maximum along the S-D symmetry line is then found by differentiating (8) with respect to u, um =
(1 + k)Vds , 2k(2Vbi + Vds − 2Vgs + 2VFB )
r vm =
1 k
− u2m .
(9)
The negative sign is added to um , so as to make its direction consistent with direction of source. Using the mapping function of (1), we obtain the position xm . The corresponding potential is ϕm = ϕ (um ). The above analysis slightly overestimates the potential in the channel and thus it also overestimates the current. This occurs due to the fact that near the source and drain contacts, the potential is relatively large, allowing a significant amount of electrons to accumulate. Thus it is necessary, to consider their electrostatic effect even under subthreshold conditions. In the subthreshold regime this effect is included simply by adjusting the boundary conditions at source and drain as follows [4]: Vcs = Vbi −
Vth 2
EO ES
2 Vcd = Vbi + Vds −
Vth 2
EO
2 (10)
ES
√
where Vcs and Vcd are new contact potentials at source and drain boundary, Eo is given as E0 = 2qNc Vth /εsi , Nc is total electron concentration at each contact, ES and ED represent the total surface electric field at the source and drain contacts.
904
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 3. Comparison of modeled (symbols) and numerically simulated (solid curves) Ids −Vgs characteristics. The upper curve for each L represents Vds = 0.5 V and the lower represents Vds = 0.1 V.
When applying (4) and (7) in the integrals of (3), and making necessary correction as described earlier, we obtain the following expression for the drain current:
h
Wqµn Vth 1 − exp − Vds
Id = √
V
i
th
π
2
σ Erf
L−2xm 2σ
nsm
L+2xm 2σ
+ Erf
.
(11)
Finally we need the value of σ , as discussed before. At x = xm + σ the charge increases by a factor of e compared to that at top of the barrier. Again resorting to parabolic functions, we use an equation similar to (5) for the S-D symmetry line. σ is then given as:
p σ =
0 +t )−ϕ ϕ(xm + σ , tox si m (L/2 − xm ) (Vbi − ϕm )
(12)
0 where the potential ϕ(xm + σ , tox + tsi /2) is obtained from the following expression:
√ n2i Na
0 +t /2) ϕ(xm +σ ,tox si
H
Vth
e
√ 0 +t /2)t ϕ(xm +σ ,tox si si √ π Vth Erf H Vth p = ensm . 0 + t /2) 2 ϕ(xm + σ , tox si
(13)
Fig. 3 shows the modeled current compared with numerical simulations for DG MOSFETS of lengths L = 25 nm and 50 nm. 3. Subthreshold slope modeling Using (2) and (7) at x = 0, the drain current in the subthreshold regime can also be expressed as:
Ids = q
n2i Na
µn W
dVF (x) dx
e
− VV(x) tf
ϕSD (0)
e
x =0
Vth
H
n √ o √ c (0)tsi π Vth Erf Hϕ√ Vth . √ 2 ϕc (0)
(14)
Since the current is constant along the x-direction, it can be evaluated at any point along x and the results will be same. Thus to evaluate the subthreshold slope, we have chosen x = 0. Since this point is at the gate-to-gate symmetry line, we obtain the model for subthreshold slope in the following comprehensible, compact form: S=
ln(10)Ids
(15)
dIds /dVgs .
Differentiating (14) with respect to Vgs results in:
∂ Ids mIds = ∂ Vgs Vth
1+
Vth ∂ F /∂ Vgs mF
(16)
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
905
Fig. 4. Comparison of modeled (symbols) and numerically simulated (solid curves) subthreshold slope with channel at Vds = 0.1 V and Vgs = VFB .
where H F =
o n √ √ c (0)tsi π Vth Erf Hϕ√ Vth . √ 2 ϕc (0)
(17)
and m=1+
2
π
−1
√
−1
k − tan
tan
1
√
k
.
(18)
From the above we get an expression of subthreshold slope as: S=
Vth ln(10)
m 1+
Vth ∂ F /∂ Vgs mF
.
(19)
This model encompasses the effects of the device dimensions, bias voltages and source–drain junction doping on the subthreshold slope. Fig. 4 shows the variation of subthreshold slope with channel length and body thickness at Vds = 0.1 V and Vgs = VFB . It is evident from Fig. 4 that subthreshold slope decreases with increase in channel length, and hence the device shows ideal characteristics around 50 nm. The improvement in subthreshold characteristics is the result of decreased short-channel effects (SCEs), as we move to longer devices. 4. Corner correction For realistic gate insulator thicknesses, (8) is a good approximation for most of the body interior, but tends to fail close to the four corners of the body. The capacitance modeling requires a more accurate description of the potential in the corner areas. Therefore a more accurate description of the potential profile across the insulator gaps of the four corners must be included in the boundary conditions of the Laplace equation. The potential distribution across the insulator gaps is analyzed considering a single corner with infinitely long gate contact and infinitely deep source or drain contacts [8]. The potential distribution in this region can be modeled by performing a new conformal mapping defined by the following Schwartz–Christoffel transformation [7]
∂ z1c = ∂w1c
√ w1c − 1 , w1c
or z1c = x1c + iy1c =
0 nhp 2tox
π
w1c − 1 − tan−1
p
w1c − 1
i
o +1 .
(20)
Here, the subscript 1c indicates that this is a one-corner transform. An analytical solution to the potential distribution in the (u1c , v1c )-plane can then be found [7],
ϕL1c (u1c , v1c ) =
1 2
Vgs − VFB + VS /D −
1
π
Vgs − VFB − VS /D tan−1
u1c
v1c
.
(21)
Mapping this potential back to the (x1c , y1c )-plane we obtain an estimate of the potential profile across the four oxide gaps in the four corners of the DG device. With some simplifications these potential profiles are applied as boundary conditions of the four-corner transform to obtain an improved potential distribution in the corner regions [2].
906
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
Fig. 5. Cross-section of DG and GAA device with equivalent circuit of charge conserving intrinsic capacitances, CXY . Extrinsic overlap capacitances, Cextr , and boundaries between the intrinsic and extrinsic capacitances are indicated.
5. Intrinsic capacitance modeling The intrinsic capacitances are calculated from charges determined by the vertical electric displacement field on the electrode surfaces. This displacement field can be derived from the inter-electrode electrostatics. Charge conservation requires that the charges associated with the inter-electrode coupling add up to zero charge when summed over all electrodes. For symmetric gate biasing, the DG MOSFET can be considered a three-terminal device. Imposing charge conservation [9], such a device will have nine capacitances CXY as indicated in the equivalent circuit in Fig. 5 [10]. Of these, four are independent. The self- (X = Y ) and trans-capacitances (X 6= Y ) reflect the change of the charge on electrode X with a small variation in the voltage on electrode Y . Using the approximate potential distribution in (8), we derive the following analytical expression for the electrode charges associated with the inter-electrode coupling [8],
∂ϕL du ∂v v→0 zmin umin umax iεsi (u − 1)(ku + 1) u+1 u−1 . = VG ln + VS ln + VD ln π (u + 1)(ku − 1) ku + 1 ku − 1 umin
QX = εsi
Z
zmax
E⊥ dz = iεsi
Z
umax
(22)
Here, VG = Vgs − VFB , VS = Vbi , VD = Vds + Vbi . The limits of integration are the appropriate dimensions of the electrodes 0 over which the intrinsic charges are distributed. For the drain and source electrodes, the integration runs from y = tox to 0 tsi + tox for x = ±L/2, respectively, or between the corresponding coordinates on the u-axis in the four-corner (u, v)-plane. The latter are obtained from the transformation in (1). An analysis of the electric field distributions in the corner regions [8] indicates that the charges located on the gate electrode inside the boundary of the extrinsic overlapping capacitance as indicated in Fig. 5, correspond to field lines that terminate on the sides of the source/drain electrodes. Therefore, in order to preserve intrinsic total charge neutrality, these charges should be excluded and assigned to the extrinsic capacitances. From the one-corner analysis, we find that the integration for the intrinsic gate charge should run between x = −L/2 + x0 and L/2 − x0 for the two gates, where 0 x0 = 0.339tox [8], or between the corresponding coordinates along the u-axis in the four-corner (u, v)-plane. The capacitances are found by taking the derivatives of the electrode charges with respect to the various electrode potentials, i.e., CXY = ±dQX /dV Y , where the plus sign is used for X = Y (self-capacitances) and the minus sign is used for X 6= Y (trans-capacitances). For the inter-electrode coupling, QX is a linear function in VY , which means that all the capacitances are bias independent. Owing to symmetry and charge conservation, we find that CGS = CGD = CSG = CDG = −CGG /2, CDS = CSD , and CSS = CDD . The modeled intrinsic inter-electrode capacitances at different device lengths are compared to numerical simulations (Silvaco Atlas) in Fig. 6. We observe that the source and drain self-capacitances,. CSS and CDD , are not noticeably affected by the change in gate length. As expected we also observe that the capacitances associated with the gate increase for increasing gate lengths until they flatten out at L ≈ 30–35 nm. This concurs with CDS and CSD approaching zero, and can be considered the long-channel limit above which short-channel effects are no longer significant. 6. Extrinsic overlap capacitance The extrinsic overlap capacitances which are indicated in Fig. 5 can be modeled by the one-corner conformal mapping procedure, which is described above. Using ϕL1c of (21), the charge associated with the gate contact, QG1c , is derived following
U. Monga et al. / Mathematical and Computer Modelling 51 (2010) 901–907
907
Fig. 6. Modeled intrinsic inter-electrode capacitances (symbols) as a function of DG device length. Oxide and substrate thicknesses are held constant, tsi = 12 nm and tox = 1.6 nm. Model is bias independent. Numerical simulations were carried out at Vds = 0 V and Vgs = −0.2 V and are plotted with solid lines.
the same procedure as in (22) QG1c = iεsi
Z
umax umin
dϕL1c
VG − VS / D umax = i ε ln . si dv v→0 π umin
(23)
Here umin and umax correspond to the boundary between the intrinsic and extrinsic capacitance as indicated in Fig. 5, and the overlap length between the gate and the source/drain contact, respectively. The extrinsic overlapping capacitance associated with one of the corners is then found by differentiating QG1c in (23) with respect to VG . 7. Conclusion We have developed a precise, compact two-dimensional subthreshold current and capacitance model for nanoscale DG MOSFETs. The two-dimensional modeling is based on a precise conformal mapping of the inter-electrode capacitive coupling, including a special treatment of the corner regions. The current and capacitances calculated by the present method show excellent agreement with numerical simulations (Silvaco Atlas). The scaling of the modeled subthreshold slope with channel length, also shows good agreement with simulated results. Acknowledgements This work was supported by the European Commission under contract no. 506844 (SINANO) and the Norwegian Research Council under contract No. 159559/130 (SMIDA). We acknowledge the donation of TCAD tools from Silvaco. References [1] The International Technology Roadmap for Semiconductors, 2003. Available: http://public.itrs.net. [Online]. [2] Jae-Joon Kim, Kaushik Roy, Double Gate-MOSFET subthreshold circuit for ultralow power applications, IEEE Tran. Electron Devices 51 (9) (2004). [3] T.A. Fjeldly, S. Kolberg, B. Iniguez, Precise 2D compact modeling of nanoscale DG MOSFETs based on conformal mapping techniques, in: Proc. NSTINanotech 2006, Boston, MA, May 7-11, 2006, vol. 3, pp. 668–673. [4] S. Kolberg, Modeling of electrostatics and drain current in nanoscale double-gate MOSFETs, Ph.D. Thesis, Norwegian University of Science and Technology, 2007. [5] S. Kolberg, T.A. Fjeldly, 2D modeling of nanoscale double gate silicon-on-insulator MOSFETs using conformal mapping, Phys. Scripta T126 (2006) 54–60. [6] H. Børli, S. Kolberg, T.A. Fjeldly, B. Iniguez, Precise modeling framework for short-channel double-gate and gate-all-around MOSFETs, IEEE Trans. Electron Devices 55 (10) (2008). [7] E. Weber, Electromagnetic Fields. 1, Wiley, New York, 1950. [8] H. Børli, K. Vinkenes, T.A. Fjeldly, Physics based capacitance modeling of short-channel double-gate MOSFETs, Phys. Status Solidi 1–4 (2008). [9] D.E. Ward, R.W. Dutton, A charge-oriented model for MOS transistor capacitances, IEEE J Solid-State Circuits 13 (1978) 703–708. [10] M. Nawaz, T.A. Fjeldly, A new charge conserving capacitance model for GaAs MESFETs, IEEE Trans. Electron Devices 44 (1997) 1813–1.