- Email: [email protected]
An efficient deterministic solution of the space-dependent Boltzmann transport equation for silicon
Solid-State ElectronicsVol. 35, No. 1, pp. 33~,2, 1992 Printed in Great Britain. All rights reserved
0038-1101/92 $5.00+ 0.00 Copyright © 1992 PergamonPress pie
AN EFFICIENT DETERMINISTIC SOLUTION OF THE SPACE-DEPENDENT BOLTZMANN TRANSPORT EQUATION FOR SILICON HONGCHIN LIN, NEIL GOLDSMAN and I. D. MAYERGOYZ Department of Electrical Engineering, University of Maryland, College Park, MD 20742, U.S.A.
(Received 17 April 1991; in revised form 27 June 1991) Abstract--The space-dependent nonequilibrium distribution function is obtained by directly solving the Boltzmann transport equation for Si. Legendre polynomials are used to express the Boltzmann equation in an analytical form which incorporates the effects of a nonparabolic band structure, inelastic phonon scattering, as well as spatial variation. The resulting linear difference/differential equation is discretized and solved using sparse-matrix methods in the energy domain, and iterative methods in the real-space domain. Values for average electron energy and electron drift velocity calculated from the resulting distribution function are in good agreement with those obtained from Monte-Carlo calculations. This method should prove to be useful in CAD tools for semiconductor devices.
1. INTRODUCTION
has also been applied to the quantum-mechanical regime to obtain degenerate distribution functions[6]. Here, the work is extended into the real-space domain by including the spatial derivative of the Boltzmann transport equation in the analysis. As in Ref. [4], our transport model[3] includes the effects of nonparabolic band structure, as well as acoustic and inelastic intervalley p h o n o n scattering. To solve the space-dependent Boltzmann equation, we first use analytical methods to obtain a linear, second-order partial differential/difference equation. This form of the Boltzmann equation is then discretized and solved numerically by employing a sparse-matrix method in the energy domain, and an iterative method in the real-space domain. Once the S D M D F is ascertained, we use it to compute average electron energy, drift velocity and carrier concentrations as functions of position. Calculated values are in excellent agreement with results from Monte-Carlo simulations.
We present here a new efficient technique for calculating the space-dependent distribution function in Si. The method requires relatively little computation time, and can therefore be useful for everyday use in the computer-aided design of semiconductor devices. The space-dependent m o m e n t u m distribution function ( S D M D F ) provides the probability density of finding an electron with a specific electron wave vector at some point in a device. Information provided by the S D M D F is essential for accurate analysis of many semiconductor phenomena, such as M O S F E T degradation, E P R O M programming and H E M T performance. Traditionally, the distribution function has often been taken as a Maxwellian[1,2]. However, the Maxwellian approximation is only accurate for very low electric fields. To obtain the distribution function for high electric fields, investigators usually use Monte-Carlo methods. Monte-Carlo techniques are stochastic, and thereby often require extremely large computation times to obtain statistically reliable values for the distribution function. To more efficiently obtain the distribution function, we have developed a deterministic method, which requires much less computational time than the traditional stochastic approach. In this paper we present the new deterministic method for ascertaining the S D M D F . The method employs the same transport model used by many Monte-Carlo calculations[3], and requires relatively little CPU time to calculate. To obtain the S D M D F , we build upon an existing technique for solving the Boltzmann transport equation which combines analytical techniques involving Legendre polynomials with numerical methods. This new method has already been used for obtaining values for the space-independent distribution function which agree with MonteCarlo simulations[4] and experiment[5]. This approach
2. ANALYTICALFORMULATION
2. I. Formulation o f transport equation
space-dependent
Boltzmann
To obtain the space-dependent m o m e n t u m distribution function ( S D M D F ) for Si, we begin with the standard Boltzmann transport equation for electrons in steady state: 1 e [Of(k, r ) ] hV~E'Vff(k'r)--hE(r)'V~f(k'r)=l_ 8t dc' (1) where k = electron wave vector, E = electron energy, r = electron position vector, f ( k , r ) = e l e c t r o n S D M D F in one of the conduction band valleys, E(r) = space-dependent electric field. 33
first
HONGCHIN LIN et al.
34
The first term on the left-hand side is often considered to give the variation i n f ( k , r) which is due to diffusion; the second term gives the effect of electric field; the term on the right-hand side is the collision integral which describes the effects of the major scattering mechanisms in Si--acoustic and intervalley phonon scattering. Therefore, we rewrite the right-hand side as:
~j~
L~ja
c L~jI
v.
(2)
In this equation the subscripts ac and iv correspond to acoustic and intervalley phonons, respectively. To solve eqn (1), while also accounting for the effects of Si's ellipsoidal energy surfaces, we have to first transform eqn (1) to the starred momentum space (k* space), where the constant energy surfaces are spherical[7]. To facilitate our space-dependent solution, we perform the inverse transformation to the position vector r. The starred electron wave vector k* and the starred electric field vector E*, as well as the starred position vector r* have the following relationship to their original values: k* = ~ l / 2 k ,
(3a)
E* = ~l/2E,
(3b)
In this well-known dispersion relation, /3 is the nonparabolicity factor which has been given the value /3 = 0 . 5 e V ~ for Si[3]. While we have used this popular expression for the band structure, it is important to note that only 7(e) will be used in our formulation. This alleviates restrictions on allowable dispersion relations to continuous analytical functions of energy. To facilitate our formulation of eqn (4) into a solvable expression, we choose the electric field to be in the symmetrical[l 11] crystallographic direction[4]. This field direction forces the S D M D F s to be identical in each valley. Furthermore, the position vector is taken to be I-D and parallel to the electric field E*(r*). Under this condition, the problem will be 1-D in real-space, and E* and r* can be replaced by the scalars E* and x*. The next step is to express the S D M D F as the sum of the first and second Legendre polynomials multiplied by their respective coefficients[4, 8]: f(k*, r*) = f ( k * , x*) =.f0(< x*) + k'g(< x*)cos 0, (6) where
(3c)
r* = ~ 12r,
0 = angle between E*(x*) and k*, f0(e,x*) = coefficient of the first Legendre polynomial which is the symmetrical part of the SDMDF, k*g(e, x*) = coefficient of the second Legendre polynomial which is the antisymmetrical portion of the SDMDF.
where the transformation matrix fi~/2 has the following form:
x/mo/m,
~1:2 =
0
.
o
Here, m 0 is the free electron mass. After transforming the Boltzmann transport equation, we obtain: 1
e
Vk,e . V,.f(k*, r*)- ~ E*(r*)- Vk,f(k*, r*) _
L
at
r*)j.I
(4)
We use the following dispersion relation to account for Si's nonparabolic band structure: h2k .2
7(E) = E +/3E ~ =
2 7x/~%~'a2Jo(e,x *) 37,(e)Xmo( ~XX,7
(5)
2m0
To ascertain the S D M D F we must now determine the two unknown functionsf0(e, x*) and g(< x*). By substituting eqn (6) into eqn (4), we can obtain a tractable form of the Boltzmann equation, which contains both symmetrical and anti-symmetrical parts. By equating symmetrical parts to symmetrical parts, and anti-symmetrical parts to anti-symmetrical parts, and performing algebraic manipulations, two equations are obtained (details in Appendix A). These two equations can be solved for the two unknown functions. The equation for the symmetrical part of the S D M D F .1~(~.x*) is:
•2Jo(e,x*)
2eE*(x*)
I , #fo(E, x*) x eE*(x ) ~
&ax~
t-[eE*(x*)]
?~,2Jo(e,x*)] 2 1 ~ ~+~ 1
efo(,, x * ) ] _ 2x/~(,)ro e aE*(x*) OJa(< x*) + ax* J 37'(e)3mo ax* a~
27 (_e)? "(e ! l roeE*(x*) ?,(e)2 jT,(e)2x/7(e)mo
4xf2m~m? '2D~¢ ~ > " ( e ) 2 rch4p
D .(mt mt) 1 +,,=, x/2=h3peo,, exp(hog.[-KoU)- 1 {x/Y(e + ho),,)7'(e + hog.)[exp(ha~,,/KoT)fo(e + hoJ.,x*) - f o (e, x *)] + x/Y (e - hog.)7 '(e - hoo,,)[fo(e - ha~., x *) - exp(h~o./KoT)fo (e, x *)]} = 0,
(7)
Boltzmann transport equation for Si where r 0 = proportionality factor which relates the p h o n o n scattering rate to energy, z0 Has units of eV m/2s and its values are given in Appendix A, 7'(e), y"(e) = first and second derivatives ofT(e) with respect to energy, respectively, mt, m~ = Si's transverse and longitudinal effective masses, respectively, Dac = acoustic phonon deformation potential, p = density, T = lattice temperature, K0 = Boltzmann's constant, D,, = intervalley phonon deformation potential and subscript n stands for different kinds of intervalley phonon scattering, o~,=intervalley phonon vibrational frequency. The second-order differential/difference equation given by eqn (7) represents a boundary value problem in two dimensions: real-space and energy space. The four sets of boundary conditions used to help define the problem are depicted in Fig. t and are described below: 1. The distribution function at x = 0, which corresponds to the region where electrons are injected from an ohmic contact, is given by a Maxwellian. 2. At x = / , which is the right edge of the device, there we can use either a Dirchlet condition where the S D M D F is specified, or a N e u m a n n condition with the normal, real-space derivative of the distribution function equal to zero. In the presented calculations, the Dirchlet condition is used, with the requirement that the distribution function is equal to its homogeneous-field limiting value. 3. The S D M D F is zero at very large energy, which for our calculations is approx. 1.5 eV.
35
4. Since we are assuming there is no generation or recombination, the product of electron concentration and drift velocity n(x)v(x), and therefore current density are constant. The value of current density depends on the doping concentration. We use the requirement of constant current density as a boundary condition. In addition to specifying boundary conditions, it is important to recall that we are interested in physical phenomena. However, eqn (7) is expressed in terms of starred-space variables. Since the vectors in eqn (7) have been reduced to scalars, we specify the scalar relationship between E*, x* and their real-space counterparts E, x. These relations, which are derived in Appendix A, are given by: E* =
m~E,
(8a)
x* =
--x,
(8b)
where the conductivity effective mass m c is given by:
1 _ 1 (2+1~. mc
3 \m t
(9) mj
To facilitate our formulation, we first transform E and x to E* and x*. Eventually, after eqn (7) is solved, we will then transform back to E and x.
2.2. Formulation of average electron energy and drift velocity While the S D M D F provides essentially all the information required for performing a comprehensive device analysis, submicron device performance is often evaluated in terms of the average values of carrier velocity, energy and concentration. Furthermore, the product of average electron concentration and drift velocity provides one of the boundary
Yo(¢ > N A ¢ , x ) = 0 = NA¢
l
fo(~, o) , -
......
~
fo(~, l)
¢=A¢
xj-1 xj xj+l
x----0
x= l
n(x)v(x) = g Fig. 1. The 2-D domain for the space-dependent Boltzmann transport equation is shown. The vertical dimension represents the energy domain whereas the horizontal dimension gives the real-space coordinates. At the left boundary the SDMDF is a Maxwellian; at the right boundary the SDMDF is given by its high-field limiting value; at the top the SDMDF is equal to zero; along the bottom the boundary condition is given by a constant value of current density at each real-space coordinate.
36
HONGCHIN L1N et al.
conditions for the numerical solution. Thus, the formulas used for calculating these average quantities from the SDMDF are given in this section. The carrier concentration in terms of the SDMDF is determined first: 1 ('~ n(x*) = 47r - - 3 J0 fo(& x*)h(e) d& (1O) where h(e) is the density of states which accounts for Si's ellipsoidal energy surfaces and nonparabolic band structure: .4f 2 2n
,
.
~9h3n
x/~m,~mT
× EeE,(x,) ~j;(< x*)
n(x*)vT(x*) = nAx*)vXx*) + n~(x*)vy(x*) + nAx*)v~(x *) =
7'(E)2
ef0(
f~ ?(e} [ t?fo(e,x*) ~?fo(&x*)~ [_eE*(x*) ~ ]ae. (15) J0 7 ' ( e ) 2 0e The product of concentration and total drift velocity for electrons in the y-valley and the z-valley are identical to eqn (15), so the average product of total carrier concentration n(x*), and average drift velocity vT(x*) for the entire conduction band is given by the vector sum of the products of concentration and total drift velocity in each valley: x
2x/-2 [m~ml 1
l, 1),0
fo ?'(e) 27(e)
(13)
In contrast to average energy, the drift velocity originates from the anti-symmetrical part of the SDMDF, and is difficult to calculate because it is a vector and depends on which valley the electron is located. Furthermore, the nonparabolic band structure and the ellipsoidal surfaces of constant energy give rise to additional complications. After performing a great deal of mathematical manipulations (details in Appendix B), the x-component of the product of concentration and drift velocity in the x-valley is: n~(x*)v~(x*) =
~mo(m~~'1 my 1_, mjl-%
(11)
It is worth noting that the anti-symmetrical part of the SDMDF does not contribute to the integral for the carrier concentration. By definition, the space-dependent average electron energy is given by: f0~ 1 ef0 (& x*)h (e) de. (12) e(x*) 4n3n(x,) As with carrier concentration, the anti-symmetrical part of the SDMDF does not contribute to the integral for the average energy, To obtain an expression for the product of average velocity and electron concentration, we begin with the definition:
f v(k)f(k, x*) dk.
x/3rn0m~j0 7,(e)2
[c'}f0(e,x*) {~f0(e,x*)-] [ (14c) x LeE*(x*) & ~xxg j d e , To obtain the product of concentration nx(X*) and total drift velocity v'(x*), for electrons in the x-valley at position x*, we take the vector sum of the three components given by the above three equations: n~(x,)v%~c,) = 2 9h3n ~ 2
.
h({) = -Y~- (m7ml)l'2?(e)127'(e).
n(x*)v(x*) = ~ 1
n~(x*)v~(x*)=T 9h3g
(14a)
where
× I eE*(x*) {~fo(e,& x*)
~fo(e, x*) 1 ~ x ; .jd&
(16)
The above equation is in the [111] direction as expected. After taking the magnitude of n (x*)v T(x*), we have: =
2,/5
x eE*(x*)
5fo(& x,) &
{~fo(&x*) 1 ~xx* J d&
(17)
Of course, by dividing eqn (17) by n(x*), which is given by eqn (10), values for average electron drift velocity over all conduction band minima are obtained. 3. NUMERICAL SOLUTION
nx (x*) = concentration of electrons in the x-valley at the position x*, v~(x*) = velocity of electrons in the x-valley for the x-component at the position x*. By using similar analyses, the y- and the z-components of the products of concentration and velocity in the x-valley are:
×
I- _,,t x ..) ~£(e, x*)
ef0(e,x*)7jae,
(148)
With the analytical expression of the Boltzmann transport equation derived and given by eqn (7), the next step is to apply numerical techniques to obtain the solution, which is the symmetrical part of the distribution function f0(& x). The basic strategy we use is to approximate eqn (7) by a system of coupled difference equations, and then solve that system. The first step is to discretize eqn (7) with respect to energy E and position x*. We choose mesh points in the energy domain which are equally spaced with increments of Ae. In the real-space domain we employ a fine mesh in regions where the
Boltzmann transport equation for Si S D M D F is rapidly varying, and a coarse mesh when its spatial variation is reduced. To accomplish the discretization we first divide the energy domain into N points, and the position domain into M intervals with M + 1 points, which results in a 2-D mesh with N × (M + 11 points. Using the central-difference method, we obtain a difference equation for each point or node. The i, jth difference equation is listed below: 2xf~i Zo [ 2(f,j ~ l_--f,;j)
32;3molsj+l(sj+sj+1)
2(f,j--fi,j-~] sj(sj+sj+ 1
37
for f0 (E, x*). The procedure is repeated iteratively for all spatial points until the system converges. Numerically, this involves solving an N x N matrix equation repeatedly. It is interesting to note that this solution technique is facilitated by the physics of the problem. In the energy domain, the various difference equations are strongly coupled over a large energy range by intervalley phonon scattering. With this
2x~i%eE*f+l,j+l-f+l,j
I Ifi-l,j+l-~f/
2.,/-~zo(eE*\ 2 2 ( _27i2;"]_ zoeE* f,j+~--f.j , +~--~-E ) ( f + , , j - 2 f o + f _ , j ) - - ~ 1 2;2 ]x/~2;2mo s,-+sj+-~l -27ieE*+,--eE*_ l 2 X - 32;
x
Sj + Sj+ 1
I+-~-.)]U+L~7;+
~
27i2;"~.-..2]f+,J-f-,,] 1 -- -'~-i2)~dff-,j ) J -2--~E
l+~T,}.-o
j
~
1,j-I
sg+sj+ 1
32;3m0 AE
-~
ro N~/i2 ;2 mo
4x/2rntm, Dac ,G2; ~'14p
+ 27 ; o
A¢2
D2"(m2'rn1)'/25 ~,/T-7 , Fexp(h°~./GT)f+, J f,;] + / ~
.
r f - , . . j - exp(ha~./Ko T)f.f]]
v7'-'.7-'° L-
j;=0,
where f j is equivalent to f0(iAE, x*) with i and j corresponding to energy and position domains; sj and sj+l are equal to x * - x * l and x * + l - x * , respectively. We also assume x* is always less than x*+~ for any j ; AE = very small energy increment step; / . = c l o s e s t integer to ho)./Ae; 7~=2(iAE); 2t+t, = 2(i A• + hco,); Y~= 2'(i Ae) = first derivative of the dispersion relation with respect to energy; 71' = 7"(i Ae) = second derivative of the dispersion relation with respect to energy; E~=E*(x*)= electric field at x*. There is an analogous difference equation for each of the N x M mesh points in the 2-D energy/ real-space domain. Since the energy domain and the real-space domain are divided into approx. 300 intervals and 60 intervals, respectively, this represents a matrix equation with approx. 18,000 rows. A direct solution of this system proved to be computationally inefficient. For this reason, we employ two different techniques simultaneously to obtain a solution. In the real-space domain we use an iterative method (method of lines), and in the energy domain we solve our system of equations directly using a sparse-matrix method which employs the minimum degree algorithm[9]. The method can be understood by first observing that in eqn (18), the distribution function at the point x* can be determined if the distribution functions at the points x * j and x*+l are already known. In other words, if we know f0 (e, x*_ l) and f0 (E, x*+l ), we can solve for f0 (E, x*). Of course we do not know f0 (E, x*+ 1), so we assume its values and then solve forf0(e, x*). The result will be improved values
(18)
strong coupling, a direct solution of the matrix equation is prudent. However, in the real-space domain, the difference equations are coupled only to their nearest neighbors, and rapid convergence is obtained using the iterative method described. Before we provide details of the numerical methods, however, we first implement boundary condition (4), which is described above, into the system of coupled equations. To implement constant-current density as a boundary condition, we replace the first difference equation in our N x N matrix, which is i = 1, by the product of concentration and average drift velocity given by eqn (17). The discretized form of eqn (17) is: 2,~z0
~ 7 ~
J = n(X*)VT(X*) = 3h3~ 2 X / ~ , = , x (eE, f+ L~-
t,J
~ 7;
f~0+-~-- -- "~ J - 1/ ~a E -, Sj"~- Sj + I ]
(19)
where J is a constant for any x*. One may now obtain a better understanding of the numerical techniques employed by considering the following discussion. For a given point x* in the real-space domain, we have N - 1 equations in the energy domain in the form of eqn (18), and one equation in the form of eqn (19), The equations analogous to eqn (18) are enumerated from i = 2 to i = N, while we place eqn (19) in the first matrix row, which corresponds to the energy point i = 1. This set
38
HONGCHINLIN et al.
of equations can be expressed by the following matrix equation:
--~j
I~_I--~,+1F~j+l
I,
.002
I
~
'
I
'
'
'
I
'
'
'
I
'
'
'
i
(20) .01,urn
.001
L,
where
.o0o
If.l% ~f = column matrix
t ¢ O.05,um
] f.i~j
o
0
2
.4
.6
.8
Energy(eV)
F:!
l~Ij, [/
1=
column matrix on the left point of x* with kth iteration, ~ $ ~ = column matrix on the right point of x* with k - lth iteration, 1, Rj ~ J = coefficient square matrices of pjk P~ l and Pf+l:, respectively.
All the values on the right-hand side o f e q n (20) are presumed to be known, so P~ can be obtained using the sparse matrix solution algorithm[9]. Once r~ is obtained, we use it to obtain an improved value for the distribution function at the next point in the space domain P~+ i. We perform this calculation for each real-space point until we reach the end of the realspace domain. The process is then repeated until the system converges• To make this system converge faster, we use the SOR method[10] with the expression given below:
=.,
0 .
-[i ,F)
,-
~i+iP~+;
,
(21)
6 where w is the relaxation parameter. We found that the system converges most rapidly when w takes on values of around 1.4. This technique gives results which agree with Monte-Carlo calculations, and require considerably less CPU time to calculate.
Fig. 2. Electron energy distribution function for electric field value of 100kVcm -t at various positions marked by arrows. S D M D F to have a carrier concentration of 10 ~vcm 3 in the homogeneous-field limit. Figure 2 shows the distribution function at different positions for an electron ensemble which is in the presence of a constant electric field with a magnitude of 100 k V c m -1. It is interesting to note that the distribution spreads over a wider energy range as the electron ensemble drifts toward the far electrode. Eventually, the distribution function reaches its non-Maxwellian homogeneous field limit. In Fig. 3, we replot the distribution function, which is given in Fig. 2, on a logarithmic scale. Figure 4 shows the product of the distribution and the density of states, f0(E, x ) h ( O , at different points throughout the device, for a field magnitude of 100 kV cm 1. From the figure it is clear that electrons obtain more and more energy as x increases, so the peak of this 3-D plot is moving from low energy at small distances to higher energy at larger distances. In Fig. 5 we give the resutls of average velocity and carrier concentration investigations. Values for the electron drift velocity (solid curve) and the carrier concentration (dashed curve) are plotted as a function of position for an electric field of 100 k V c m -1. As expected, the product of velocity and concentration is constant, which satisfies the continuity equation ~ - [n(x ) v ( x )] = O. CX
4. R E S U L T S
To analyze the accuracy of our new method for determining the S D M D F , we performed a series of calculations. In these calculations our semiconductor device is an Si bar of length 0.2 gm, at the lattice temperature of 300 K with a constant applied potential. We assume that the electric field within the bar is spatially constant. Electrons are injected into the bar with an equilibrium Maxwellian distribution at x = 0, as they would be from an ohmic contact. Electrons then gain energy from the field, and the S D M D F varies along the device until it reaches its homogeneous-field limiting value. We normalize the
In addition, we compare our values for drift velocity with those calculated by Monte-Carlo simulations[11] (open circles). Agreement between the two methods is usually within 5%, and both techniques predict velocity overshoot which is characteristic of nonequilibrium electron transport. In Fig. 6 we show values for average energy which were calculated using the S D M D F , As is expected, average electron energy values (solid curves) increase gradually with distance until they reach the homogeneous-field limit. The figure shows the results of calculations for three different electric field magnitudes: 100, 60 and 20 kV cm -1. We performed similar
Boltzmann transport equation for Si 1 0 -~
' ' I ' ' ' I ' ' ' I ......
10-~
I ' ' "_---
?
~
°o
1.5
10 -5
1 0 -e
o~ ~ "~
1.5
-~°°Oooooooooooo
~e C ]9,urr
x
\\
1 0 -~
39
...----
10-8
/"
i
I 0 -g
o o
1 0 -10
.5, co r.)
>
0
,,~1 .... 0
i0 -~
~ ~ ~ ] ~ 0 .2
~ I t 4 6 Energy(eV)
.8
I.E
Fig. 3. Electron energy distribution function on a logarithmic scale for electric field value of 100 kV cm t at various positions marked by arrows.
calculations using the M o n t e - C a r l o p r o g r a m which are plotted as open circles. The figure indicates t h a t the two m e t h o d s are in excellent agreement.
5. C O N C L U S I O N
W e have solved the space-dependent B o l t z m a n n t r a n s p o r t e q u a t i o n to o b t a i n the space-dependent m o m e n t u m distribution function. The m e t h o d , which expands u p o n an earlier space-independent technique[4], uses Legendre polynomials to formulate a space-dependent B o l t z m a n n equation. T h e resulting partial differential/difference form o f the B o l t z m a n n e q u a t i o n is then solved numerically. The numerical
.05
I .... .1
I,,, .15
0 .2
Distance(#m) Fig. 5. The space-dependent average electron drift velocity calculated by the SDMDF (--) is compared with electron drift velocity calculated by Monte-Carlo method (O) for the electric field value of 100 kVcm-L The dashed line is the space-dependent concentration obtained from the SDMDF. technique uses a novel c o m b i n a t i o n of a direct sparse matrix a l g o r i t h m a n d a n iterative technique. In addition to space dependence, the effects o f inelastic p h o n o n scattering a n d n o n p a r a b o l i c b a n d structure are included in o u r formulation. W i t h this a p p r o a c h , the S D M D F s are calculated at least 10 times faster t h a n those o b t a i n e d from M o n t e - C a r l o calculations. W i t h accurate S D M D F s determined, space-dependent values for average electron energy, electron drift velocity a n d carrier c o n c e n t r a t i o n s were o b t a i n e d which agree well with those o b t a i n e d from M o n t e Carlo simulations. In addition to calculating these average quantities, the technique gives the entire distribution function which is i m p o r t a n t for modeling
Fig. 4. Electron distribution function vs energy and distance forms this 3-D picture for the electric field value of 100 kV cm-t. The distribution function has been multiplied by the density of states.
HONGCHIN LIN et al.
40
real-space variable from r* to the I-D position variable x*. The direction of x*, with the unit vector denoted as 2", is parallel to the transformed electric field E*(x*). As was mentioned in Section 2, we divide our solution into symmetrical and anti-symmetrical components. In order to obtain a closed set of equations in this analysis, we decompose k* into k* and k~ which are the wave vectors perpendicular and parallel to E*(x*), respectively. By using the dispersion relation given by eqn (5) we then make the following approximation which has been justified in the space-independent case by comparison with Monte-Carlo calculation[4]:
100kV/cm
60kV/cm E2
.1 0-~-'
V/cm ~ ~_i ~ ~ ~ 1 ~ 05
7(E) = h2k~Z- . 3 2rn0
; ~ ~ i ~ t ~ tiff
.1
.15
.2
Distance(#m) Fig. 6. The average electron energy calculated by the S D M D F (--) is compared with average energy calculated by Monte-Carlo simulations (O) for different electric field values of 20, 60 and 100 kV cm L
(A1)
From the dispersion relation, the derivative of energy with respect to k* can be expressed as another form which is: h2k * Vk. E 7 '(E )m 0 (A2) With the above equation, as well as the chain rule, eqn (4) can be written as an expression in 4-D phase-space:
M O S F E T reliability, E P R O M p r o g r a m m i n g a n d H E M T p e r f o r m a n c e . In s u m m a r y , we h a v e p r e s e n t e d a n a c c u r a t e m e t h o d for solving the B o l t z m a n n t r a n s p o r t e q u a t i o n w h i c h is sufficiently e c o n o m i c a l for daily use in device simulation.
Acknowledgements--The authors are grateful to ShiuhLuen Wang and Ken Hennacy for helpful discussions. This work was supported in part by the Semiconductor Research Corporation under Contract No. 91-DJ-130, and the National Science Foundation under Contract No. ECS9010457. REFERENCES
1. J. Frey and N. Goldsman, IEEE Electron Device Lett. EDL-6, 28 (1985). 2. M. Fukuma and W. W. Lui, IEEE Electron Device Len. EDL-8, 214 (1987). 3. C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983). 4. N. Goldsman, L. Henrickson and J. Frey, Solid-St. Electron. 34, 389 (1991). 5. N. Goldsman, Y. J. Wu and J. Frey, J. appl. Phys. 68, 1075 (1990). 6. H. Lin and N. Goldsman, Solid-St. Electron. 34, 1035 (1991). 7. C. Herring and E. Vogt, Phys. Rev. 101, 944 (1956). 8. E. M. Conwell and M. O. Vassel, Phys. Rev. 166, 797 (1968). 9. A. George, J. Liu and E. Ng, SPARSPAK User Guide: Waterloo Sparse Linear Equation Package. Department of Computer Science, University of Waterloo, Ontario, Canada (1980). 10. W. F. Ames, Numerical Methods for Partial Differential Equations, 2nd Edn. Academic Press, New York (1977). 1l. N. Goldsman and J. Frey, IEEE Trans. Electron Devices 35, 1524 (1988). APPENDIX A
Derivation of Boltzmann Transport Equation To derive the space-dependent Boltzmann equation [eqn (7)], we start from eqn (4) which is the Boltzmann equation transformed into the starred phase-space, which has six dimensions: three in momentum space, and three in real space. A 6-D analysis is indeed cumbersome, so we limit our real-space analysis to one dimension, and we change our
hk* "2* 0 f ( k * , x * ) 7'(E)m 0 ~x*"
x
eE*(x*)Y:*
hk* ~ , [-~f(k*,x*) 1.~t y,(,)m0~J(k*,x*j=| J
. (A3)
L
To facilitate our formulation, we limit our specification of momentum space to directions perpendicular and parallel to the electric field. We then use Legendre polynomials to express the distribution function as the sum of its symmetrical and anti-symmetrical components, and then substitute this expression, which is eqn (6) with k * = k * cos 0, into eqn (A3). After some mathematical manipulations, the following expression is obtained:
hk* [afo(e,x*) +k,~g(E,x*)l y '(E )mo L- Ox *
" ?~X
*--J
eE*(x*) [ h2k* cJo!~: x*) h LT'(Qmo & +g(~, x*) h 2k~':
+ 't'(~o
Dg(E, x*)l = [g.fo(<, x*) ] k ' g ( , , x*) <~ J L at j~ ,({)
(A4)
where the collision term on the right-hand side of eqn (A4) has been separated into its symmetrical and anti-symmetrical components; k~ and E*(x*) are the magnitudes of k* and E*(x*), respectively; and l/r (E) is the summation of all acoustic and intervalley scattering rates. In Ref. [4], r(E) is shown to be well approximated by: r(Q =
-
TO
~ ,
(A5)
;"(EICg(E) where z0 represents all the coupling factors for deformation potential scattering. The approximate values of l/z 0 are the slopes o f the curve obtained by plotting total scattering rate vs y'(E),~314,51: r 0 = 4 . 2 x 1 0 - t 4 s . e V t'2, ( <0.055eV, % = 2 . 1 x 10 las.eWl'2.
E >0.055eV.
By equating the terms of the same Legendre polynomial order, and using eqn (A5), we eliminate the independent variable k~, and obtain the following two equations from which to determine the unknown functions fo(E,x*) and g(E, x*):
~'(Qm0 L ax*
eE*(x*)
0T:, ]
Y'(Q 7 ~ g ( E , x * ) , %0
(A6a)
Boltzmann transport equation for Si 1 ~'27(e)
Any vector x*2* can be transformed to the original space by eqn (3b). The result is:
Og(E,x*) eE*(x*)[g(e,x*) .
1 1%
(A6b>
+ 3y'(e) & J J L Ot Jc To eliminate g(E, x*) from the above two equations, we substitute eqn (A6a) into eqn (A6b) and obtain an equation which only contains the symmetrical part of the distribution function: 27 x / ~ % _a fOfo(e,x*) 3y'(Q3m0 Ox* L Ox*
eE* x* a [
+ 27(Q%
41
'
which is not in the [111] direction and depends on which valley the electrons are located; for example, if the
eE*(x*)~] I
eE*(x*)z o f. of° (E' x*) +7'(Qzx/y(Qm0l-- ax*
\ m x'my'm~/
~?fo(,,x*) ee*(x*) Of0(,,x*)]
eE*(x*)
= dE
J
(A7) I
ot
_1~
Upon taking the derivatives of eqn (A7), the final form is: 3y,(e)3m° [ +~[1
2eE*(x*) & O~xT - + [eE.(x.)]2 O
~XO2
27(,,7"(e!]
%eE*(x*> ~Ofo(,,x* , e E ' ( x ' , ~ ]
~'(E)= JT,(,)~moL COX*
+ 3Y'(E)3m0
&
Ox*
=foeo¢,,x*q L
at
J¢"
(A8)
The acoustic and the intervalley phonon collision terms have already been given in a previous work[4], and are rewritten below to facilitate our presentation: c~t
nh4p
Jac
~y'(e)
2
1+
f0(E, x*)
+f 7(')+(1
L~
+Y(')Y"(')XIK T] cof°(''x*) 7(QK°TCO2f°(e'x*)'I ~ ) o ~ -~ ~ ~ -j (A9)
and
Ofo(,,x*) at
_ ~ D2.(rnZtm,)',2f
1 .,/2~h3p~o. L e x p ( ~ ) - I J
3,-.-%
+ he).) [exp(R.)f0 (E + h~o., x*)
-- fo(E, x*)]
+ x / ~ - he).)y '(E - hco.)[fo(E -- he)., x*) - exp(R.)f0(E , x*)]}, (A10) where R. = ho).lKo T. Finally, eqn (7) is obtained by taking the summation of these scattering mechanisms, and moving the left-hand side of eqn (A8) to the right-hand side. After eqn (7) is solved numerically, the distribution function, which is in the starred space, has to be transformed to the original space. First of all, E(x*) can be expressed as (I/x/3)(1, 1, 1)E(x*) since the electric field in the [Ill] crystallographic direction, so by using eqn (3b), the relation between E*(x*) and E(x*) is: E*(x*) =
,
-~f]E(x*)
,m/
electrons are residing in the x-valley, the position direction will be: ' /T/t '
However, electron transport in Si is governed by the combined effect of motion within tile six valleys in momentum space. By taking the vector average of the x-direction in each valley, the value becomes: m.,/~om~{ 1 2", ° "/-- +--]x*.
(A1 la)
3
\m I mt/
Thereby, the relation between x and x* is:
E*(x*) = / m° E(x*),
(Allb)
where E(x*) and E*(x*) are the magnitudes of E(x*) and E*(x*), respectively. As was discussed above, the unit vector 2* is parallel to E*(x*), thus 2" is given by: ., X
SSE 35/1--D
/-~¢f =
~-
1 -
1 ,
1
x="
m,/ vom" [r- 1- + - - l x * = 3
k,m,
mr,/
mUo
/ Vx*. ~/m c
(A13)
APPENDIX B
Derivation of Average Electron Drift Velocity To formulate the expression of average drift velocity, we begin with the definition of product of concentration and
HONGCHIN LIN el
42
therefore
drift velocity from semi-classical kinetic theory: . (x *)v(x *) =
1 4-~
fv(k)f(k,x*)dk.
(B1)
By substituting the instantaneous electron velocity v(k) = 1/h V ~ into the above equation and transforming k to the starred space, the preceding expression becomes: 1 ('~1,2
n(x*)v(x*) = 4rc~ J ~
V~.~f(k*, x*)]fi ~'21dk*,
(B2)
where I~t t21 = x/mxm,m:/m 3 is the Jacobian determinant. Upon using eqn (A2) and the distribution function expressed as the sum of the symmetrical and the antisymmetrical parts, we have: 1 ,, . I~l/2hk * n(x*)v(x*) = 47z~v' m, m~m:/m~ j ~ [fo(~,x*)
~\/'mcnvm./rn ~ I
x k* ' g(E, x*) dk*,
(B3)
where k*g(E, x*)cos 0 has been expressed as a dot product, and the vector g(E, x*) is parallel to the electric field E*(x*) in starred space. Furthermore, the integral of k*f0(E,x*) goes to zero because k* is odd and f0(E, x*) is even with respect to k*. To avoid the writing of lengthy formulas, the vector form of eqn (B3) is considered by its x-, y- or z-components separately. We write the x-component of n(x*)v(x*) as:
h x/m~.m. (" k, n ( r * ) G ( x * ) - 47r 3 m ~
)
+ k~g~.(~,x*) + k:g:(~, x*)] dk* -
(B6) &
0x*
(B4)
ing to spherical coordinates, we obtain:
2T0
CvJ,( ,
n(x*)G(x*)=~4
3m, " J T ' ( ~ ) ~
x[eE*(x*)Of°(~ x*)
gf°("x*)~k*2dk*~x*J
f
x leE*(x*) L _2x/2mtml
~/;(e, x*) al0(E,x*)l - "~ I d~ a~ ~x* J
/~
,, '(e )2x/~' (~:)m0
J'
cE
9h3x 2
F
g(~, x*) = g(~, x*)2*
....
×[_e~-tx*;
0S0(,,x*q ~x*
J'
n(x*)
~ .
(B9)
where G(x*), n~.(x*) and n:(x*) are the electron concentration in x-, y- and z-valleys, respectively. Upon modifying eqn (B8) with G(x*)=n(x*)/3, the x-component of the product of concentration and velocity in the x-valley is:
Since g(~, x*) is equal to 2* • [g(~, x*)2*], the vector g@, x*) can be expressed as:
&
7(0 ./,(~)2
?,/;(~ x*)l :"~'.', / d e ax ]
(BIO)
Using a similar approach gives the /-component of the product of concentration and velocity in the ./-valley:
C~fo(e,x*)l ax*
S
r0
4 3morn'
L
•
. (B7)
~ f~ 7(E_!. 3h3x 2 ~/ m, ~/3moJo ?"(e) 2
F ?Jo(,E x*) ×/eE*(x*) ~i
hz0
[eE*(x*)ifo(~x*)
J"
Next we substitute gx(E,x*) into eqn (B4), then using k 2 = 2moy(E)/3h2, which is identical to eqn (A1), and chang-
*
h x/Im~mz P b.~, |",'g4,,x*)dk*, 4 rr~ m2 j 7 (E)
x
tx ;
G(x*)=nr(x* )=n:(x )=
where k,, k,. and k: are the components of k*; the integrals of k~k~g~(Eix*) and k~k:g:(E,x*) have vanished due to the orthogonal properties of k~ and ky, as well as k~ and k~. Next, the appropriate expression of gx(E, x*) has to be determined. By recalling eqn (A6a), we can write g(E, x*) in terms of the symmetrical part of the distribution function:
g(~, x*)
~ / e ~ z
q3m~L
0
Since we chose the electric field to be in the [111] crystallographic direction, the distribution functions are identical for each valley. Hence:
*
"iT'(,) [k~g4''x
x
hz o
'(E )2y ~ m
n(x*)v~(x*)= 2x/F2z° ~
Jy'(Om0
-
gx(E,x*) =
By making use of the dispersion relation, the integration with respect to k* can be changed to the following integration over energy ~:
+ k* . g@, x*)] dk* , . ~"~l'2hk* 4u
al.
(B5)
~[ 3mom, Jo )"(E) 2 0J'°(e' x *)
?-
~/;)(" x * ) l d,
~x*
]
(BI1)
where i = x , y or z a n d j = x , y or z. Finally, the average drift velocity can be obtained from this formula by dividing by the electron concentration as described in Section 2.2.
0038-1101/92 $5.00+ 0.00 Copyright © 1992 PergamonPress pie
AN EFFICIENT DETERMINISTIC SOLUTION OF THE SPACE-DEPENDENT BOLTZMANN TRANSPORT EQUATION FOR SILICON HONGCHIN LIN, NEIL GOLDSMAN and I. D. MAYERGOYZ Department of Electrical Engineering, University of Maryland, College Park, MD 20742, U.S.A.
(Received 17 April 1991; in revised form 27 June 1991) Abstract--The space-dependent nonequilibrium distribution function is obtained by directly solving the Boltzmann transport equation for Si. Legendre polynomials are used to express the Boltzmann equation in an analytical form which incorporates the effects of a nonparabolic band structure, inelastic phonon scattering, as well as spatial variation. The resulting linear difference/differential equation is discretized and solved using sparse-matrix methods in the energy domain, and iterative methods in the real-space domain. Values for average electron energy and electron drift velocity calculated from the resulting distribution function are in good agreement with those obtained from Monte-Carlo calculations. This method should prove to be useful in CAD tools for semiconductor devices.
1. INTRODUCTION
has also been applied to the quantum-mechanical regime to obtain degenerate distribution functions[6]. Here, the work is extended into the real-space domain by including the spatial derivative of the Boltzmann transport equation in the analysis. As in Ref. [4], our transport model[3] includes the effects of nonparabolic band structure, as well as acoustic and inelastic intervalley p h o n o n scattering. To solve the space-dependent Boltzmann equation, we first use analytical methods to obtain a linear, second-order partial differential/difference equation. This form of the Boltzmann equation is then discretized and solved numerically by employing a sparse-matrix method in the energy domain, and an iterative method in the real-space domain. Once the S D M D F is ascertained, we use it to compute average electron energy, drift velocity and carrier concentrations as functions of position. Calculated values are in excellent agreement with results from Monte-Carlo simulations.
We present here a new efficient technique for calculating the space-dependent distribution function in Si. The method requires relatively little computation time, and can therefore be useful for everyday use in the computer-aided design of semiconductor devices. The space-dependent m o m e n t u m distribution function ( S D M D F ) provides the probability density of finding an electron with a specific electron wave vector at some point in a device. Information provided by the S D M D F is essential for accurate analysis of many semiconductor phenomena, such as M O S F E T degradation, E P R O M programming and H E M T performance. Traditionally, the distribution function has often been taken as a Maxwellian[1,2]. However, the Maxwellian approximation is only accurate for very low electric fields. To obtain the distribution function for high electric fields, investigators usually use Monte-Carlo methods. Monte-Carlo techniques are stochastic, and thereby often require extremely large computation times to obtain statistically reliable values for the distribution function. To more efficiently obtain the distribution function, we have developed a deterministic method, which requires much less computational time than the traditional stochastic approach. In this paper we present the new deterministic method for ascertaining the S D M D F . The method employs the same transport model used by many Monte-Carlo calculations[3], and requires relatively little CPU time to calculate. To obtain the S D M D F , we build upon an existing technique for solving the Boltzmann transport equation which combines analytical techniques involving Legendre polynomials with numerical methods. This new method has already been used for obtaining values for the space-independent distribution function which agree with MonteCarlo simulations[4] and experiment[5]. This approach
2. ANALYTICALFORMULATION
2. I. Formulation o f transport equation
space-dependent
Boltzmann
To obtain the space-dependent m o m e n t u m distribution function ( S D M D F ) for Si, we begin with the standard Boltzmann transport equation for electrons in steady state: 1 e [Of(k, r ) ] hV~E'Vff(k'r)--hE(r)'V~f(k'r)=l_ 8t dc' (1) where k = electron wave vector, E = electron energy, r = electron position vector, f ( k , r ) = e l e c t r o n S D M D F in one of the conduction band valleys, E(r) = space-dependent electric field. 33
first
HONGCHIN LIN et al.
34
The first term on the left-hand side is often considered to give the variation i n f ( k , r) which is due to diffusion; the second term gives the effect of electric field; the term on the right-hand side is the collision integral which describes the effects of the major scattering mechanisms in Si--acoustic and intervalley phonon scattering. Therefore, we rewrite the right-hand side as:
~j~
L~ja
c L~jI
v.
(2)
In this equation the subscripts ac and iv correspond to acoustic and intervalley phonons, respectively. To solve eqn (1), while also accounting for the effects of Si's ellipsoidal energy surfaces, we have to first transform eqn (1) to the starred momentum space (k* space), where the constant energy surfaces are spherical[7]. To facilitate our space-dependent solution, we perform the inverse transformation to the position vector r. The starred electron wave vector k* and the starred electric field vector E*, as well as the starred position vector r* have the following relationship to their original values: k* = ~ l / 2 k ,
(3a)
E* = ~l/2E,
(3b)
In this well-known dispersion relation, /3 is the nonparabolicity factor which has been given the value /3 = 0 . 5 e V ~ for Si[3]. While we have used this popular expression for the band structure, it is important to note that only 7(e) will be used in our formulation. This alleviates restrictions on allowable dispersion relations to continuous analytical functions of energy. To facilitate our formulation of eqn (4) into a solvable expression, we choose the electric field to be in the symmetrical[l 11] crystallographic direction[4]. This field direction forces the S D M D F s to be identical in each valley. Furthermore, the position vector is taken to be I-D and parallel to the electric field E*(r*). Under this condition, the problem will be 1-D in real-space, and E* and r* can be replaced by the scalars E* and x*. The next step is to express the S D M D F as the sum of the first and second Legendre polynomials multiplied by their respective coefficients[4, 8]: f(k*, r*) = f ( k * , x*) =.f0(< x*) + k'g(< x*)cos 0, (6) where
(3c)
r* = ~ 12r,
0 = angle between E*(x*) and k*, f0(e,x*) = coefficient of the first Legendre polynomial which is the symmetrical part of the SDMDF, k*g(e, x*) = coefficient of the second Legendre polynomial which is the antisymmetrical portion of the SDMDF.
where the transformation matrix fi~/2 has the following form:
x/mo/m,
~1:2 =
0
.
o
Here, m 0 is the free electron mass. After transforming the Boltzmann transport equation, we obtain: 1
e
Vk,e . V,.f(k*, r*)- ~ E*(r*)- Vk,f(k*, r*) _
L
at
r*)j.I
(4)
We use the following dispersion relation to account for Si's nonparabolic band structure: h2k .2
7(E) = E +/3E ~ =
2 7x/~%~'a2Jo(e,x *) 37,(e)Xmo( ~XX,7
(5)
2m0
To ascertain the S D M D F we must now determine the two unknown functionsf0(e, x*) and g(< x*). By substituting eqn (6) into eqn (4), we can obtain a tractable form of the Boltzmann equation, which contains both symmetrical and anti-symmetrical parts. By equating symmetrical parts to symmetrical parts, and anti-symmetrical parts to anti-symmetrical parts, and performing algebraic manipulations, two equations are obtained (details in Appendix A). These two equations can be solved for the two unknown functions. The equation for the symmetrical part of the S D M D F .1~(~.x*) is:
•2Jo(e,x*)
2eE*(x*)
I , #fo(E, x*) x eE*(x ) ~
&ax~
t-[eE*(x*)]
?~,2Jo(e,x*)] 2 1 ~ ~+~ 1
efo(,, x * ) ] _ 2x/~(,)ro e aE*(x*) OJa(< x*) + ax* J 37'(e)3mo ax* a~
27 (_e)? "(e ! l roeE*(x*) ?,(e)2 jT,(e)2x/7(e)mo
4xf2m~m? '2D~¢ ~ > " ( e ) 2 rch4p
D .(mt mt) 1 +,,=, x/2=h3peo,, exp(hog.[-KoU)- 1 {x/Y(e + ho),,)7'(e + hog.)[exp(ha~,,/KoT)fo(e + hoJ.,x*) - f o (e, x *)] + x/Y (e - hog.)7 '(e - hoo,,)[fo(e - ha~., x *) - exp(h~o./KoT)fo (e, x *)]} = 0,
(7)
Boltzmann transport equation for Si where r 0 = proportionality factor which relates the p h o n o n scattering rate to energy, z0 Has units of eV m/2s and its values are given in Appendix A, 7'(e), y"(e) = first and second derivatives ofT(e) with respect to energy, respectively, mt, m~ = Si's transverse and longitudinal effective masses, respectively, Dac = acoustic phonon deformation potential, p = density, T = lattice temperature, K0 = Boltzmann's constant, D,, = intervalley phonon deformation potential and subscript n stands for different kinds of intervalley phonon scattering, o~,=intervalley phonon vibrational frequency. The second-order differential/difference equation given by eqn (7) represents a boundary value problem in two dimensions: real-space and energy space. The four sets of boundary conditions used to help define the problem are depicted in Fig. t and are described below: 1. The distribution function at x = 0, which corresponds to the region where electrons are injected from an ohmic contact, is given by a Maxwellian. 2. At x = / , which is the right edge of the device, there we can use either a Dirchlet condition where the S D M D F is specified, or a N e u m a n n condition with the normal, real-space derivative of the distribution function equal to zero. In the presented calculations, the Dirchlet condition is used, with the requirement that the distribution function is equal to its homogeneous-field limiting value. 3. The S D M D F is zero at very large energy, which for our calculations is approx. 1.5 eV.
35
4. Since we are assuming there is no generation or recombination, the product of electron concentration and drift velocity n(x)v(x), and therefore current density are constant. The value of current density depends on the doping concentration. We use the requirement of constant current density as a boundary condition. In addition to specifying boundary conditions, it is important to recall that we are interested in physical phenomena. However, eqn (7) is expressed in terms of starred-space variables. Since the vectors in eqn (7) have been reduced to scalars, we specify the scalar relationship between E*, x* and their real-space counterparts E, x. These relations, which are derived in Appendix A, are given by: E* =
m~E,
(8a)
x* =
--x,
(8b)
where the conductivity effective mass m c is given by:
1 _ 1 (2+1~. mc
3 \m t
(9) mj
To facilitate our formulation, we first transform E and x to E* and x*. Eventually, after eqn (7) is solved, we will then transform back to E and x.
2.2. Formulation of average electron energy and drift velocity While the S D M D F provides essentially all the information required for performing a comprehensive device analysis, submicron device performance is often evaluated in terms of the average values of carrier velocity, energy and concentration. Furthermore, the product of average electron concentration and drift velocity provides one of the boundary
Yo(¢ > N A ¢ , x ) = 0 = NA¢
l
fo(~, o) , -
......
~
fo(~, l)
¢=A¢
xj-1 xj xj+l
x----0
x= l
n(x)v(x) = g Fig. 1. The 2-D domain for the space-dependent Boltzmann transport equation is shown. The vertical dimension represents the energy domain whereas the horizontal dimension gives the real-space coordinates. At the left boundary the SDMDF is a Maxwellian; at the right boundary the SDMDF is given by its high-field limiting value; at the top the SDMDF is equal to zero; along the bottom the boundary condition is given by a constant value of current density at each real-space coordinate.
36
HONGCHIN L1N et al.
conditions for the numerical solution. Thus, the formulas used for calculating these average quantities from the SDMDF are given in this section. The carrier concentration in terms of the SDMDF is determined first: 1 ('~ n(x*) = 47r - - 3 J0 fo(& x*)h(e) d& (1O) where h(e) is the density of states which accounts for Si's ellipsoidal energy surfaces and nonparabolic band structure: .4f 2 2n
,
.
~9h3n
x/~m,~mT
× EeE,(x,) ~j;(< x*)
n(x*)vT(x*) = nAx*)vXx*) + n~(x*)vy(x*) + nAx*)v~(x *) =
7'(E)2
ef0(
f~ ?(e} [ t?fo(e,x*) ~?fo(&x*)~ [_eE*(x*) ~ ]ae. (15) J0 7 ' ( e ) 2 0e The product of concentration and total drift velocity for electrons in the y-valley and the z-valley are identical to eqn (15), so the average product of total carrier concentration n(x*), and average drift velocity vT(x*) for the entire conduction band is given by the vector sum of the products of concentration and total drift velocity in each valley: x
2x/-2 [m~ml 1
l, 1),0
fo ?'(e) 27(e)
(13)
In contrast to average energy, the drift velocity originates from the anti-symmetrical part of the SDMDF, and is difficult to calculate because it is a vector and depends on which valley the electron is located. Furthermore, the nonparabolic band structure and the ellipsoidal surfaces of constant energy give rise to additional complications. After performing a great deal of mathematical manipulations (details in Appendix B), the x-component of the product of concentration and drift velocity in the x-valley is: n~(x*)v~(x*) =
~mo(m~~'1 my 1_, mjl-%
(11)
It is worth noting that the anti-symmetrical part of the SDMDF does not contribute to the integral for the carrier concentration. By definition, the space-dependent average electron energy is given by: f0~ 1 ef0 (& x*)h (e) de. (12) e(x*) 4n3n(x,) As with carrier concentration, the anti-symmetrical part of the SDMDF does not contribute to the integral for the average energy, To obtain an expression for the product of average velocity and electron concentration, we begin with the definition:
f v(k)f(k, x*) dk.
x/3rn0m~j0 7,(e)2
[c'}f0(e,x*) {~f0(e,x*)-] [ (14c) x LeE*(x*) & ~xxg j d e , To obtain the product of concentration nx(X*) and total drift velocity v'(x*), for electrons in the x-valley at position x*, we take the vector sum of the three components given by the above three equations: n~(x,)v%~c,) = 2 9h3n ~ 2
.
h({) = -Y~- (m7ml)l'2?(e)127'(e).
n(x*)v(x*) = ~ 1
n~(x*)v~(x*)=T 9h3g
(14a)
where
× I eE*(x*) {~fo(e,& x*)
~fo(e, x*) 1 ~ x ; .jd&
(16)
The above equation is in the [111] direction as expected. After taking the magnitude of n (x*)v T(x*), we have: =
2,/5
x eE*(x*)
5fo(& x,) &
{~fo(&x*) 1 ~xx* J d&
(17)
Of course, by dividing eqn (17) by n(x*), which is given by eqn (10), values for average electron drift velocity over all conduction band minima are obtained. 3. NUMERICAL SOLUTION
nx (x*) = concentration of electrons in the x-valley at the position x*, v~(x*) = velocity of electrons in the x-valley for the x-component at the position x*. By using similar analyses, the y- and the z-components of the products of concentration and velocity in the x-valley are:
×
I- _,,t x ..) ~£(e, x*)
ef0(e,x*)7jae,
(148)
With the analytical expression of the Boltzmann transport equation derived and given by eqn (7), the next step is to apply numerical techniques to obtain the solution, which is the symmetrical part of the distribution function f0(& x). The basic strategy we use is to approximate eqn (7) by a system of coupled difference equations, and then solve that system. The first step is to discretize eqn (7) with respect to energy E and position x*. We choose mesh points in the energy domain which are equally spaced with increments of Ae. In the real-space domain we employ a fine mesh in regions where the
Boltzmann transport equation for Si S D M D F is rapidly varying, and a coarse mesh when its spatial variation is reduced. To accomplish the discretization we first divide the energy domain into N points, and the position domain into M intervals with M + 1 points, which results in a 2-D mesh with N × (M + 11 points. Using the central-difference method, we obtain a difference equation for each point or node. The i, jth difference equation is listed below: 2xf~i Zo [ 2(f,j ~ l_--f,;j)
32;3molsj+l(sj+sj+1)
2(f,j--fi,j-~] sj(sj+sj+ 1
37
for f0 (E, x*). The procedure is repeated iteratively for all spatial points until the system converges. Numerically, this involves solving an N x N matrix equation repeatedly. It is interesting to note that this solution technique is facilitated by the physics of the problem. In the energy domain, the various difference equations are strongly coupled over a large energy range by intervalley phonon scattering. With this
2x~i%eE*f+l,j+l-f+l,j
I Ifi-l,j+l-~f/
2.,/-~zo(eE*\ 2 2 ( _27i2;"]_ zoeE* f,j+~--f.j , +~--~-E ) ( f + , , j - 2 f o + f _ , j ) - - ~ 1 2;2 ]x/~2;2mo s,-+sj+-~l -27ieE*+,--eE*_ l 2 X - 32;
x
Sj + Sj+ 1
I+-~-.)]U+L~7;+
~
27i2;"~.-..2]f+,J-f-,,] 1 -- -'~-i2)~dff-,j ) J -2--~E
l+~T,}.-o
j
~
1,j-I
sg+sj+ 1
32;3m0 AE
-~
ro N~/i2 ;2 mo
4x/2rntm, Dac ,G2; ~'14p
+ 27 ; o
A¢2
D2"(m2'rn1)'/25 ~,/T-7 , Fexp(h°~./GT)f+, J f,;] + / ~
.
r f - , . . j - exp(ha~./Ko T)f.f]]
v7'-'.7-'° L-
j;=0,
where f j is equivalent to f0(iAE, x*) with i and j corresponding to energy and position domains; sj and sj+l are equal to x * - x * l and x * + l - x * , respectively. We also assume x* is always less than x*+~ for any j ; AE = very small energy increment step; / . = c l o s e s t integer to ho)./Ae; 7~=2(iAE); 2t+t, = 2(i A• + hco,); Y~= 2'(i Ae) = first derivative of the dispersion relation with respect to energy; 71' = 7"(i Ae) = second derivative of the dispersion relation with respect to energy; E~=E*(x*)= electric field at x*. There is an analogous difference equation for each of the N x M mesh points in the 2-D energy/ real-space domain. Since the energy domain and the real-space domain are divided into approx. 300 intervals and 60 intervals, respectively, this represents a matrix equation with approx. 18,000 rows. A direct solution of this system proved to be computationally inefficient. For this reason, we employ two different techniques simultaneously to obtain a solution. In the real-space domain we use an iterative method (method of lines), and in the energy domain we solve our system of equations directly using a sparse-matrix method which employs the minimum degree algorithm[9]. The method can be understood by first observing that in eqn (18), the distribution function at the point x* can be determined if the distribution functions at the points x * j and x*+l are already known. In other words, if we know f0 (e, x*_ l) and f0 (E, x*+l ), we can solve for f0 (E, x*). Of course we do not know f0 (E, x*+ 1), so we assume its values and then solve forf0(e, x*). The result will be improved values
(18)
strong coupling, a direct solution of the matrix equation is prudent. However, in the real-space domain, the difference equations are coupled only to their nearest neighbors, and rapid convergence is obtained using the iterative method described. Before we provide details of the numerical methods, however, we first implement boundary condition (4), which is described above, into the system of coupled equations. To implement constant-current density as a boundary condition, we replace the first difference equation in our N x N matrix, which is i = 1, by the product of concentration and average drift velocity given by eqn (17). The discretized form of eqn (17) is: 2,~z0
~ 7 ~
J = n(X*)VT(X*) = 3h3~ 2 X / ~ , = , x (eE, f+ L~-
t,J
~ 7;
f~0+-~-- -- "~ J - 1/ ~a E -, Sj"~- Sj + I ]
(19)
where J is a constant for any x*. One may now obtain a better understanding of the numerical techniques employed by considering the following discussion. For a given point x* in the real-space domain, we have N - 1 equations in the energy domain in the form of eqn (18), and one equation in the form of eqn (19), The equations analogous to eqn (18) are enumerated from i = 2 to i = N, while we place eqn (19) in the first matrix row, which corresponds to the energy point i = 1. This set
38
HONGCHINLIN et al.
of equations can be expressed by the following matrix equation:
--~j
I~_I--~,+1F~j+l
I,
.002
I
~
'
I
'
'
'
I
'
'
'
I
'
'
'
i
(20) .01,urn
.001
L,
where
.o0o
If.l% ~f = column matrix
t ¢ O.05,um
] f.i~j
o
0
2
.4
.6
.8
Energy(eV)
F:!
l~Ij, [/
1=
column matrix on the left point of x* with kth iteration, ~ $ ~ = column matrix on the right point of x* with k - lth iteration, 1, Rj ~ J = coefficient square matrices of pjk P~ l and Pf+l:, respectively.
All the values on the right-hand side o f e q n (20) are presumed to be known, so P~ can be obtained using the sparse matrix solution algorithm[9]. Once r~ is obtained, we use it to obtain an improved value for the distribution function at the next point in the space domain P~+ i. We perform this calculation for each real-space point until we reach the end of the realspace domain. The process is then repeated until the system converges• To make this system converge faster, we use the SOR method[10] with the expression given below:
=.,
0 .
-[i ,F)
,-
~i+iP~+;
,
(21)
6 where w is the relaxation parameter. We found that the system converges most rapidly when w takes on values of around 1.4. This technique gives results which agree with Monte-Carlo calculations, and require considerably less CPU time to calculate.
Fig. 2. Electron energy distribution function for electric field value of 100kVcm -t at various positions marked by arrows. S D M D F to have a carrier concentration of 10 ~vcm 3 in the homogeneous-field limit. Figure 2 shows the distribution function at different positions for an electron ensemble which is in the presence of a constant electric field with a magnitude of 100 k V c m -1. It is interesting to note that the distribution spreads over a wider energy range as the electron ensemble drifts toward the far electrode. Eventually, the distribution function reaches its non-Maxwellian homogeneous field limit. In Fig. 3, we replot the distribution function, which is given in Fig. 2, on a logarithmic scale. Figure 4 shows the product of the distribution and the density of states, f0(E, x ) h ( O , at different points throughout the device, for a field magnitude of 100 kV cm 1. From the figure it is clear that electrons obtain more and more energy as x increases, so the peak of this 3-D plot is moving from low energy at small distances to higher energy at larger distances. In Fig. 5 we give the resutls of average velocity and carrier concentration investigations. Values for the electron drift velocity (solid curve) and the carrier concentration (dashed curve) are plotted as a function of position for an electric field of 100 k V c m -1. As expected, the product of velocity and concentration is constant, which satisfies the continuity equation ~ - [n(x ) v ( x )] = O. CX
4. R E S U L T S
To analyze the accuracy of our new method for determining the S D M D F , we performed a series of calculations. In these calculations our semiconductor device is an Si bar of length 0.2 gm, at the lattice temperature of 300 K with a constant applied potential. We assume that the electric field within the bar is spatially constant. Electrons are injected into the bar with an equilibrium Maxwellian distribution at x = 0, as they would be from an ohmic contact. Electrons then gain energy from the field, and the S D M D F varies along the device until it reaches its homogeneous-field limiting value. We normalize the
In addition, we compare our values for drift velocity with those calculated by Monte-Carlo simulations[11] (open circles). Agreement between the two methods is usually within 5%, and both techniques predict velocity overshoot which is characteristic of nonequilibrium electron transport. In Fig. 6 we show values for average energy which were calculated using the S D M D F , As is expected, average electron energy values (solid curves) increase gradually with distance until they reach the homogeneous-field limit. The figure shows the results of calculations for three different electric field magnitudes: 100, 60 and 20 kV cm -1. We performed similar
Boltzmann transport equation for Si 1 0 -~
' ' I ' ' ' I ' ' ' I ......
10-~
I ' ' "_---
?
~
°o
1.5
10 -5
1 0 -e
o~ ~ "~
1.5
-~°°Oooooooooooo
~e C ]9,urr
x
\\
1 0 -~
39
...----
10-8
/"
i
I 0 -g
o o
1 0 -10
.5, co r.)
>
0
,,~1 .... 0
i0 -~
~ ~ ~ ] ~ 0 .2
~ I t 4 6 Energy(eV)
.8
I.E
Fig. 3. Electron energy distribution function on a logarithmic scale for electric field value of 100 kV cm t at various positions marked by arrows.
calculations using the M o n t e - C a r l o p r o g r a m which are plotted as open circles. The figure indicates t h a t the two m e t h o d s are in excellent agreement.
5. C O N C L U S I O N
W e have solved the space-dependent B o l t z m a n n t r a n s p o r t e q u a t i o n to o b t a i n the space-dependent m o m e n t u m distribution function. The m e t h o d , which expands u p o n an earlier space-independent technique[4], uses Legendre polynomials to formulate a space-dependent B o l t z m a n n equation. T h e resulting partial differential/difference form o f the B o l t z m a n n e q u a t i o n is then solved numerically. The numerical
.05
I .... .1
I,,, .15
0 .2
Distance(#m) Fig. 5. The space-dependent average electron drift velocity calculated by the SDMDF (--) is compared with electron drift velocity calculated by Monte-Carlo method (O) for the electric field value of 100 kVcm-L The dashed line is the space-dependent concentration obtained from the SDMDF. technique uses a novel c o m b i n a t i o n of a direct sparse matrix a l g o r i t h m a n d a n iterative technique. In addition to space dependence, the effects o f inelastic p h o n o n scattering a n d n o n p a r a b o l i c b a n d structure are included in o u r formulation. W i t h this a p p r o a c h , the S D M D F s are calculated at least 10 times faster t h a n those o b t a i n e d from M o n t e - C a r l o calculations. W i t h accurate S D M D F s determined, space-dependent values for average electron energy, electron drift velocity a n d carrier c o n c e n t r a t i o n s were o b t a i n e d which agree well with those o b t a i n e d from M o n t e Carlo simulations. In addition to calculating these average quantities, the technique gives the entire distribution function which is i m p o r t a n t for modeling
Fig. 4. Electron distribution function vs energy and distance forms this 3-D picture for the electric field value of 100 kV cm-t. The distribution function has been multiplied by the density of states.
HONGCHIN LIN et al.
40
real-space variable from r* to the I-D position variable x*. The direction of x*, with the unit vector denoted as 2", is parallel to the transformed electric field E*(x*). As was mentioned in Section 2, we divide our solution into symmetrical and anti-symmetrical components. In order to obtain a closed set of equations in this analysis, we decompose k* into k* and k~ which are the wave vectors perpendicular and parallel to E*(x*), respectively. By using the dispersion relation given by eqn (5) we then make the following approximation which has been justified in the space-independent case by comparison with Monte-Carlo calculation[4]:
100kV/cm
60kV/cm E2
.1 0-~-'
V/cm ~ ~_i ~ ~ ~ 1 ~ 05
7(E) = h2k~Z- . 3 2rn0
; ~ ~ i ~ t ~ tiff
.1
.15
.2
Distance(#m) Fig. 6. The average electron energy calculated by the S D M D F (--) is compared with average energy calculated by Monte-Carlo simulations (O) for different electric field values of 20, 60 and 100 kV cm L
(A1)
From the dispersion relation, the derivative of energy with respect to k* can be expressed as another form which is: h2k * Vk. E 7 '(E )m 0 (A2) With the above equation, as well as the chain rule, eqn (4) can be written as an expression in 4-D phase-space:
M O S F E T reliability, E P R O M p r o g r a m m i n g a n d H E M T p e r f o r m a n c e . In s u m m a r y , we h a v e p r e s e n t e d a n a c c u r a t e m e t h o d for solving the B o l t z m a n n t r a n s p o r t e q u a t i o n w h i c h is sufficiently e c o n o m i c a l for daily use in device simulation.
Acknowledgements--The authors are grateful to ShiuhLuen Wang and Ken Hennacy for helpful discussions. This work was supported in part by the Semiconductor Research Corporation under Contract No. 91-DJ-130, and the National Science Foundation under Contract No. ECS9010457. REFERENCES
1. J. Frey and N. Goldsman, IEEE Electron Device Lett. EDL-6, 28 (1985). 2. M. Fukuma and W. W. Lui, IEEE Electron Device Len. EDL-8, 214 (1987). 3. C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983). 4. N. Goldsman, L. Henrickson and J. Frey, Solid-St. Electron. 34, 389 (1991). 5. N. Goldsman, Y. J. Wu and J. Frey, J. appl. Phys. 68, 1075 (1990). 6. H. Lin and N. Goldsman, Solid-St. Electron. 34, 1035 (1991). 7. C. Herring and E. Vogt, Phys. Rev. 101, 944 (1956). 8. E. M. Conwell and M. O. Vassel, Phys. Rev. 166, 797 (1968). 9. A. George, J. Liu and E. Ng, SPARSPAK User Guide: Waterloo Sparse Linear Equation Package. Department of Computer Science, University of Waterloo, Ontario, Canada (1980). 10. W. F. Ames, Numerical Methods for Partial Differential Equations, 2nd Edn. Academic Press, New York (1977). 1l. N. Goldsman and J. Frey, IEEE Trans. Electron Devices 35, 1524 (1988). APPENDIX A
Derivation of Boltzmann Transport Equation To derive the space-dependent Boltzmann equation [eqn (7)], we start from eqn (4) which is the Boltzmann equation transformed into the starred phase-space, which has six dimensions: three in momentum space, and three in real space. A 6-D analysis is indeed cumbersome, so we limit our real-space analysis to one dimension, and we change our
hk* "2* 0 f ( k * , x * ) 7'(E)m 0 ~x*"
x
eE*(x*)Y:*
hk* ~ , [-~f(k*,x*) 1.~t y,(,)m0~J(k*,x*j=| J
. (A3)
L
To facilitate our formulation, we limit our specification of momentum space to directions perpendicular and parallel to the electric field. We then use Legendre polynomials to express the distribution function as the sum of its symmetrical and anti-symmetrical components, and then substitute this expression, which is eqn (6) with k * = k * cos 0, into eqn (A3). After some mathematical manipulations, the following expression is obtained:
hk* [afo(e,x*) +k,~g(E,x*)l y '(E )mo L- Ox *
" ?~X
*--J
eE*(x*) [ h2k* cJo!~: x*) h LT'(Qmo & +g(~, x*) h 2k~':
+ 't'(~o
Dg(E, x*)l = [g.fo(<, x*) ] k ' g ( , , x*) <~ J L at j~ ,({)
(A4)
where the collision term on the right-hand side of eqn (A4) has been separated into its symmetrical and anti-symmetrical components; k~ and E*(x*) are the magnitudes of k* and E*(x*), respectively; and l/r (E) is the summation of all acoustic and intervalley scattering rates. In Ref. [4], r(E) is shown to be well approximated by: r(Q =
-
TO
~ ,
(A5)
;"(EICg(E) where z0 represents all the coupling factors for deformation potential scattering. The approximate values of l/z 0 are the slopes o f the curve obtained by plotting total scattering rate vs y'(E),~314,51: r 0 = 4 . 2 x 1 0 - t 4 s . e V t'2, ( <0.055eV, % = 2 . 1 x 10 las.eWl'2.
E >0.055eV.
By equating the terms of the same Legendre polynomial order, and using eqn (A5), we eliminate the independent variable k~, and obtain the following two equations from which to determine the unknown functions fo(E,x*) and g(E, x*):
~'(Qm0 L ax*
eE*(x*)
0T:, ]
Y'(Q 7 ~ g ( E , x * ) , %0
(A6a)
Boltzmann transport equation for Si 1 ~'27(e)
Any vector x*2* can be transformed to the original space by eqn (3b). The result is:
Og(E,x*) eE*(x*)[g(e,x*) .
1 1%
(A6b>
+ 3y'(e) & J J L Ot Jc To eliminate g(E, x*) from the above two equations, we substitute eqn (A6a) into eqn (A6b) and obtain an equation which only contains the symmetrical part of the distribution function: 27 x / ~ % _a fOfo(e,x*) 3y'(Q3m0 Ox* L Ox*
eE* x* a [
+ 27(Q%
41
'
which is not in the [111] direction and depends on which valley the electrons are located; for example, if the
eE*(x*)~] I
eE*(x*)z o f. of° (E' x*) +7'(Qzx/y(Qm0l-- ax*
\ m x'my'm~/
~?fo(,,x*) ee*(x*) Of0(,,x*)]
eE*(x*)
= dE
J
(A7) I
ot
_1~
Upon taking the derivatives of eqn (A7), the final form is: 3y,(e)3m° [ +~[1
2eE*(x*) & O~xT - + [eE.(x.)]2 O
~XO2
27(,,7"(e!]
%eE*(x*> ~Ofo(,,x* , e E ' ( x ' , ~ ]
~'(E)= JT,(,)~moL COX*
+ 3Y'(E)3m0
&
Ox*
=foeo¢,,x*q L
at
J¢"
(A8)
The acoustic and the intervalley phonon collision terms have already been given in a previous work[4], and are rewritten below to facilitate our presentation: c~t
nh4p
Jac
~y'(e)
2
1+
f0(E, x*)
+f 7(')+(1
L~
+Y(')Y"(')XIK T] cof°(''x*) 7(QK°TCO2f°(e'x*)'I ~ ) o ~ -~ ~ ~ -j (A9)
and
Ofo(,,x*) at
_ ~ D2.(rnZtm,)',2f
1 .,/2~h3p~o. L e x p ( ~ ) - I J
3,-.-%
+ he).) [exp(R.)f0 (E + h~o., x*)
-- fo(E, x*)]
+ x / ~ - he).)y '(E - hco.)[fo(E -- he)., x*) - exp(R.)f0(E , x*)]}, (A10) where R. = ho).lKo T. Finally, eqn (7) is obtained by taking the summation of these scattering mechanisms, and moving the left-hand side of eqn (A8) to the right-hand side. After eqn (7) is solved numerically, the distribution function, which is in the starred space, has to be transformed to the original space. First of all, E(x*) can be expressed as (I/x/3)(1, 1, 1)E(x*) since the electric field in the [Ill] crystallographic direction, so by using eqn (3b), the relation between E*(x*) and E(x*) is: E*(x*) =
,
-~f]E(x*)
,m/
electrons are residing in the x-valley, the position direction will be: ' /T/t '
However, electron transport in Si is governed by the combined effect of motion within tile six valleys in momentum space. By taking the vector average of the x-direction in each valley, the value becomes: m.,/~om~{ 1 2", ° "/-- +--]x*.
(A1 la)
3
\m I mt/
Thereby, the relation between x and x* is:
E*(x*) = / m° E(x*),
(Allb)
where E(x*) and E*(x*) are the magnitudes of E(x*) and E*(x*), respectively. As was discussed above, the unit vector 2* is parallel to E*(x*), thus 2" is given by: ., X
SSE 35/1--D
/-~¢f =
~-
1 -
1 ,
1
x="
m,/ vom" [r- 1- + - - l x * = 3
k,m,
mr,/
mUo
/ Vx*. ~/m c
(A13)
APPENDIX B
Derivation of Average Electron Drift Velocity To formulate the expression of average drift velocity, we begin with the definition of product of concentration and
HONGCHIN LIN el
42
therefore
drift velocity from semi-classical kinetic theory: . (x *)v(x *) =
1 4-~
fv(k)f(k,x*)dk.
(B1)
By substituting the instantaneous electron velocity v(k) = 1/h V ~ into the above equation and transforming k to the starred space, the preceding expression becomes: 1 ('~1,2
n(x*)v(x*) = 4rc~ J ~
V~.~f(k*, x*)]fi ~'21dk*,
(B2)
where I~t t21 = x/mxm,m:/m 3 is the Jacobian determinant. Upon using eqn (A2) and the distribution function expressed as the sum of the symmetrical and the antisymmetrical parts, we have: 1 ,, . I~l/2hk * n(x*)v(x*) = 47z~v' m, m~m:/m~ j ~ [fo(~,x*)
~\/'mcnvm./rn ~ I
x k* ' g(E, x*) dk*,
(B3)
where k*g(E, x*)cos 0 has been expressed as a dot product, and the vector g(E, x*) is parallel to the electric field E*(x*) in starred space. Furthermore, the integral of k*f0(E,x*) goes to zero because k* is odd and f0(E, x*) is even with respect to k*. To avoid the writing of lengthy formulas, the vector form of eqn (B3) is considered by its x-, y- or z-components separately. We write the x-component of n(x*)v(x*) as:
h x/m~.m. (" k, n ( r * ) G ( x * ) - 47r 3 m ~
)
+ k~g~.(~,x*) + k:g:(~, x*)] dk* -
(B6) &
0x*
(B4)
ing to spherical coordinates, we obtain:
2T0
CvJ,( ,
n(x*)G(x*)=~4
3m, " J T ' ( ~ ) ~
x[eE*(x*)Of°(~ x*)
gf°("x*)~k*2dk*~x*J
f
x leE*(x*) L _2x/2mtml
~/;(e, x*) al0(E,x*)l - "~ I d~ a~ ~x* J
/~
,, '(e )2x/~' (~:)m0
J'
cE
9h3x 2
F
g(~, x*) = g(~, x*)2*
....
×[_e~-tx*;
0S0(,,x*q ~x*
J'
n(x*)
~ .
(B9)
where G(x*), n~.(x*) and n:(x*) are the electron concentration in x-, y- and z-valleys, respectively. Upon modifying eqn (B8) with G(x*)=n(x*)/3, the x-component of the product of concentration and velocity in the x-valley is:
Since g(~, x*) is equal to 2* • [g(~, x*)2*], the vector g@, x*) can be expressed as:
&
7(0 ./,(~)2
?,/;(~ x*)l :"~'.', / d e ax ]
(BIO)
Using a similar approach gives the /-component of the product of concentration and velocity in the ./-valley:
C~fo(e,x*)l ax*
S
r0
4 3morn'
L
•
. (B7)
~ f~ 7(E_!. 3h3x 2 ~/ m, ~/3moJo ?"(e) 2
F ?Jo(,E x*) ×/eE*(x*) ~i
hz0
[eE*(x*)ifo(~x*)
J"
Next we substitute gx(E,x*) into eqn (B4), then using k 2 = 2moy(E)/3h2, which is identical to eqn (A1), and chang-
*
h x/Im~mz P b.~, |",'g4,,x*)dk*, 4 rr~ m2 j 7 (E)
x
tx ;
G(x*)=nr(x* )=n:(x )=
where k,, k,. and k: are the components of k*; the integrals of k~k~g~(Eix*) and k~k:g:(E,x*) have vanished due to the orthogonal properties of k~ and ky, as well as k~ and k~. Next, the appropriate expression of gx(E, x*) has to be determined. By recalling eqn (A6a), we can write g(E, x*) in terms of the symmetrical part of the distribution function:
g(~, x*)
~ / e ~ z
q3m~L
0
Since we chose the electric field to be in the [111] crystallographic direction, the distribution functions are identical for each valley. Hence:
*
"iT'(,) [k~g4''x
x
hz o
'(E )2y ~ m
n(x*)v~(x*)= 2x/F2z° ~
Jy'(Om0
-
gx(E,x*) =
By making use of the dispersion relation, the integration with respect to k* can be changed to the following integration over energy ~:
+ k* . g@, x*)] dk* , . ~"~l'2hk* 4u
al.
(B5)
~[ 3mom, Jo )"(E) 2 0J'°(e' x *)
?-
~/;)(" x * ) l d,
~x*
]
(BI1)
where i = x , y or z a n d j = x , y or z. Finally, the average drift velocity can be obtained from this formula by dividing by the electron concentration as described in Section 2.2.