- Email: [email protected]
Transport coefficients for a silicon hydrodynamic model extracted from inhomogeneous Monte-Carlo calculations
Solid-State ElectronicsVol. 35, No. 4. pp. 561-569, 1992 Printed in Great Britain. All rights reserved
0038-1101/92 $5.00+ 0.00 Copyright © 1992 Pergamon Press plc
T R A N S P O R T COEFFICIENTS FOR A SILICON HYDRODYNAMIC MODEL EXTRACTED FROM INHOMOGENEOUS MONTE-CARLO CALCULATIONS SHIN-CHI LEEr and TING-WEI TANG Department of Electrical and Computer Engineering, University of Massachusetts, Amherst, MA 01003, U.S.A.
(Received 12 March 1991; in revised form 4 September 1991) Abstract--A set of hydrodynamic equations is derived for a more accurate description of the hot-electron and nonstationary transport phenomena in Si devices containing local inhomogeneities in the submicron range. The transport coefficients are extracted from a Monte-Carlo calculation of a modeled test device which includes a rapidly varying electric field. The mobility is found to depend on the ratio of the average energy flux density to the current density (rather than the average energy) and the components of the electronic heat flux density are identified separately. This transport model can accurately predict velocity overshoot and carrier heating phenomena in realistic devices.
NOTATION
a,b
c~ Cw
C~p
D B° B~
D(E)
f h
L
k~ K
L(E)
#o n
b P
~
conv diff
R
Si
parameters of the mobility model of the electron gas electron random velocity ( = f , - I:) electron generation rate due to scatterings rate of change of the crystal momentum with time due to scattering rate of change of the kinetic energy with time due to scattering rate of change of quantity ~b with time due to scattering in the ith bin rate of change of 6 b with time due to scattering diffusivity tensor of the electron gas energy-dependent particle diffusivity energy-dependent energy diffusivity electric field-dependent particle diffusivity applied electric field quasi-electric field due to bandgap variation electron energy electron density distribution function in phase space total force exerted on the electron gas Planck constant divided by 27t current density of the electron gas Boltzmann constant thermal conductivity of the electron gas electric field-dependent characteristic length kinetic energy effective mass ( = 2~/v 2) mobility tensor of the electron gas energy flow mobility tensor of the electron gas low-field mobility directional component of the mobility energy-dependent mobility electron gas concentration crystal momentum ( = hk) average crystal momentum of the electron gas magnitude of the electron charge average heat flow vector of the electron gas convective part of O diffusive part of Q fourth order (gb/5) average tensor moment o f f average energy flow vector of the electron gas
tPresent address: $300, ITRI/ERSO, Hsinchu, Taiwan. 561
Sh t0, t~, te Ati 7",~ rp U U~ Ut
projection of S in the 7 direction average energy flow for homogeneous field condition parameters of the energy relaxation time model of the electron gas total time particle spent in the ith bin of the MC calculation average temperature tensor moment (m,b~/kB) off momentum relaxation time of the electron gas energy relaxation time of the electron gas second-order (t3,b) average tensor moment o f f longitudinal component of the diagonal terms of U transversal component of the diagonal terms of
g:
electron velocity high-field saturation velocity average velocity of the electron gas projection of V in the 7 direction average velocity for homogeneous field condition average kinetic energy of the electron gas average kinetic energy of the electron gas at thermal equilibrium average kinetic energy for homogeneous field condition parameter of the energy relaxation time model of the electron gas parameters of the average energy flux density model of the electron gas
v~
V V, Vh W W0 Wh
/~ ~',/~
I6f] ~
co,
rate of change o f f with time due to scattering
l. INTRODUCTION In m o d e r n - d a y submicron Si devices, high electric field ( I E r > ~ 1 0 0 k V c m -~) and high field-gradient (PEI/IVEI <~200 ]~) conditions are routinely encountered during operation and the carrier dynamics are far from thermal equilibrium. Hence, the driftdiffusion equations ( D D E ) model can be a p o o r a p p r o x i m a t i o n for these devices. The M o n t e - C a r l o
562
SH1N-CHILEE and TING-WEITANG
solution of Boltzmann's transport equation (BTE)[1], on the other hand, provides a more accurate description of carrier transport even under such highly nonequilibrium conditions. This is because various scattering mechanisms and band structure model are taken into consideration explicitly. However, the extensive computation required by the Monte-Carlo method makes it impractical for device design on a regular basis. Moreover, for regions of the device where diffusive transport is dominant, the result obtained from the Monte-Carlo method is often very noisy. The near-thermal equilibrium assumption of the DDE can be relaxed by including an additional energy balance equation[2]. These macroscopic transport equations are originally obtained from taking different moments of the BTE. This set of moment equations of the BTE, usually referred to as the hydrodynamic equations (HDE), have been proposed by Stratton[3] and Blotekjaer[4] to study the highfield transport properties in semiconductors. Different derivations have been proposed over the years[5] and the formulation has also been generalized to include more complicated semiconductor band structures[6] as well as heterostructure devices[7]. Many different attempts have been made to analyze devices such as n ~ n n ' diodes[8], IMPATT diodes[9], Si MOSFETs[10], Si BJTs[11], GaAs TEDs[12], MESFETs[13], HBTs[14] and MODFETs[15], using this model. As expected, this hydrodynamic equation approach can describe non-stationary phenomena such as velocity overshoot and carrier heating, which the DDE model fails to predict. However, the average velocity and energy calculated from this model do not always agree with that predicted by the MC method [16], even though the transport coefficients are obtained from the homogeneous Monte-Carlo calculations. In order to minimize this discrepancy, several difficult approaches have been proposed, For example, different values of energy relaxation time, ranging from 0.08 ps[17] to 0.5 ps[18] have been used. Recently, it has been reported that this disagreement can be reduced by including a nonparabolic formulation in the hydrodynamic model[19]. However, our study[20] shows that including the nonparabolicity of energy band structure alone cannot eliminate this discrepancy. It is the purpose of this study that through better extraction of the non-stationary transport information provided by the MC calculation, we can "calibrate" the macroscopic models and achieve a better agreement between the results predicted by the MC and HDE approaches. In Section 2, a particular set of moment equations is derived and the corresponding transport coefficients are chosen accordingly. A Monte-Carlo extraction method is described in Section 3, emphasizing the nonlocal effect due to a rapidly varying electric field. The calculated results and model construction are described in Section 4.
Some discussions and conclusions are given in Section 5. A new method developed to calculate the spatial variation of collision moments from the Monte-Carlo scheme is described in the Appendix. 2. H D E A N D T R A N S P O R T
COEFFICIENTS
The Boltzmann transport equation is given by:
dt Ot
df _ ~ +
b- V~f + ~" V-af = & coil'
(1)
wherefis the particle density distribution function in phase space (f¢, t) at a time t, b is the particle velocity and F(F) is the force exerted on the particle at location F which includes the external electric f i e l d - qE and the quasi-field induced by the bandgap variation. The collision term on the right-hand side of eqn (1) includes all the scattering mechanisms. The HDE usually consists of the first three moments of the BTE, representing the particle, momentum and energy conservations. However, there can be more than one form for these conservation equations, depending on the moments we choose. In this work, we select 1, the crystal momentum hf¢, and the particle kinetic energy o¢. Multiplying the BTE with these three quantities and integrating over the 1¢space yields[21]: On --
c~t
+ V "(nh
= nC.,
(n/~)
~--7- + v. (~rO - ~P= riG,
~(nW)
-
-
~
+ V ' (nS) -
-
n V. F = nC.,
(2)
(3) (4)
where n, V, P and W are the concentration, average velocity, average crystal momentum and average kinetic energy of the particles, respectively. They are defined as: r/
a__j ' f d£',
,a__! f hkf dr;
(5)
(7)
and
The tensor moment ~7, which represents the average flow of crystal momentum and the average energy flux density nS are defined as:
ngaf~fdf.
(10)
Silicon hydrodynamic model The terms on the right-hand side of eqns (2-4) are defined as: C . - ~ ! f[6~]¢o dk,
Cp--njI
/ ~ / //coil L
(11)
L ~ - q n ~ = -~7~ ' P + qb . (Vn)
+ ~ . (v. r_0, (17) is readily obtained from eqn (3) and a generalized Einstein relationship takes the form:
dL
(13,
C~ is zero here since no generation-recombination process is considered in this work. ~?p is the average rate of change of crystal momentum due to scattering, or the average friction force due to scattering, and C.. is the average rate of change of carrier kinetic energy due to scattering, or the average power absorption due to scattering. One well-known characteristic of the moment equations is that the ith order of the moment equation always contains a term of spatial gradient of the (i 4-1)th order moment and a term with the (i - 1)th order moment multiplied by the force. For example, in eqn (3), a second-order moment ~] appears on the left-hand side and eqn (4) includes a third-order moment S. Thus, it is necessary to express ~J and S in terms of P (or P) and W in order to close the hierarchy. Accurate models for these two terms are crucial for the accuracy of the HDE. Since the focus of this work is on steady-state behavior, all ~3/dt terms are subsequently neglected. Only an electron gas is discussed here, although the method can be applied to a hole gas as well. One popular expression for Cp is - P / r p , where zp is the momentum relaxation time. However, since it is k', not P, that appears in (2), it is more convenient to express Cp as - V / # , where mobility p is introduced, in order to couple the two eqns (2) and (3). For an isotropic parabolic band model, the difference between the two is simply a constant effective mass rn*. For a general band structure, however, m* is no longer a constant and the advantage of adopting the mobility expression is thus obvious: zp and m* can be lumped into p and only one transport coefficient p needs to be modeled. The band structure, or m*, does not appear in the HDE explicitly. With this in mind, we express Cp as:
~.p = _q#_+zS. . ~/,
(14)
where ~'--*~is the reciprocal electron mobility tensor, For a system with a homogeneous electric field, qJ~ = 6"p and eqn (14) becomes: ~ ' = --p • ~',
where ~?x is the transpose of ~', then a familiar phenomenological form of the electron current density:
(12)
and
c.~-~ ~lfw[~]co,
563
(15)
where ~ . ff~-~ equals identity tensor I. If the particle diffusivity tensor b is defined as:
y l . ~ = ! ~fT
(18) q For the expression of C~, we take the conventional energy relaxation time approach: C,-4
W-
W0,
(19)
~w where W o is the average kinetic energy of electrons in thermal equilibrium. Similar to #, zw is a functional of the distribution function and should not be treated as a constant. The energy flux density nS, defined in eqn (10), can be separated into two parts: the convective energy flux density and the heat flux density. If the electron energy ~ can be expressed as msv2/2, where rn~ = m~(~') for a general band structure, and the random velocity ? is defined as: --4~
-
~',
(20)
then the average heat flux density nO is defined as:
n(2~-f ~m~c~fdL
(21)
Equation (10) can then be rewritten as:
nS=nO+n(WT+k~r,)
• ['+-~
m, Of dk,
(22)
where the temperature tensor ~s is defined as S maO?fdk/nkB. The last term in eqn (22) is usually negligible even if m, is not a constant. Conventionally, the heat flux density nO is assumed to be proportional to the gradient of a scalar temperature: nO = - x VT,
(23)
where the thermal conductivity x is defined. Equation (23), however, cannot include the fact that (2 may not be zero even in a homogeneous system. Applying a uniform electric field in semiconductor, for example, always give rise to an asymmetric velocity distribution function and, according to eqn (21), yields a nonzero (2. Another approach to derive an expression for the energy flux density nS is to take one more higherorder moment equation of 8h# as follows[21,22]: f 8h~(BTE)dfc=~-V. (nR)
b - ~ l ~ - 8 T, q
(16/
--n(WI + ~J)" F=nfgp,
(24)
564
SmN-Cm LEE and T1NG-WEITANG
where
k~!fhe~gfd£
(25)
In analogy with eqn (14), we may express ~Tgp as: +~-i
C'gp = - qPs ' S.
(26)
Using eqns (3), (14), (24) and (26), the following expression for S can be derived:
,~=[~s. (w7+8) . . ~m] . p+l-
qn
~ . {(wT+ ~)
• [V" (nU)] - V. (n~)}.
(27)
Consequently, an expression for the average heat flux density nQ can be obtained from eqn (22) and (27) as: nO = nolo.,, + nO~fr,
(28)
where 0 ..... = [ ~ . (WI + U)- p~' - (WI + ku 7"a)] P (29) and 1 _
()~
- - p ~ - { ( w T + ~)" [ V ( n ~ ] - V" (nk)}. (30)
=
qn
Equation (28) generalizes the phenomenological Wiedemann-Franz law[23] which is actually an approximation for n0~i~ only, and not for n(), since the preassumption of the Wiedemann-Franz law is = 0. However, most devices are operated under a nonzero current flow condition and Q~on, must be included in the heat transport. In summary, to arrive at a set of close HDE equations, we must express the four terms Cp, C,., and ~ in eqn (3,4) in terms of the variables to be solved, i.e. n, P and W. The expression for (Tp is transformed to the problem of modeling the reciprocal mobility ~*-h and the modeling of C,, reduces to
!5
100 i~ , ,\ Q )
so !°--o E o
60 (
/
"o
~,~ ,~0 ~
P
~iI
~'~/°
L 20 !
~ 2
~
i
3
//f
........
~
! o
~
L oM 0 0
/ 0 2
/
C 4
/
//
) 6 }* ! ,j.r'-! ;
the modeling the energy relaxation time zw. It is worth mentioning that all the derivations in this section are independent of the material band structure. The modeling of #, %, U and S in terms of n, V and W for the I-D inhomogenous field case is described in the next section.
O8
' ,S
Fig. 1. (a) The electric field profile of test device (left axis) and the calculated electron concentration profile (right axis). (b) The corresponding potential profile.
3. M O N T E - C A R L O
EXTRACTION
The electron transport properties in Si, or any other semiconductor, are determined basically by the band structure and scattering mechanisms. One major problem of the HDE formulation is how to incorporate this microscopic information into the macroscopic transport coefficients, A Monte-Carlo code is developed in this work for that purpose. A nonparabolic band structure with a nonparabolicity factor of 0.5, acoustic phonon scattering, f,g-type intervalley phonon scattering and Conwell Weisskopf ionized impurity scattering are included. The expressions and parameters of the scattering rates are adopted from Ref. [24]. Besides statistical error, it has been proved[l] that the Monte-Carlo method solves the BTE and hence the solution can be used to extract the relationship between the transport coefficients and the moments of the distribution function. Conventionally, transport coefficients such as ~p, % are assumed to be functions of the average kinetic energy and the relationships are constructed by performing Monte-Carlo calculations under homogeneous field conditions[25]. These expressions are often extended to the case of an inhomogeneous field without any justification. In this work, an inhomogeneous electric field profile, as shown in Fig. la, is specified as a test device. Also shown in Fig. 1b is the corresponding potential profile. This device has a 0.2~m wide high-field (80kVcm 1) region sandwiched between two low-field (1.3 kV cm -]) regions. The two transition regions between uniform field regions are both chosen to be 0.1 pm, about 4-10 electron mean free paths. Significant nonstationary transport phenomena are expected to be present in these transition regions. The transport coefficients and the moments are calculated independently and their relationships under such conditions are examined. The impurity doping concentration is chosen to be 10~scm -3 everywhere for simplicity. This is obviously not a realistic device since the applied field and the calculated carrier concentration (shown in Fig. la) do not satisfy Poission's equation. This choice, however, is very useful for a general calibration purpose, since self-consistent solutions tend to smooth out the field gradient. The field gradient can be arbitrarily adjusted in this way and the impurity doping dependence of the nonstationary transport properties can be identified separately. Selfconsistent calculations, on the other hand, were performed for verification and similar conclusions were also obtained[26].
Silicon hydrodynamic model
--~ o
565
1.4
06
10
V
1.2
w
0
0 4 r"o
~ 1
%" 7.0 / I ~ c ~ -~ 08 /
T
,
0,8
06
~~o_
q
05 ca
<--
I 0.6
2>
t
3, ,t-
02
<]
I
.~
/3
04-
CP I
04L
0.2
0 C) £ 0
0.4
02
0.8
0.6 x (~)
! .0
On
O0
Fig. 2. The average velocity V (left axis), average kinetic energy W (right axis) and average energy flow S (right axis) from the Monte-Carlo calculation of the test device. The velocity anticipatory and overshoot phenomena are significant. 4. C A L C U L A T I O N
100
I
E 0.eL
I
/
E o
i
6o
.~.
40
o
2O -> 0.0 ~ oo
0.2
o.~
o.6 x (ff~)
~ " ~ "- "~ ~ -
0.8
04
06
0.8
1.0
Fig. 4. The energy relaxation time % (left axis) and the average power dissipation -C~ (right axis) of the test device. Energy relaxation time is calculated from % = - ( W - Wo)/Cw and the Monte-Carlo data are noisy only in the low-field regions.
4.1, Mobility and energy relaxation time
RESULTS
A well-known nonlocal phenomenon is described by the average velocity profile in Fig. 2. Not only does the velocity increase before the electric field starts to ramp up (anticipatory effect[27]), but an obvious peak (overshoot) and a less significant dip (undershoot) are also observed. Also shown in Fig. 2 are the average energy and the average energy flow. No overshoot phenomena for these two moments are observed. It is worth mentioning that although n and V are not uniform individually, the product nV is constant throughout the region, since no generation or recombination of electrons is considered.
]>
02
I .,0
0
Fig. 3. The mobility # (left axis) and the average frictional force Cp (right axis) of the test device. Mobility profile is calculated from expression p = - V/Cp and the MonteCarlo data are noisy only in the low field regions.
The two collision moments Cp and Cw are calculated directly by the Monte-Carlo method. Details are discussed in the Appendix. Using Cp, C,. and the calculated V, W, the spatial variation o f p and % can be obtained. Figure 3 shows the profiles of Cp and p and Fig. 4 shows the calculated values of C,, and %. Although C, and Cw follow the electric field shape in general, a spatial delay is observed, especially in the region where the electric field increases rapidly. This clearly demonstrates the nonlocal transport property of the electron gas where the response lags behind the driving field. It is more complicated to examine the profiles of/~ and ~,. Although the data are noisy in the low-field regions, excellent results are obtained in the high-field and field-transition regions. These data are used subsequently to examine the possible models for p and ~,. Three possible mobility models are compared in Fig. 5. Figure 5a shows the reciprocal mobility p - ~vs electric field. The hysteresis loop clearly indicates that p is not a single-valued function of the local electric field, Compared to homogeneous field data, kt is greater if dJE[/dx is positive and smaller if d[E[/dx is negative. If/~ -~ is plotted against W, as shown in Fig. 5b, there is still a hysteresis loop. Although the homogeneous-field data approximate the data in the positive dlEI/dx region, they still overestimate the value of/~ in the region with a negative dlEWdx. This implies that the nonstationary transport phenomena cannot be properly described by any p =/~(W) model. If, however, the moment expansion suggested in Ref. [28] is followed, i.e. to express C'p as:
Cp = - q ( a V + bS),
(31)
where a, b are the expansion coefficients, then:
p -~ = a + b(S/V).
(32)
566
S H I N - C H I LEE a n d
'~
•
homog
c,'
inhomog
o
a
6
°
•
~
ij
"
oe
2">
(C) "
20
r}
4tt~
1+~-2 E 2
1/2 E2
(36)
and the impurity-concentration dependence can be included conveniently by adopting any empirical form of P0 (NI). A similar approach can be applied to the modeling of the energy relaxation time. If %. is plotted against W, as shown in Fig. 6, only one path is followed for all data. Homogeneous field calculations also show a similar relationship. Hence, an empirical expression r , = z,,(W) can be applied to fit all data. A fourparameters empirical form:
J
o
_o o
1+
o o
o
W=W°+q%
t
o
7~
TING-WEI TANG
4C, 60 - £ (kV,/cm)
80
! '}
0
8
o9!
.....
q'~omog
o o
o0
!
• o
I
o a _J o)
T~ 2 t
i
o
L,(W)=to+fi(W/Wo-l)
o
o,O
i
,. I' U
+t=exp(-fl(W/Wo-1)),
i
O0
01
0 2
is proposed here, where to, tl, t2 and/~ are chosen to be 0.28 ps, 3 fs, 2.2 ps and 10, respectively. This result basically supports the conventional assumption that % is a single-valued function of the average energy, even in the presence of a rapidly changing field. This % - W relationship is very sensitive to the band structure model assumed in the calculation and a parabolic band model would be an oversimplification [29].
0
(ev) F
6 k
~3
>
fs
o - -
]
=a¢o(S.,'\.,',
~
4
J
O~ O0
O1
02
O}
4.2. Energy tensor and heat flow
04
s / v (ev) Fig. 5. The reciprocal mobility is plotted against: (a) electric field; (b) average kinetic energy W; and (c) ratio between average energy flow S and average velocity V. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous calculation. The arrows in (a), (b) indicate the directions of the hysteresis loops and the solid line in (c) is a linear fit of the data.
Figure 5c demonstrates that this is a much more accurate model of #, since all data follow this linear relationship independent of field variations. The parameters a and b can be determined directly from Fig. 5c, or they can be expressed by the more familiar parameters: low-field mobility P0 and saturation velocity %, as: a
1
w0
(33)
l~o q~:,,v~'
In general, the energy tensor U has six components; however, for the I-D case, ~/is a diagonal matrix and only the longitudinal component Ut and the transverse component Ut need to be considered. Figure 7 shows that UI is greater than Ut, especially in the high-field region, indicating that the distribution function is elongated along the field direction and the
b = ~ ,
(34)
•
homog
I
o
inhomoq,
i
---G? D_ .J
fittir, g
£;6
'2 I
[
1
r
'., ~
1
(37)
v _ r ~ _ 4~.~q~_~_ ~
02
i t
where S is a parameter to be discussed in the next section. This transformation gives the electric field dependence o f the average velocity and energy of bulk Si as: =
- 2/~0 4,2
1+
l+~-E2)
\1.2
~,
(35)
00 00
0 '
02
I 0.3
w (ev) Fig. 6. The energy relaxation time plotted against the average kinetic energy. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous field calculation. The solid line is an empirical fit of the data.
Silicon hydrodynamic model
• 02
U
o
Ut
--
2w/5
.7 o 0~
~J
,l ° °>
<~ !
,,/o o/o Ic,
i
% [ O0
i
I O0
O2
O~
0 6
08
567
observed in a uniform-field condition, where Q is always in the opposite direction of S. This is a direct manifestation of the change in the shape of the distribution function. The contribution of Q to S is significant in the high-field region, mainly due to Q .... term, and neglecting Q would give an almost 30% overestimation of S in the region. The closure requirement of the H D E still needs an expression of the energy flux density nS by the lower-order moments. If we can assume that the fourth-order moment R has an quadratic dependence on W and ~q is proportional to #[20], then Qdiff may be approximated by a form similar to the Wiedem a n n - F r a n z law: 7~
10
< ,:,~m)
Qdiff =
Fig. 7. The longitudinal component Ul and the transverse component Ut of the tensor U. The solid line is two-thirds of the average energy W which approximates U~fairly well.
equipartition assumption of the energy is invalid. However, if 2 W/3 is also plotted in Fig. 7, we can see that this quantity approximates G fairly well. This is because the effective mass m, used in the definition of Wis greater than the effective mass used in ~, which effectively simulates the effect of the elongated distribution function on U~. Other empirical models of G are also possible[20]. No attempt is made to model U,, however, since it does not affect the electron transport for 1-D problems. The total Q in the device is plotted in Fig. 8, in which the convective part Qconvand the diffusive part Qaifr are identified seperately. Q .... is the dominant term in the uniform field regions where all the spatial gradients are negligible while Qdierdominates over the field ramp regions. It is interesting to see that Q is not always negative. Near x = 0.7 #m, Q is actually in the same direction as S. No such phenomenon is
(38)
dx
S = ( W + k B Te, ) V + Q ..... +
Qdiff,
(39)
,~ g V W "}- Qdiff,
(40) dW
~-SVW-
~W--
q
(41)
dx
Figure 9 shows S plotted against V W taken from the Monte-Carlo data for both homogeneous and in-
!i '3'.}
F ?
•
homcg
o
nbomoc
/~o,@/~" Y. ,61" /
%-
O0
q
The parameter A differs from the conventionally defined Lorenz number by a factor of (3/2) 2, since an energy formulation, rather than a temperature formulation[3], is adopted here. The homogeneous field calculations suggest that ( 2 ~ - ( 1 / 3 ) W V [ 2 0 ] and hence, a simple expression of Q ..... proportional to VW, is assumed. If the convection energy flow ( W + k B T~)V is also assumed to be proportional to VW, then it can be combined with the Q .... term and the energy flow S may now be expressed as:
01
52 E o
dW
-- --/Z W - - .
. . . . }'vw
/
,
,xs~t.." ,
02 /
o
• ,/
/ Dq
o
t/
© {
] ¢I
/
SI - 0®~
Q ~ ;
I
....
0.0 0.0
01
0.2
03
VW , ' 1 0 7 e V * c m / s )
I
02 ]
O0
02
0.4
0.6 x (~)
0.8
10
Fig. 8. The total average heat flow Q ( ) and its two components: convective heat flow Q.... ( 0 ) and diffusive heat flow Qai~r(0).
Fig. 9. The energy flow S plotted against the product of average velocity V and the average kinetic energy W. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous field calculation. Dashed line is 6 VW, calculated from the homogeneous field data, and the solid line is o~VW- z~#W V W/q, calculated from the inhomogeneous field data. ~" and z~ are 1.29 and 0.2, respectively.
568
SHIN-CHI LEE and TING-WEI TANG
homogeneous-field conditions. For the uniform field case, S is predicted fairly well by selecting ~ = 1.29. This 3 V W term combines the convection energy flow and the convective part of the heat flow. For the nonuniform-field case, eqn (41) qualitatively predicts the characteristics of S with a A of 0.2. These values of the parameters S and ~ correspond to values of 1.95 and 0.45 for the parameters 6 and A defined in Ref. [3], respectively. However, there are still slight deviations of the model from the Monte-Carlo calculated S, indicating that the Wiedemann-Franz law is only a first-order approximation for Qain. Other empirical models of Qd~ are also possible [20]. This set of closed HDE with transport coefficients given by (32), (37) and (41) is applied to simulate real devices with significant nonstationary transport phenomena. Good agreement is obtained between the Monte-Carlo calculated data and the HDE solutions [31]. 5. DISCUSSIONS AND CONCLUSIONS
The derivation of the HDE and associated transport coefficients is given in this study. In order to close the moment equation hierarchy, emphasizing the validity in devices operated under rapidly changing electric field conditions, an inhomogeneous field Monte-Carlo extraction scheme is developed to examine the appropriateness of the transport coefficient models. The final HDE model proposed for device simulations can be summarized as follows:
different definitions of the effective mass. The advantage of this approach is that the form of the equations is simpler and it is more adaptable to any generalization to a real band structure, where no analytical form of d'(Yc) can be found. All the band structure information is included in the empirical expressions of the transport coefficients. It is worthwhile to discuss the relationships between this model and other HDE formulations. If expressions (31), (46) are adopted, instead of using mobility, then (43) takes the form: 2
V(n W) - nf" = - qn (a ~>+ bS) = aJ - qbnS.
Combined with (47), we can obtain: _2V(n W) - nF = (a + b6 W)J + qnpb A W A W. (51) 3 This suggests that an equivalent current density equation can be expressed as[32]: ~)= - n/iF + q/9,, Vn + n/),~ V W,
V
(42)
/~V(nU) - n p F = ),,,
(43)
1_ q
(nS) + - J,,. P =
W - Wo r,
- n -
(44)
where
P= -q(;~+ ;~), 2 nS=gi=
q
(45)
U = - W. 3
(46)
W)~--npWVW, q
(47)
b qnSi Ia--(~-,,,)J
1
,
/i ~ (a + bgW) 1
(53)
2W ~,, ~/~ - - 3q'
(54)
p \ ~ - qb Al~ W).
(55)
~zx'2(.
Equations (43), (52) are equivalent only if the transport coefficients are defined according to the relationships described above. Another form of the I-D current equation that has often been referred to as the augmented drift diffusion current equation[33] can also be approximated from (50). If the field-dependent average velocity, average energy and average energy flow for a homogeneous system are designated as Vh, W h and Sh, respectively, then: - E = aV h + bSh. (56)
(48)
If the nonstationary phenomena of the energy flow are neglected and S is approximated by Sh in (50) then: 2 d V = Vh(E ) (nW). (57) 3qan dx
(49)
If W is further approximated by Wh, (57) becomes:
r,, = to+ t t ( W / W o - 1) + t2exp(-fl(W/Wo -1)).
(52)
if an equivalent energy-dependent mobility /i, an equivalent energy-dependent particle diffusivity 15. and an equivalent energy diffusivity/), are defined as:
D V . ]1,,= qn( R - G),
(50)
Values of parameters 3", A, to, tl, t: and fl are selected as 1.29, 0.2, 0.28 ps, 3 fs, 2.2 ps and 10, respectively. It should be pointed out that models of U, S and r,. are empirical approximations only. Further improvements are possible. One of the considerations for choosing the form derived here is to avoid any explicit average mass terms appearing in the equations. This reduces the number of transport coefficients needed in the equations to a minimum, since there are several
V = Vh(E)[1
L(E) dE dx
D(E)n dx'an
(58)
where the field-dependent diffusivity D(E) and the characteristic length L ( E ) are defined as: D(E)~_ 2Wh(E) ' 3qa L(E) ~
2E dWh 3qaVh(E ) dF
Silicon hydrodynamic model The a d v a n t a g e of this f o r m u l a t i o n is t h a t energy i n f o r m a t i o n is no longer needed for ocity calculation. However, the a s s u m p t i o n p e r t u r b a t i o n from the h o m o g e n e o u s state valid for the weak field gradient regions applicable range o f this expression is very
average the velof small is only a n d the limited.
Acknowledgements--This work was supported in part by
National Science Foundation Grant No. ECS- 8611456. We have also benefited from the General Technology Division, East Fishkill Facility, of the IBM corporation. The authors also wish to thank Professor David Navon for carefully reading the manuscript. REFERENCES
1. W. Fawcett, A. D. Boardman and S. Swain, J. Phys. Chem. Solids 31, 1963 (1970). 2. B. Meinerzhagen and W. L. Engl, IEEE Trans. Electron Devices ED-35, 689 (1988). 3. R. Stratton, Phys. Rev. 126, 2002 (1962). 4. K. Blgtekjaer, IEEE Trans. Electron Devices ED-17, 38 (1970). 5. M. Sanchez, Solid-St. Electron. 16, 549 (1973); W. Htinsch, In Proc. Two-Dimensional Systems: Physics and Devices, p. 296. Springer, Berlin (1986); E. M. Azoff, Solid-St. Electron. 30, 913 (1987); C. C. McAndrew, E. L. Heasell and K. Singhal, Semicond. Sei. Technol. 2, 643 (1987). 6. D. L. Woolard, R. J. Trew and M. A. Littlejohn, Solid-St. Electron. 31, 571 (1988); R. Thoma, A. Emunds, B. Meinerzhagen, H. Peifer and W. L. Engl, I E D M Tech. Dig., p. 139 (1989). 7. C, C. McAndrew, E. L. Heasell and K. Singhal, Semicond. Sei. Technol. 2, 643 (1987); E. M. Azoff, J. appl. Phys. 64, 2439 (1988). 8. M. Rudan, F. Odeh and J. White, C O M P E L 6, 151 (1987). 9. R, K. Froelich, Ph.D. Dissertation, University of Michigan, Ann Arbor (1982). 10. M. Fukuma and R. H. Uebbing, I E D M Tech. Dig., p. 621 (1984); M. Tomizawa, K. Yokoyama and A. Yoshii, IEEE Trans. Comput. Aided-Des. CAD-7, 254 (1988). 1I. R. K. Cook, IEEE Trans. Electron Devices ED-30, 1103 (1983); G. Baccarani and M. R. Wordeman, Solid-St. Electron. 28, 407 (1985). 12. R. Bosch and H. W. Thim, 1EEE Trans. Electron Devices ED-21, 16 (1974). 13. W. R. Curtice and Y. H. Yun, IEEE Trans. Electron Devices ED-28, 954 (1981); R. K. Cook and J. Frey, IEEE Trans. Electron P - ~ ~,9, 970 (1982); C. M. Snowden and D ,s. Electron Devices ED-34, 21 ~ 14. E. ~"
evices ED-36, 609
t, IEEE Trans. "ctron Devices
koyama and TAD-7, 254 ai, IEEE IEEE
569
19. T. J. Bordelon, X. L. Wang, C. M. Maziar and A. F. Tasch, I E D M Tech. Dig., p. 353 (1990). 20. S.-C. Lee, Ph.D. Dissertation, University of Massachusetts, Amherst (1990). 21. A. Bringer and G. Sch6n, J. appl. Phys. 64, 2447 (1988). 22. R. Thoma, A. Emunds, B. Meinerzhagen, H. J. Peifer and W. L. Engl, 1EEE Trans. Electron Devices ED-38, 1343 (1991). 23. N. W. Ashcroft and N. D. Mermin, Solid State Physics. Saunders College, Philadelphia (1976). 24. C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983). 25. M. Shur, Electron Lett. 12, 615 (1976); R. Thoma, A. Emunds, B. Meinerzhagen, H. J. Peifer and W. L. Engl, IEEE Trans. Electron Devices ED-38, 1343 (1991). 26. Y. J. Park, Private Communication. To be published. 27. S. Bandyopadhyay, M. E. Brown, C. M. Maziar, S. Datta and M. S. Lundstrom, IEEE Trans. Electron Devices ED-34, 392 (1987). 28. W. H~insch and M. Miura-Mattausch, J. appl. Phys. 60, 650 (1986). 29. S.-C. Lee and T.-W. Tang, Computational Electronics, p. 127, (Edited by K. Hess, J. P. Leburton and U. Ravaioli). Kluwer Academic, Boston (1991). 30. P. J. Price, I B M J. Res. Dev. 14, 12 (1970). 31. S.-C. Lee and T.-W. Tang, Proc. 1991 VPAD, p. 47 (1991). 32. R. K. Cook, IEEE Trans. Electron Devices ED-30, 1103 (1983). 33. K. K. Thornber, IEEE Electron Device Lett. EDL-3, 69 (1982).
APPENDIX
Collision Moment Calculation
In Monte-Carlo calculations, each spatial bin is treated as one slab of bulk material with a constant electric field. There are two fundamental processes taking place in each bin: free flight and scattering. Noticing that the collision integral on the right-hand side of the BTE is basically the difference between the A-state distribution and the B-state distribution[30], a direct approach is developed in the following to calculate the average collision moments of each bin. The ¢(fQ-collision moment in ith bin, C~j, is defined as the rate of change of function ¢ due to collisions taking place in this bin. If the states of the particle before and after the j t h scattering event are indicated by f% and ~aj, respectively, then the change of ~ due to this j t h scattering is simply ¢(k~j) - ¢(/%). Hence, the value of C~ at this ith bin can be calculated by:
C,~,,-
Ati
(61)
where the summation runs over all the scattering events taking place in the ith bin and A6 is the total flight time the particle spends in this bin. If function ¢ is substituted by hfc or g, 6"p or C~ of each bin can be calculated directly, and the spatial variation of ~p and C, are obtained. Since Monte-Carlo results automatically satisfy all the moment equations, the spatialgradient terms of the moment equations can also be calculated from the differences between the collision terms and the field terms. This approach avoids the use of the conventional finite-difference scheme on the noisy Monte-Carlo data to calculate the spatial-gradients and a much smoother function over the coordinate space can be produced.
0038-1101/92 $5.00+ 0.00 Copyright © 1992 Pergamon Press plc
T R A N S P O R T COEFFICIENTS FOR A SILICON HYDRODYNAMIC MODEL EXTRACTED FROM INHOMOGENEOUS MONTE-CARLO CALCULATIONS SHIN-CHI LEEr and TING-WEI TANG Department of Electrical and Computer Engineering, University of Massachusetts, Amherst, MA 01003, U.S.A.
(Received 12 March 1991; in revised form 4 September 1991) Abstract--A set of hydrodynamic equations is derived for a more accurate description of the hot-electron and nonstationary transport phenomena in Si devices containing local inhomogeneities in the submicron range. The transport coefficients are extracted from a Monte-Carlo calculation of a modeled test device which includes a rapidly varying electric field. The mobility is found to depend on the ratio of the average energy flux density to the current density (rather than the average energy) and the components of the electronic heat flux density are identified separately. This transport model can accurately predict velocity overshoot and carrier heating phenomena in realistic devices.
NOTATION
a,b
c~ Cw
C~p
D B° B~
D(E)
f h
L
k~ K
L(E)
#o n
b P
~
conv diff
R
Si
parameters of the mobility model of the electron gas electron random velocity ( = f , - I:) electron generation rate due to scatterings rate of change of the crystal momentum with time due to scattering rate of change of the kinetic energy with time due to scattering rate of change of quantity ~b with time due to scattering in the ith bin rate of change of 6 b with time due to scattering diffusivity tensor of the electron gas energy-dependent particle diffusivity energy-dependent energy diffusivity electric field-dependent particle diffusivity applied electric field quasi-electric field due to bandgap variation electron energy electron density distribution function in phase space total force exerted on the electron gas Planck constant divided by 27t current density of the electron gas Boltzmann constant thermal conductivity of the electron gas electric field-dependent characteristic length kinetic energy effective mass ( = 2~/v 2) mobility tensor of the electron gas energy flow mobility tensor of the electron gas low-field mobility directional component of the mobility energy-dependent mobility electron gas concentration crystal momentum ( = hk) average crystal momentum of the electron gas magnitude of the electron charge average heat flow vector of the electron gas convective part of O diffusive part of Q fourth order (gb/5) average tensor moment o f f average energy flow vector of the electron gas
tPresent address: $300, ITRI/ERSO, Hsinchu, Taiwan. 561
Sh t0, t~, te Ati 7",~ rp U U~ Ut
projection of S in the 7 direction average energy flow for homogeneous field condition parameters of the energy relaxation time model of the electron gas total time particle spent in the ith bin of the MC calculation average temperature tensor moment (m,b~/kB) off momentum relaxation time of the electron gas energy relaxation time of the electron gas second-order (t3,b) average tensor moment o f f longitudinal component of the diagonal terms of U transversal component of the diagonal terms of
g:
electron velocity high-field saturation velocity average velocity of the electron gas projection of V in the 7 direction average velocity for homogeneous field condition average kinetic energy of the electron gas average kinetic energy of the electron gas at thermal equilibrium average kinetic energy for homogeneous field condition parameter of the energy relaxation time model of the electron gas parameters of the average energy flux density model of the electron gas
v~
V V, Vh W W0 Wh
/~ ~',/~
I6f] ~
co,
rate of change o f f with time due to scattering
l. INTRODUCTION In m o d e r n - d a y submicron Si devices, high electric field ( I E r > ~ 1 0 0 k V c m -~) and high field-gradient (PEI/IVEI <~200 ]~) conditions are routinely encountered during operation and the carrier dynamics are far from thermal equilibrium. Hence, the driftdiffusion equations ( D D E ) model can be a p o o r a p p r o x i m a t i o n for these devices. The M o n t e - C a r l o
562
SH1N-CHILEE and TING-WEITANG
solution of Boltzmann's transport equation (BTE)[1], on the other hand, provides a more accurate description of carrier transport even under such highly nonequilibrium conditions. This is because various scattering mechanisms and band structure model are taken into consideration explicitly. However, the extensive computation required by the Monte-Carlo method makes it impractical for device design on a regular basis. Moreover, for regions of the device where diffusive transport is dominant, the result obtained from the Monte-Carlo method is often very noisy. The near-thermal equilibrium assumption of the DDE can be relaxed by including an additional energy balance equation[2]. These macroscopic transport equations are originally obtained from taking different moments of the BTE. This set of moment equations of the BTE, usually referred to as the hydrodynamic equations (HDE), have been proposed by Stratton[3] and Blotekjaer[4] to study the highfield transport properties in semiconductors. Different derivations have been proposed over the years[5] and the formulation has also been generalized to include more complicated semiconductor band structures[6] as well as heterostructure devices[7]. Many different attempts have been made to analyze devices such as n ~ n n ' diodes[8], IMPATT diodes[9], Si MOSFETs[10], Si BJTs[11], GaAs TEDs[12], MESFETs[13], HBTs[14] and MODFETs[15], using this model. As expected, this hydrodynamic equation approach can describe non-stationary phenomena such as velocity overshoot and carrier heating, which the DDE model fails to predict. However, the average velocity and energy calculated from this model do not always agree with that predicted by the MC method [16], even though the transport coefficients are obtained from the homogeneous Monte-Carlo calculations. In order to minimize this discrepancy, several difficult approaches have been proposed, For example, different values of energy relaxation time, ranging from 0.08 ps[17] to 0.5 ps[18] have been used. Recently, it has been reported that this disagreement can be reduced by including a nonparabolic formulation in the hydrodynamic model[19]. However, our study[20] shows that including the nonparabolicity of energy band structure alone cannot eliminate this discrepancy. It is the purpose of this study that through better extraction of the non-stationary transport information provided by the MC calculation, we can "calibrate" the macroscopic models and achieve a better agreement between the results predicted by the MC and HDE approaches. In Section 2, a particular set of moment equations is derived and the corresponding transport coefficients are chosen accordingly. A Monte-Carlo extraction method is described in Section 3, emphasizing the nonlocal effect due to a rapidly varying electric field. The calculated results and model construction are described in Section 4.
Some discussions and conclusions are given in Section 5. A new method developed to calculate the spatial variation of collision moments from the Monte-Carlo scheme is described in the Appendix. 2. H D E A N D T R A N S P O R T
COEFFICIENTS
The Boltzmann transport equation is given by:
dt Ot
df _ ~ +
b- V~f + ~" V-af = & coil'
(1)
wherefis the particle density distribution function in phase space (f¢, t) at a time t, b is the particle velocity and F(F) is the force exerted on the particle at location F which includes the external electric f i e l d - qE and the quasi-field induced by the bandgap variation. The collision term on the right-hand side of eqn (1) includes all the scattering mechanisms. The HDE usually consists of the first three moments of the BTE, representing the particle, momentum and energy conservations. However, there can be more than one form for these conservation equations, depending on the moments we choose. In this work, we select 1, the crystal momentum hf¢, and the particle kinetic energy o¢. Multiplying the BTE with these three quantities and integrating over the 1¢space yields[21]: On --
c~t
+ V "(nh
= nC.,
(n/~)
~--7- + v. (~rO - ~P= riG,
~(nW)
-
-
~
+ V ' (nS) -
-
n V. F = nC.,
(2)
(3) (4)
where n, V, P and W are the concentration, average velocity, average crystal momentum and average kinetic energy of the particles, respectively. They are defined as: r/
a__j ' f d£',
,a__! f hkf dr;
(5)
(7)
and
The tensor moment ~7, which represents the average flow of crystal momentum and the average energy flux density nS are defined as:
ngaf~fdf.
(10)
Silicon hydrodynamic model The terms on the right-hand side of eqns (2-4) are defined as: C . - ~ ! f[6~]¢o dk,
Cp--njI
/ ~ / //coil L
(11)
L ~ - q n ~ = -~7~ ' P + qb . (Vn)
+ ~ . (v. r_0, (17) is readily obtained from eqn (3) and a generalized Einstein relationship takes the form:
dL
(13,
C~ is zero here since no generation-recombination process is considered in this work. ~?p is the average rate of change of crystal momentum due to scattering, or the average friction force due to scattering, and C.. is the average rate of change of carrier kinetic energy due to scattering, or the average power absorption due to scattering. One well-known characteristic of the moment equations is that the ith order of the moment equation always contains a term of spatial gradient of the (i 4-1)th order moment and a term with the (i - 1)th order moment multiplied by the force. For example, in eqn (3), a second-order moment ~] appears on the left-hand side and eqn (4) includes a third-order moment S. Thus, it is necessary to express ~J and S in terms of P (or P) and W in order to close the hierarchy. Accurate models for these two terms are crucial for the accuracy of the HDE. Since the focus of this work is on steady-state behavior, all ~3/dt terms are subsequently neglected. Only an electron gas is discussed here, although the method can be applied to a hole gas as well. One popular expression for Cp is - P / r p , where zp is the momentum relaxation time. However, since it is k', not P, that appears in (2), it is more convenient to express Cp as - V / # , where mobility p is introduced, in order to couple the two eqns (2) and (3). For an isotropic parabolic band model, the difference between the two is simply a constant effective mass rn*. For a general band structure, however, m* is no longer a constant and the advantage of adopting the mobility expression is thus obvious: zp and m* can be lumped into p and only one transport coefficient p needs to be modeled. The band structure, or m*, does not appear in the HDE explicitly. With this in mind, we express Cp as:
~.p = _q#_+zS. . ~/,
(14)
where ~'--*~is the reciprocal electron mobility tensor, For a system with a homogeneous electric field, qJ~ = 6"p and eqn (14) becomes: ~ ' = --p • ~',
where ~?x is the transpose of ~', then a familiar phenomenological form of the electron current density:
(12)
and
c.~-~ ~lfw[~]co,
563
(15)
where ~ . ff~-~ equals identity tensor I. If the particle diffusivity tensor b is defined as:
y l . ~ = ! ~fT
(18) q For the expression of C~, we take the conventional energy relaxation time approach: C,-4
W-
W0,
(19)
~w where W o is the average kinetic energy of electrons in thermal equilibrium. Similar to #, zw is a functional of the distribution function and should not be treated as a constant. The energy flux density nS, defined in eqn (10), can be separated into two parts: the convective energy flux density and the heat flux density. If the electron energy ~ can be expressed as msv2/2, where rn~ = m~(~') for a general band structure, and the random velocity ? is defined as: --4~
-
~',
(20)
then the average heat flux density nO is defined as:
n(2~-f ~m~c~fdL
(21)
Equation (10) can then be rewritten as:
nS=nO+n(WT+k~r,)
• ['+-~
m, Of dk,
(22)
where the temperature tensor ~s is defined as S maO?fdk/nkB. The last term in eqn (22) is usually negligible even if m, is not a constant. Conventionally, the heat flux density nO is assumed to be proportional to the gradient of a scalar temperature: nO = - x VT,
(23)
where the thermal conductivity x is defined. Equation (23), however, cannot include the fact that (2 may not be zero even in a homogeneous system. Applying a uniform electric field in semiconductor, for example, always give rise to an asymmetric velocity distribution function and, according to eqn (21), yields a nonzero (2. Another approach to derive an expression for the energy flux density nS is to take one more higherorder moment equation of 8h# as follows[21,22]: f 8h~(BTE)dfc=~-V. (nR)
b - ~ l ~ - 8 T, q
(16/
--n(WI + ~J)" F=nfgp,
(24)
564
SmN-Cm LEE and T1NG-WEITANG
where
k~!fhe~gfd£
(25)
In analogy with eqn (14), we may express ~Tgp as: +~-i
C'gp = - qPs ' S.
(26)
Using eqns (3), (14), (24) and (26), the following expression for S can be derived:
,~=[~s. (w7+8) . . ~m] . p+l-
qn
~ . {(wT+ ~)
• [V" (nU)] - V. (n~)}.
(27)
Consequently, an expression for the average heat flux density nQ can be obtained from eqn (22) and (27) as: nO = nolo.,, + nO~fr,
(28)
where 0 ..... = [ ~ . (WI + U)- p~' - (WI + ku 7"a)] P (29) and 1 _
()~
- - p ~ - { ( w T + ~)" [ V ( n ~ ] - V" (nk)}. (30)
=
qn
Equation (28) generalizes the phenomenological Wiedemann-Franz law[23] which is actually an approximation for n0~i~ only, and not for n(), since the preassumption of the Wiedemann-Franz law is = 0. However, most devices are operated under a nonzero current flow condition and Q~on, must be included in the heat transport. In summary, to arrive at a set of close HDE equations, we must express the four terms Cp, C,., and ~ in eqn (3,4) in terms of the variables to be solved, i.e. n, P and W. The expression for (Tp is transformed to the problem of modeling the reciprocal mobility ~*-h and the modeling of C,, reduces to
!5
100 i~ , ,\ Q )
so !°--o E o
60 (
/
"o
~,~ ,~0 ~
P
~iI
~'~/°
L 20 !
~ 2
~
i
3
//f
........
~
! o
~
L oM 0 0
/ 0 2
/
C 4
/
//
) 6 }* ! ,j.r'-! ;
the modeling the energy relaxation time zw. It is worth mentioning that all the derivations in this section are independent of the material band structure. The modeling of #, %, U and S in terms of n, V and W for the I-D inhomogenous field case is described in the next section.
O8
' ,S
Fig. 1. (a) The electric field profile of test device (left axis) and the calculated electron concentration profile (right axis). (b) The corresponding potential profile.
3. M O N T E - C A R L O
EXTRACTION
The electron transport properties in Si, or any other semiconductor, are determined basically by the band structure and scattering mechanisms. One major problem of the HDE formulation is how to incorporate this microscopic information into the macroscopic transport coefficients, A Monte-Carlo code is developed in this work for that purpose. A nonparabolic band structure with a nonparabolicity factor of 0.5, acoustic phonon scattering, f,g-type intervalley phonon scattering and Conwell Weisskopf ionized impurity scattering are included. The expressions and parameters of the scattering rates are adopted from Ref. [24]. Besides statistical error, it has been proved[l] that the Monte-Carlo method solves the BTE and hence the solution can be used to extract the relationship between the transport coefficients and the moments of the distribution function. Conventionally, transport coefficients such as ~p, % are assumed to be functions of the average kinetic energy and the relationships are constructed by performing Monte-Carlo calculations under homogeneous field conditions[25]. These expressions are often extended to the case of an inhomogeneous field without any justification. In this work, an inhomogeneous electric field profile, as shown in Fig. la, is specified as a test device. Also shown in Fig. 1b is the corresponding potential profile. This device has a 0.2~m wide high-field (80kVcm 1) region sandwiched between two low-field (1.3 kV cm -]) regions. The two transition regions between uniform field regions are both chosen to be 0.1 pm, about 4-10 electron mean free paths. Significant nonstationary transport phenomena are expected to be present in these transition regions. The transport coefficients and the moments are calculated independently and their relationships under such conditions are examined. The impurity doping concentration is chosen to be 10~scm -3 everywhere for simplicity. This is obviously not a realistic device since the applied field and the calculated carrier concentration (shown in Fig. la) do not satisfy Poission's equation. This choice, however, is very useful for a general calibration purpose, since self-consistent solutions tend to smooth out the field gradient. The field gradient can be arbitrarily adjusted in this way and the impurity doping dependence of the nonstationary transport properties can be identified separately. Selfconsistent calculations, on the other hand, were performed for verification and similar conclusions were also obtained[26].
Silicon hydrodynamic model
--~ o
565
1.4
06
10
V
1.2
w
0
0 4 r"o
~ 1
%" 7.0 / I ~ c ~ -~ 08 /
T
,
0,8
06
~~o_
q
05 ca
<--
I 0.6
2>
t
3, ,t-
02
<]
I
.~
/3
04-
CP I
04L
0.2
0 C) £ 0
0.4
02
0.8
0.6 x (~)
! .0
On
O0
Fig. 2. The average velocity V (left axis), average kinetic energy W (right axis) and average energy flow S (right axis) from the Monte-Carlo calculation of the test device. The velocity anticipatory and overshoot phenomena are significant. 4. C A L C U L A T I O N
100
I
E 0.eL
I
/
E o
i
6o
.~.
40
o
2O -> 0.0 ~ oo
0.2
o.~
o.6 x (ff~)
~ " ~ "- "~ ~ -
0.8
04
06
0.8
1.0
Fig. 4. The energy relaxation time % (left axis) and the average power dissipation -C~ (right axis) of the test device. Energy relaxation time is calculated from % = - ( W - Wo)/Cw and the Monte-Carlo data are noisy only in the low-field regions.
4.1, Mobility and energy relaxation time
RESULTS
A well-known nonlocal phenomenon is described by the average velocity profile in Fig. 2. Not only does the velocity increase before the electric field starts to ramp up (anticipatory effect[27]), but an obvious peak (overshoot) and a less significant dip (undershoot) are also observed. Also shown in Fig. 2 are the average energy and the average energy flow. No overshoot phenomena for these two moments are observed. It is worth mentioning that although n and V are not uniform individually, the product nV is constant throughout the region, since no generation or recombination of electrons is considered.
]>
02
I .,0
0
Fig. 3. The mobility # (left axis) and the average frictional force Cp (right axis) of the test device. Mobility profile is calculated from expression p = - V/Cp and the MonteCarlo data are noisy only in the low field regions.
The two collision moments Cp and Cw are calculated directly by the Monte-Carlo method. Details are discussed in the Appendix. Using Cp, C,. and the calculated V, W, the spatial variation o f p and % can be obtained. Figure 3 shows the profiles of Cp and p and Fig. 4 shows the calculated values of C,, and %. Although C, and Cw follow the electric field shape in general, a spatial delay is observed, especially in the region where the electric field increases rapidly. This clearly demonstrates the nonlocal transport property of the electron gas where the response lags behind the driving field. It is more complicated to examine the profiles of/~ and ~,. Although the data are noisy in the low-field regions, excellent results are obtained in the high-field and field-transition regions. These data are used subsequently to examine the possible models for p and ~,. Three possible mobility models are compared in Fig. 5. Figure 5a shows the reciprocal mobility p - ~vs electric field. The hysteresis loop clearly indicates that p is not a single-valued function of the local electric field, Compared to homogeneous field data, kt is greater if dJE[/dx is positive and smaller if d[E[/dx is negative. If/~ -~ is plotted against W, as shown in Fig. 5b, there is still a hysteresis loop. Although the homogeneous-field data approximate the data in the positive dlEI/dx region, they still overestimate the value of/~ in the region with a negative dlEWdx. This implies that the nonstationary transport phenomena cannot be properly described by any p =/~(W) model. If, however, the moment expansion suggested in Ref. [28] is followed, i.e. to express C'p as:
Cp = - q ( a V + bS),
(31)
where a, b are the expansion coefficients, then:
p -~ = a + b(S/V).
(32)
566
S H I N - C H I LEE a n d
'~
•
homog
c,'
inhomog
o
a
6
°
•
~
ij
"
oe
2">
(C) "
20
r}
4tt~
1+~-2 E 2
1/2 E2
(36)
and the impurity-concentration dependence can be included conveniently by adopting any empirical form of P0 (NI). A similar approach can be applied to the modeling of the energy relaxation time. If %. is plotted against W, as shown in Fig. 6, only one path is followed for all data. Homogeneous field calculations also show a similar relationship. Hence, an empirical expression r , = z,,(W) can be applied to fit all data. A fourparameters empirical form:
J
o
_o o
1+
o o
o
W=W°+q%
t
o
7~
TING-WEI TANG
4C, 60 - £ (kV,/cm)
80
! '}
0
8
o9!
.....
q'~omog
o o
o0
!
• o
I
o a _J o)
T~ 2 t
i
o
L,(W)=to+fi(W/Wo-l)
o
o,O
i
,. I' U
+t=exp(-fl(W/Wo-1)),
i
O0
01
0 2
is proposed here, where to, tl, t2 and/~ are chosen to be 0.28 ps, 3 fs, 2.2 ps and 10, respectively. This result basically supports the conventional assumption that % is a single-valued function of the average energy, even in the presence of a rapidly changing field. This % - W relationship is very sensitive to the band structure model assumed in the calculation and a parabolic band model would be an oversimplification [29].
0
(ev) F
6 k
~3
>
fs
o - -
]
=a¢o(S.,'\.,',
~
4
J
O~ O0
O1
02
O}
4.2. Energy tensor and heat flow
04
s / v (ev) Fig. 5. The reciprocal mobility is plotted against: (a) electric field; (b) average kinetic energy W; and (c) ratio between average energy flow S and average velocity V. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous calculation. The arrows in (a), (b) indicate the directions of the hysteresis loops and the solid line in (c) is a linear fit of the data.
Figure 5c demonstrates that this is a much more accurate model of #, since all data follow this linear relationship independent of field variations. The parameters a and b can be determined directly from Fig. 5c, or they can be expressed by the more familiar parameters: low-field mobility P0 and saturation velocity %, as: a
1
w0
(33)
l~o q~:,,v~'
In general, the energy tensor U has six components; however, for the I-D case, ~/is a diagonal matrix and only the longitudinal component Ut and the transverse component Ut need to be considered. Figure 7 shows that UI is greater than Ut, especially in the high-field region, indicating that the distribution function is elongated along the field direction and the
b = ~ ,
(34)
•
homog
I
o
inhomoq,
i
---G? D_ .J
fittir, g
£;6
'2 I
[
1
r
'., ~
1
(37)
v _ r ~ _ 4~.~q~_~_ ~
02
i t
where S is a parameter to be discussed in the next section. This transformation gives the electric field dependence o f the average velocity and energy of bulk Si as: =
- 2/~0 4,2
1+
l+~-E2)
\1.2
~,
(35)
00 00
0 '
02
I 0.3
w (ev) Fig. 6. The energy relaxation time plotted against the average kinetic energy. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous field calculation. The solid line is an empirical fit of the data.
Silicon hydrodynamic model
• 02
U
o
Ut
--
2w/5
.7 o 0~
~J
,l ° °>
<~ !
,,/o o/o Ic,
i
% [ O0
i
I O0
O2
O~
0 6
08
567
observed in a uniform-field condition, where Q is always in the opposite direction of S. This is a direct manifestation of the change in the shape of the distribution function. The contribution of Q to S is significant in the high-field region, mainly due to Q .... term, and neglecting Q would give an almost 30% overestimation of S in the region. The closure requirement of the H D E still needs an expression of the energy flux density nS by the lower-order moments. If we can assume that the fourth-order moment R has an quadratic dependence on W and ~q is proportional to #[20], then Qdiff may be approximated by a form similar to the Wiedem a n n - F r a n z law: 7~
10
< ,:,~m)
Qdiff =
Fig. 7. The longitudinal component Ul and the transverse component Ut of the tensor U. The solid line is two-thirds of the average energy W which approximates U~fairly well.
equipartition assumption of the energy is invalid. However, if 2 W/3 is also plotted in Fig. 7, we can see that this quantity approximates G fairly well. This is because the effective mass m, used in the definition of Wis greater than the effective mass used in ~, which effectively simulates the effect of the elongated distribution function on U~. Other empirical models of G are also possible[20]. No attempt is made to model U,, however, since it does not affect the electron transport for 1-D problems. The total Q in the device is plotted in Fig. 8, in which the convective part Qconvand the diffusive part Qaifr are identified seperately. Q .... is the dominant term in the uniform field regions where all the spatial gradients are negligible while Qdierdominates over the field ramp regions. It is interesting to see that Q is not always negative. Near x = 0.7 #m, Q is actually in the same direction as S. No such phenomenon is
(38)
dx
S = ( W + k B Te, ) V + Q ..... +
Qdiff,
(39)
,~ g V W "}- Qdiff,
(40) dW
~-SVW-
~W--
q
(41)
dx
Figure 9 shows S plotted against V W taken from the Monte-Carlo data for both homogeneous and in-
!i '3'.}
F ?
•
homcg
o
nbomoc
/~o,@/~" Y. ,61" /
%-
O0
q
The parameter A differs from the conventionally defined Lorenz number by a factor of (3/2) 2, since an energy formulation, rather than a temperature formulation[3], is adopted here. The homogeneous field calculations suggest that ( 2 ~ - ( 1 / 3 ) W V [ 2 0 ] and hence, a simple expression of Q ..... proportional to VW, is assumed. If the convection energy flow ( W + k B T~)V is also assumed to be proportional to VW, then it can be combined with the Q .... term and the energy flow S may now be expressed as:
01
52 E o
dW
-- --/Z W - - .
. . . . }'vw
/
,
,xs~t.." ,
02 /
o
• ,/
/ Dq
o
t/
© {
] ¢I
/
SI - 0®~
Q ~ ;
I
....
0.0 0.0
01
0.2
03
VW , ' 1 0 7 e V * c m / s )
I
02 ]
O0
02
0.4
0.6 x (~)
0.8
10
Fig. 8. The total average heat flow Q ( ) and its two components: convective heat flow Q.... ( 0 ) and diffusive heat flow Qai~r(0).
Fig. 9. The energy flow S plotted against the product of average velocity V and the average kinetic energy W. Dark circles are from homogeneous field calculations and open circles are from the inhomogeneous field calculation. Dashed line is 6 VW, calculated from the homogeneous field data, and the solid line is o~VW- z~#W V W/q, calculated from the inhomogeneous field data. ~" and z~ are 1.29 and 0.2, respectively.
568
SHIN-CHI LEE and TING-WEI TANG
homogeneous-field conditions. For the uniform field case, S is predicted fairly well by selecting ~ = 1.29. This 3 V W term combines the convection energy flow and the convective part of the heat flow. For the nonuniform-field case, eqn (41) qualitatively predicts the characteristics of S with a A of 0.2. These values of the parameters S and ~ correspond to values of 1.95 and 0.45 for the parameters 6 and A defined in Ref. [3], respectively. However, there are still slight deviations of the model from the Monte-Carlo calculated S, indicating that the Wiedemann-Franz law is only a first-order approximation for Qain. Other empirical models of Qd~ are also possible [20]. This set of closed HDE with transport coefficients given by (32), (37) and (41) is applied to simulate real devices with significant nonstationary transport phenomena. Good agreement is obtained between the Monte-Carlo calculated data and the HDE solutions [31]. 5. DISCUSSIONS AND CONCLUSIONS
The derivation of the HDE and associated transport coefficients is given in this study. In order to close the moment equation hierarchy, emphasizing the validity in devices operated under rapidly changing electric field conditions, an inhomogeneous field Monte-Carlo extraction scheme is developed to examine the appropriateness of the transport coefficient models. The final HDE model proposed for device simulations can be summarized as follows:
different definitions of the effective mass. The advantage of this approach is that the form of the equations is simpler and it is more adaptable to any generalization to a real band structure, where no analytical form of d'(Yc) can be found. All the band structure information is included in the empirical expressions of the transport coefficients. It is worthwhile to discuss the relationships between this model and other HDE formulations. If expressions (31), (46) are adopted, instead of using mobility, then (43) takes the form: 2
V(n W) - nf" = - qn (a ~>+ bS) = aJ - qbnS.
Combined with (47), we can obtain: _2V(n W) - nF = (a + b6 W)J + qnpb A W A W. (51) 3 This suggests that an equivalent current density equation can be expressed as[32]: ~)= - n/iF + q/9,, Vn + n/),~ V W,
V
(42)
/~V(nU) - n p F = ),,,
(43)
1_ q
(nS) + - J,,. P =
W - Wo r,
- n -
(44)
where
P= -q(;~+ ;~), 2 nS=gi=
q
(45)
U = - W. 3
(46)
W)~--npWVW, q
(47)
b qnSi Ia--(~-,,,)J
1
,
/i ~ (a + bgW) 1
(53)
2W ~,, ~/~ - - 3q'
(54)
p \ ~ - qb Al~ W).
(55)
~zx'2(.
Equations (43), (52) are equivalent only if the transport coefficients are defined according to the relationships described above. Another form of the I-D current equation that has often been referred to as the augmented drift diffusion current equation[33] can also be approximated from (50). If the field-dependent average velocity, average energy and average energy flow for a homogeneous system are designated as Vh, W h and Sh, respectively, then: - E = aV h + bSh. (56)
(48)
If the nonstationary phenomena of the energy flow are neglected and S is approximated by Sh in (50) then: 2 d V = Vh(E ) (nW). (57) 3qan dx
(49)
If W is further approximated by Wh, (57) becomes:
r,, = to+ t t ( W / W o - 1) + t2exp(-fl(W/Wo -1)).
(52)
if an equivalent energy-dependent mobility /i, an equivalent energy-dependent particle diffusivity 15. and an equivalent energy diffusivity/), are defined as:
D V . ]1,,= qn( R - G),
(50)
Values of parameters 3", A, to, tl, t: and fl are selected as 1.29, 0.2, 0.28 ps, 3 fs, 2.2 ps and 10, respectively. It should be pointed out that models of U, S and r,. are empirical approximations only. Further improvements are possible. One of the considerations for choosing the form derived here is to avoid any explicit average mass terms appearing in the equations. This reduces the number of transport coefficients needed in the equations to a minimum, since there are several
V = Vh(E)[1
L(E) dE dx
D(E)n dx'an
(58)
where the field-dependent diffusivity D(E) and the characteristic length L ( E ) are defined as: D(E)~_ 2Wh(E) ' 3qa L(E) ~
2E dWh 3qaVh(E ) dF
Silicon hydrodynamic model The a d v a n t a g e of this f o r m u l a t i o n is t h a t energy i n f o r m a t i o n is no longer needed for ocity calculation. However, the a s s u m p t i o n p e r t u r b a t i o n from the h o m o g e n e o u s state valid for the weak field gradient regions applicable range o f this expression is very
average the velof small is only a n d the limited.
Acknowledgements--This work was supported in part by
National Science Foundation Grant No. ECS- 8611456. We have also benefited from the General Technology Division, East Fishkill Facility, of the IBM corporation. The authors also wish to thank Professor David Navon for carefully reading the manuscript. REFERENCES
1. W. Fawcett, A. D. Boardman and S. Swain, J. Phys. Chem. Solids 31, 1963 (1970). 2. B. Meinerzhagen and W. L. Engl, IEEE Trans. Electron Devices ED-35, 689 (1988). 3. R. Stratton, Phys. Rev. 126, 2002 (1962). 4. K. Blgtekjaer, IEEE Trans. Electron Devices ED-17, 38 (1970). 5. M. Sanchez, Solid-St. Electron. 16, 549 (1973); W. Htinsch, In Proc. Two-Dimensional Systems: Physics and Devices, p. 296. Springer, Berlin (1986); E. M. Azoff, Solid-St. Electron. 30, 913 (1987); C. C. McAndrew, E. L. Heasell and K. Singhal, Semicond. Sei. Technol. 2, 643 (1987). 6. D. L. Woolard, R. J. Trew and M. A. Littlejohn, Solid-St. Electron. 31, 571 (1988); R. Thoma, A. Emunds, B. Meinerzhagen, H. Peifer and W. L. Engl, I E D M Tech. Dig., p. 139 (1989). 7. C, C. McAndrew, E. L. Heasell and K. Singhal, Semicond. Sei. Technol. 2, 643 (1987); E. M. Azoff, J. appl. Phys. 64, 2439 (1988). 8. M. Rudan, F. Odeh and J. White, C O M P E L 6, 151 (1987). 9. R, K. Froelich, Ph.D. Dissertation, University of Michigan, Ann Arbor (1982). 10. M. Fukuma and R. H. Uebbing, I E D M Tech. Dig., p. 621 (1984); M. Tomizawa, K. Yokoyama and A. Yoshii, IEEE Trans. Comput. Aided-Des. CAD-7, 254 (1988). 1I. R. K. Cook, IEEE Trans. Electron Devices ED-30, 1103 (1983); G. Baccarani and M. R. Wordeman, Solid-St. Electron. 28, 407 (1985). 12. R. Bosch and H. W. Thim, 1EEE Trans. Electron Devices ED-21, 16 (1974). 13. W. R. Curtice and Y. H. Yun, IEEE Trans. Electron Devices ED-28, 954 (1981); R. K. Cook and J. Frey, IEEE Trans. Electron P - ~ ~,9, 970 (1982); C. M. Snowden and D ,s. Electron Devices ED-34, 21 ~ 14. E. ~"
evices ED-36, 609
t, IEEE Trans. "ctron Devices
koyama and TAD-7, 254 ai, IEEE IEEE
569
19. T. J. Bordelon, X. L. Wang, C. M. Maziar and A. F. Tasch, I E D M Tech. Dig., p. 353 (1990). 20. S.-C. Lee, Ph.D. Dissertation, University of Massachusetts, Amherst (1990). 21. A. Bringer and G. Sch6n, J. appl. Phys. 64, 2447 (1988). 22. R. Thoma, A. Emunds, B. Meinerzhagen, H. J. Peifer and W. L. Engl, 1EEE Trans. Electron Devices ED-38, 1343 (1991). 23. N. W. Ashcroft and N. D. Mermin, Solid State Physics. Saunders College, Philadelphia (1976). 24. C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983). 25. M. Shur, Electron Lett. 12, 615 (1976); R. Thoma, A. Emunds, B. Meinerzhagen, H. J. Peifer and W. L. Engl, IEEE Trans. Electron Devices ED-38, 1343 (1991). 26. Y. J. Park, Private Communication. To be published. 27. S. Bandyopadhyay, M. E. Brown, C. M. Maziar, S. Datta and M. S. Lundstrom, IEEE Trans. Electron Devices ED-34, 392 (1987). 28. W. H~insch and M. Miura-Mattausch, J. appl. Phys. 60, 650 (1986). 29. S.-C. Lee and T.-W. Tang, Computational Electronics, p. 127, (Edited by K. Hess, J. P. Leburton and U. Ravaioli). Kluwer Academic, Boston (1991). 30. P. J. Price, I B M J. Res. Dev. 14, 12 (1970). 31. S.-C. Lee and T.-W. Tang, Proc. 1991 VPAD, p. 47 (1991). 32. R. K. Cook, IEEE Trans. Electron Devices ED-30, 1103 (1983). 33. K. K. Thornber, IEEE Electron Device Lett. EDL-3, 69 (1982).
APPENDIX
Collision Moment Calculation
In Monte-Carlo calculations, each spatial bin is treated as one slab of bulk material with a constant electric field. There are two fundamental processes taking place in each bin: free flight and scattering. Noticing that the collision integral on the right-hand side of the BTE is basically the difference between the A-state distribution and the B-state distribution[30], a direct approach is developed in the following to calculate the average collision moments of each bin. The ¢(fQ-collision moment in ith bin, C~j, is defined as the rate of change of function ¢ due to collisions taking place in this bin. If the states of the particle before and after the j t h scattering event are indicated by f% and ~aj, respectively, then the change of ~ due to this j t h scattering is simply ¢(k~j) - ¢(/%). Hence, the value of C~ at this ith bin can be calculated by:
C,~,,-
Ati
(61)
where the summation runs over all the scattering events taking place in the ith bin and A6 is the total flight time the particle spends in this bin. If function ¢ is substituted by hfc or g, 6"p or C~ of each bin can be calculated directly, and the spatial variation of ~p and C, are obtained. Since Monte-Carlo results automatically satisfy all the moment equations, the spatialgradient terms of the moment equations can also be calculated from the differences between the collision terms and the field terms. This approach avoids the use of the conventional finite-difference scheme on the noisy Monte-Carlo data to calculate the spatial-gradients and a much smoother function over the coordinate space can be produced.