Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://lizard.phys.msu.su/home/science/Masyukov-Dmitriev-11-JMathSci-InN.pdf Äàòà èçìåíåíèÿ: Thu Feb 10 10:31:55 2011 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:07:08 2012 Êîäèðîâêà:

Journal of Mathematical Sciences, Vol. 172, No. 6, 2011

A NEW NUMERICAL METHOD FOR THE SOLUTION OF THE BOLTZMANN EQUATION IN THE SEMICONDUCTOR NONLINEAR ELECTRON TRANSPORT PROBLEM N. A. Masyukov and A. V. Dmitriev UDC 538.935+517.958+537.84

Abstract. A new iterative method for solving the Boltzmann transport equation in the space-uniform case is introduced. The method is based on the use of a mesh in the momentum space to represent the distribution function. In intermediate points the function is found with the help of interpolation. The method is used to study the hot-electron transport in bulk n-InN, which is a promising material for optoelectronic applications.

1. Intro duction The Boltzmann transport equation (BTE) was initially developed to study gas kinetics [6, 7, 10]. It is also used for studying charge and phonon transport in solids in the semiclassical approximation [1, 2, 5, 17, 26]. A few analytical methods for solving the BTE were developed. However, they can be applied only in simple cases, where all relevant scattering mechanisms can be treated as elastic processes and the BTE can be linearized with respect to the charge-carrier distribution function. This approach is called the -approximation. It lies at the base of the standard theory of galvanomagnetic and thermoelectric effects in solids [1, 2, 17, 26]. Problems where one has to deal with the nonlinear BTE and nonelastic scattering processes are much more complicated. In some particular cases, a set of algebraic balance equations could be solved instead of the nonlinear integro-differential BTE [3]. But generally the BTE can be solved only numerically. In this paper, we study the problem of charge transport in semiconductors in strong electric fields. This problem is called the hot-electron problem or the problem of nonlinear charge transport, because the mean energy of electrons in the field usually exceeds the lattice temperature and the current-voltage characteristics become nonlinear. When one is looking for a suitable approach to the solution of the hot-electron transport problem, several factors should be taken into account. The band structure of the material and the relevant scattering mechanisms are the most important of them. Thus, the hot-electron transport problem should be considered separately for any particular substance or for a group of substances with similar properties. Here we consider the wurtzite indium nitride InN and its alloys with GaN and AlN. These semiconductors were synthesized not long ago and are very promising for optoelectronic applications, including their use in general light sources [14]. So the problem of charge transport in these semiconductors in high electric fields is important. It is worth noting that the method we have developed is not limited to these materials but can also be used for hot-electron transport calculations in other semiconductors. In the space-uniform case, the BTE for electrons has the form eE f (k) f (k) =- + St(f ), t k with the normalizing condition f (k) = NV .
k

(1)

(2)

Translated from Fundamentalnaya i Prikladnaya Matematika, Vol. 15, No. 6, pp. 77­97, 2009. 1072­3374/11/1726­0811 c 2011 Springer Science+Business Media, Inc. 811


Here E is the electric field, f (k) is the electron distribution function, k is the electron wave vector, e denotes the electron charge, V is the crystal volume, St(f ) =
k

{W (k k)f (k )[1 - f (k)] - W (k k )f (k)[1 - f (k )]}

(3)

is the collision integral, and W (k k ) is the electron transition rate from the state k into k . The relevant scattering mechanisms in nitride semiconductors do not affect the spin of electrons. Thus, there is no summation with respect to spin in the collision integral. The well-established Monte Carlo (MC) approach [15] is a powerful tool for solving the BTE. It makes it possible to study charge transport in a wide range of materials and physical situations with rather high accuracy. But in the case of nonlinear BTE, the method is very time-consuming from the computational point of view. The problems that appear when one tries to study electron transport in the transient regime seem to be its second significant disadvantage. Deterministic methods [9, 12] based on the collision-integral linearization with respect to the distribution function assuming f (1 - f ) f and on the distribution-function expansion in terms of a suitable basis function set have been developed in recent years. Thus, a linear equation set is to be solved instead the BTE. These methods are quite accurate and effective from the computational point of view. But they become inappropriate when the linearization cannot be carried out, e.g., when the electron gas is degenerate. It is worth noting that in InN the free electron concentration is usually rather high, so that the electron gas is degenerate even at room temperature. The hot-electron transport in the bulk n-InN was numerically studied in [4, 11, 18­21, 23] using the MC approach and in [22] using the quasi-equilibrium thermodynamic approach. In all these publications the free electron concentration was assumed to be about N = 1017 cm-3 . However, no comparison with experimental data was presented in these papers. Indeed, such a comparison could not be performed because the technology at that time allowed one to grow InN crystals with the free electron concentration not less than N 1018 cm-3 [8]. Respectively, in an experimental study [25], samples with a free electron concentration of the order of N = 9 â 1018 cm-3 were used, and the lattice temperature was T = 77 K. Under such circumstances, the electron gas is significantly degenerate; thus, the deterministic methods mentioned above become inapplicable, and the MC method becomes very time-consuming. The iterative method that we develop in this paper requires only the smoothness of the distribution function. Using this method, one can find the stationary distribution function as well as its time evolution. The hot-electron transport in the bulk n-InN was studied, as an important example, using this method. The results we obtained numerically agree well with the available experimental data [25]. 2. Physical Mo del We consider the bulk n-InN crystals lattice temperature T = 77 K. The doping concentration. Therefore, we neglect holes equilibrium. Thus, the lattice temperature with free electron concentration N = 9 â 1018 cm-3 at the concentration Ni is assumed to be equal to the free electron in the valence band. The phonon system is assumed to be in remains constant.

2.1. Material Parameters. It is worth noting that at the present time there is a certain disagreement about the values of InN material parameters. In the simulation, we will use the parameters recommended in [24] that are presented in Table 1. A two-valley model is used to approximate the band structure of InN. The lower conduction band minimum is assumed to be at the -point of the Brillouin zone and the energy dispersion relation in the neighborhood of the -point is assumed to be isotropic:
2 k2

2m 812

= (1 + );

(4)


Table 1. Material parameters of InN Quantity Symb ol Unit Value effective mass in the -valley m m0 0.12 effective mass in the A-valley mA m0 1.0 band gap g eV 2.0 energy separation between -valley and A-valley minima - A eV 2.2 LO-phonon energy LO meV 73 TO-phonon energy TO meV 57 HF dielectric constant 8.4 static dielectric constant 0 15.3 5 cm/s longitudinal sound velocity sLA 10 6.24 5 cm/s transversal sound velocity sTA 10 2.55 acoustic deformation potential DA eV 7.1 optical deformation potential DO eV/cm 109 piezoelectric modulus e14 C/m2 0.375 3 mass density g/cm 6.81 here is the nonparabolicity parameter. It is defined as = 1 g 1- m m0
2

,

(5)

where m0 is the free electron mass. The next nearest energy minimum in the conductivity band is assumed to be at the A-point. The dispersion relation of the electron energy in the neighborhood of the A-point is taken to be isotropic and parabolic: 2 k2 , (6) = 2mA where k is counted from the A-point and from the A-valley minimum. 2.2. Scattering Mechanisms. In this section, transition rates for the electron scattering mechanisms important in InN are presented and their properties briefly discussed. Ionized impurity (i), deformation (DA) and polarization (PA) acoustical and polarization (PO) optical phonons are taken into account in our model as intravalley scattering mechanisms in the -valley; i-, DA-, PA-, DO-, and PO-scattering processes are intravalley scattering mechanisms in the A-valley; DO-scattering is assumed to be responsible for intervalley transitions. The following transition rate was used to describe the ionized impurity scattering: Wi (k k ) = 2 N V
i

4e2 0

2

(k, k )S (|k - k |) (k) - (k ) . |k - k | 4

(7)

The overlap integral (k, k ) of the Bloch parts of the electron wave functions appears in the case of nonparabolic dispersion relations. For normal processes it is (k, k ) =
cell

d3 ru (r)uk (r) , k

2

813


where u (r) and u (r) are the initial and final electron wave function Bloch amplitudes, respectively. The k k probability of the umclapp processes is significantly smaller, so they are not considered. For intravalley transitions in the -valley with the nonparabolic dispersion relation (4) the overlap integral equals (k, k ) = [1 + (k )][1 + (k )] + (k )(k )
kk kk 2

[1 + 2(k )][1 + 2(k )]

,

(8)

whereas for the intravalley transitions in the parabolic A-valley and for the intervalley transitions it is assumed to be unity. The term S (|k - k |) appears when one takes into account the screening of the impurity Coulomb potential 2 2 q 2 , S (q ) = 1+ 2 q 2 where the screening length can be found with the help of the expression 1 = 4e2 0 f (k) . (k) (9)

k

One can see that depends on the electron distribution function so that the BTE is nonlinear even in the case of nondegenerate electron gas statistics. Formula (9) reduces to the Debye and Thomas­Fermi screening lengths in the cases of nondegenerate and degenerate statistics, respectively. For the phonon scattering the following generic form of the transition rates is used: W± (k k ) = 2 1 (k, k )S (|k - k |) V â
j

Bj (|k - k |) Nj (k - k )+

11 ± (k) - (k ) j (k - k ) 22

,

"+" corresponds to the scattering processes with emission of a phonon and "-" to the processes with its absorption; summation with respect to j runs over all phonon branches. The dispersion relations of longitudinal and transversal acoustic phonons are assumed to be linear:
LA

(q) = sLA q, (q) = ,



TA

(q) = sTA q, .

Optical phonons are assumed to be dispersionless:
LO LO

TO (q) =

TO

The number of phonons Nj (q) = Nj (q, T ) with the wave vector q is given by the Bose­Einstein distribution with the lattice temperature. For the different types of phonons, B (q ) should be determined from the following formulas: BDA (q ) = BPA BDO BPO where D
PZ 2 DA q, 2s 2 e2 DPZ , (q ) = 2sq 2 DO , (q ) = 2 2e2 (q ) = , q2 ¯

= 4e14 /0 is the piezoelectric constant, 1 1 1 = - . ¯ 0

814


Note that only longitudinal phonons take part in DA- and PO-scattering, whereas PA- and DO-scattering are caused by both longitudinal and transversal phonons. Usually AP-scattering can be treated as an elastic process. In this case, the transition rate becomes W
AP

(k k ) =

2 1 (k, k )S (|k - k |) (k) - (k ) V

Bj (|k - k |)[2Nj (k - k )+1],
j

j runs over all branches of acoustic phonons. 3. The Numerical Metho d In a semiconductor with a multivalley band structure each band is described by its own Boltzmann equation, so a coupled system of Boltzmann equations should be solved: fj (k) eE fj (k) =- + St(fj )+ t k
l

StIV (fj ,fi ),
i=j

j = 1,...,l,

(10)

where l is the number of valleys and fj (k) is the electron distribution function in the j th valley. The normalizing condition becomes
l

fi (k) = NV .
i=0 k

(11)

In the case of InN with its - and A-valleys, we should take l = 2. 3.1. Collision Integral. The collision integral can be simplified taking into account the distribution-function axis symmetry f (k) = f (k ,k ) where k is parallel and k is perpendicular to the applied field. Let us write down the collision integral for the first valley in the form St = Sti + StOP + StAP + StIV (1) (and similarly in the second valley) and consider each term in the spherical coordinate system (z It is convenient to use the auxiliary functions k cos + k sin cos , k q 2 (k ,k ,k ,,) = k 2 + k 2 - 2kk (k ,k ,,), 1 D(k ) = 2 dk k 2 (k ) - (k ) . (k ,k ,,) =
2 The pair (k ,k ) determines the electron state k before the scattering event, k 2 = k 2 + k , and the set (k ,,) determines the electron state k after the scattering, k = |k |,

(12) k ).

kk , = (k sin , k), q (k ,k ,k ,,) = |k - k |, kk D(k ) is the density of states (DOS). For the dispersion relation (4) we have = (k ,z ), (k ,k ,,) = D(k ) = m k 2 2 1+ 2 2 k 2 . m
2

The overlap integral (8) in the new variables takes the form (k ,k ,k ,,) = [1 + (k )][1 + (k )] + (k )(k )(k ,k ,,) . 815

[1 + 2(k )][1 + 2(k )]


Now the first term in (12) can be represented by the following expression: Sti = 4e2 0
2 1

ND(k ) 2

2

d cos [f (k cos, k sin ) - f (k ,k )]
-1 0

d

(k ,k ,k ,,)S q (k ,k ,k ,,) . q 4 (k ,k ,k ,,)

So the term Sti contains double integration in finite limits. If the overlap integral is unity (intravalley transitions in the A-valley and intervalley transitions), integration over can be performed analytically:
2

d
0

S q (k ,k ,k ,,) 2k 2 - 2k k cos + -2 = 2 q 4 (k ,k ,k ,,) [(2k 2 - 2k k cos + -2 )2 - (2k k sin )2 ] 1 2
+ -

3/2

.

Further, St
OP

=

StOP,j + StOP
j

,j

,

where the summation is over all optical phonon branches,
1 + StOP,j

= D(kj,+ )
-1

d cos {f (kj,+ cos , kj,+ sin )[1 - f (k ,k )](Nj +1)

- f (k ,k )[1 - f (kj,+ cos , kj,+ sin )]Nj }
2

â
0 - StOP,j

d (k ,k ,kj,+ ,,)S q (k ,k ,kj,+ ,,) Bj q (k ,k ,kj,+ ,,) ,
1

= D(kj,- )
-1

d cos {f (kj,- cos , kj,- sin )[1 - f (k ,k )]N

j

- f (k ,k )[1 - f (kj,- cos , kj,- sin )](Nj +1)}
2

â
0

d (k ,k ,kj,- ,,)S q (k ,k ,kj,- ,,) B q (k ,k ,kj,- ,,) ,

and kj,± is the solution of the following equation: (kj,± ) = (k ) ± j . Of course, if (k ) - j < 0 for any j , then StOP,j = 0. Thus, one should calculate several double integrals over finite areas in order to evaluate collision integral terms describing OP-scattering. For AP-scattering in elastic approximation we have
1 -

Stqel

AP

D(k ) = 2
2

d cos [f (k cos , k sin ) - f (k ,k )]
-1

â
0

d (k ,k ,k ,,)S q (k ,k ,k ,,)
j

Bj q (k ,k ,k ,,) 2Nj q (k ,k ,k ,,) +1 ,

j runs over all branches of acoustical phonons. Thus, to evaluate collision integral terms describing AP-scattering one has to calculate several double integrals in finite limits. It is worth noting that if the inelasticity of AP-scattering should be taken into account, the integration in the AP-terms of the collision integral becomes more complicated. To simplify things, the following 816


quasi-elastic approximation can be introduced. will be calculated, the electron energy is much can still use the same phonon wave vector dep approximation, 1 Stqel AP = 2
j 1 + Stqel AP,j 2

In most of the k-space area where the collision integral greater than the phonon energy. Correspondingly, one endence q (k ,k ,k ,,) as in the elastic case. In this Stqel
+ AP,j

+ Stqel

-

AP,j

,

where the summation with respect to j runs over all branches of AP, = D(kj,+ )
-1

d cos
0

d S q (k ,k ,k ,,) Bj q (k ,k ,k ,,) (k ,k ,kj,+ ,,)

â f (kj,+ cos , kj,+ sin )[1 - f (k ,k )] Nj q (k ,k ,k ,,) +1 - f (k ,k )[1 - f (kj,+ cos , kj,+ sin )]Nj q (k ,k ,k ,,)
1

,

Stqel

-

2

AP,j

= D(kj,- )
-1

d cos
0

d S q (k ,k ,k ,,) Bj q (k ,k ,k ,,) (k ,k ,kj,- ,,)

â f (kj,- cos , kj,- sin )[1 - f (k ,k )]Nj q (k ,k ,k ,,) - f (k ,k )[1 - f (kj,- cos , kj,- sin )] Nj q (k ,k ,k ,,) +1 . where kj,± (k ,k ,,) is to be determined from kj,± (k ,k ,,) = k ± k
-1

sj q (k ,k ,k ,,).

One can see that the AP-scattering terms in the quasi-elastic approximation described above are analogous to the OP-scattering terms. But the first one is much more expansive from the computational point of view, because there all functions are under a double integral. Regarding the intervalley processes, one has here + - 1 2 BDO,j Stj (1) + Stj (1) , StIV (1) = 2
j

where the summation with respect to j runs over all branches of DO-phonons,
1 + Stj

(1) = D(kj,+,2 )
-1

d cos {f (kj,+,2 cos , kj,+,2 sin )[1 - f (k
,1 1

,1

,k,1 )]

â (Nj +1) - f (k St (1) = D(kj,-,2 )
- j

,k,1 )[1 - f (kj,+,2 cos , kj,+,2 sin )]Nj },
,1

d cos {f (kj,-,2 cos , kj,-,2 sin )[1 - f (k

,k,1 )]N

j

-1

- f (k

,1

,k,1 )[1 - f (kj,-,2 cos , kj,-,2 sin )][Nj +1]}, (kj,±,2 ) = 1 (k1 ) - ± j .

kj,±,2 is the solution of the following equation:
12

Of course, if 1 (k1 ) - 12 ± j < 0 for any j , the corresponding term is assumed to be equal to zero. As for the second valley, the subscripts 1 and 2 should be interchanged, and the equation to determine kj,±,2 will be (kj,±,1 ) = 2 (k2 )+12 ± j . 817


3.2. Basics of the Numerical Metho d. First, let us consider one valley (l = 1). The distribution function values at the grid nodes k i = -kmax + i +
j k = j 0 ,

1 2

0 ,

i = 0, 1,...,

2kmax - 1, 0 kmax j = 0, 1,..., , 0

are unknowns to be found. The grid can be sparse if the distribution function is expected to be smooth. The value of kmax should be chosen in correspondence to the applied field strength E so that the distribution function is well localized within the simulation area, so on the one hand kmax should be large enough. On the other hand, the computational complexity increases for large kmax , hence it should not be taken greater than required. j 0 At the initial moment of time t = 0, initial values of the unknowns fij = f 0 (k i ,k ) are taken according to the equilibrium Fermi distribution function. Then the (n + 1)th iteration step is organized in the following way. The distribution function is interpolated from the grid nodes in the (k ,k )-space using the bicubic spline interpolation. The interpolated distribution function f n (k ,k ) is used to calculate the screening length (n) = (f n ) according to formula (9) and the collision integral according to formulas from Sec. 3.1. With the interpolated distribution function, the numerical integrations in the collision integral members can be performed with as fine a step as necessary for a desired precision. Note that the collision integral values are calculated on the grid nodes only; thus we have Stij . Now the new distribution-function values in the grid nodes are evaluated according to the BTE (10):
n fij+1 = f n

k-

eE

t, k

+t Stij ,

t being the time step. Then the renormalization (11) is performed. This iterative procedure can be considered as the evolution of the distribution function according to the BTE. It is continued until the convergence is obtained. In the multivalley case (l > 1), distribution functions at l grids can be evaluated similarly to the one-valley case (l = 1) considered above. But one should note that in the multivalley case the distribution functions in different valleys are coupled via collision integral terms that correspond to the intervalley transitions. So, the iterative procedure should be interrupted only when the convergence in all valleys appears. 4. Selected Results 4.1. Fields Less Than 15 kV/cm. In the experimental research [25], the field dependence of the electron drift velocity in electric fields less than 15 kV/cm was reported at the lattice temperature T = 77 K. In this section, we will compare the results of the numerical simulation performed using our method and the experimental data from [25]. Within the specified electric-field range the one-valley model can be used because the electron energy is insufficient for the intervalley transitions. The simulation-area size was taken to be kmax = 1.2 â 107 cm-1 . One can see that the distribution function is localized within the simulation area (Fig. 1). The convergence of the method was verified in test runs. Figure 2 illustrates the time dependence of the drift velocity for different numerical-method parameters. Unfortunately, we could not find an easy way to show the convergence of the distribution function itself because of the 3-dimensional nature of the corresponding graphs. For the fixed field value E = 15 kV/cm, different time steps t and grid sizes were examined. We found t = 10-15 s and the grid size 15 â 30 to be optimal. One can see in Fig. 2 that a decrease in t and an increase in the grid size do not change the result. Furthermore, one can see that 818


1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.8 0.6 0.4 || 0.2 0 -0.2 -0.4 -0.6 -0.8 1 0.5 -1 0

k , 107 cm-1

k, 10 cm

7

-1

Fig. 1. The distribution function in the -valley, E = 15 kV/cm, T = 77 K. Grid size 15 â 30, time increment t = 10-15 s, simulation time tsim = 0.5 â 10-12 s. One can see that the field shifts the distribution function along the k axis. the results obtained with the initial approximation in the form Fn
/eq

(k ) = F (k, T , 0 ) + cos(10

-6

k)

that differ from the equilibrium distribution (n/eq) are the same as for the equilibrium initial distribution, thus being independent of the choice of the initial distribution shape.
1

0.9

0.8

0.7

vd, 10 cm s

-1 7

0.6

0.5 15â30, t =1e-15s: tCPU = 90 s 0.4 20â 40,t =1e-15s: t
CPU

= 160 s = 1030 s

15â30, t =5e-16s: tCPU=185 s 0.3 15â30, t =1e-15s, q/el: t
CPU

15â30, t =1e-15s, n/eq: tCPU= 90 s 0.2

0.1

0 0

0.5

1

1.5

2

tsim, s

2.5

3

3.5

4

4.5 x 10

5
-13

Fig. 2. Illustration of the convergence. The calculation time dependence of the drift velocity vd , E = 15 kV/cm, T = 77 K. 819


0.9 numerical result 0.8 experimantal data 0.7

0.6

vd, 107 cm s-

1 0.5 0.4 0.3 0.2 0.1 0

2

4

6

E, kV cm-

8

1

10

12

14

16

Fig. 3. The calculated field dependence of the drift velocity vd and the experimental data [25]. Grid size is 15 â 30, simulation time tsim = 0.5 â 10-12 s, time increment t = 10-15 s. The use of the quasi-elastic approximation (q/el) for AP-scattering (described above in Sec. 3.1) does not alter the result. Thus, the simpler elastic approximation can be used to save computation time. Figure 3 illustrates that the calculated field dependence of the drift velocity fits the experimental data well. We consider it as a validation of the method proposed. 4.2. Fields More Than 15 kV/cm. In this section, the results of the simulation in higher electric fields are presented. The simulation is performed as before for N = 9 â 1018 cm-3 and T = 77 K; AP-scattering is treated as an elastic process. Initially, optimal values of the numerical method parameters (time increment t, grid size (GS), time to convergence tsim , simulation area kmax , and applicability of the one-valley model) for different field ranges were found (Table 2). For this purpose the convergence was tested as in the previous section. Table 2. Numerical-method parameters used in different field ranges. E , kV/cm 15 30 50 > 50 -1 kmax , cm 1.2 1.7 2.5 4.4 (kmax ), eV 0.4 0.72 1.32 2.89 -1 kmaxA , cm 4.0 A (kmaxA )+-A , eV 2.80 grid size in the -valley 15 â 30 20 â 40 grid size in the A-valley 10 â 20 -15 s t, 10 1.0 0.5 -12 s tsim , 10 0.5 8.0 The field evolution of the integrated over k distribution function f (k ) = 2 820 dk k f (k ,k )


14

x 10

13

12

10

f(k ), cm-

8

2 ||

6

4

2

0

-4

-3

-2

-1

k||, 10 cm

0 7

-1

1

2

3

4

Fig. 4. The field evolution of the distribution function f (k ) averaged over k in the valley. The solid line corresponds to the field E = 15 kV/cm, dashed line to the field 45 kV/cm, dash-dotted line to the field 70 kV/cm. Calculations were performed using numerical-method parameters from Table 2. is presented in Fig. 4. The distribution function f (k ) is not so informative as f (k ,k ) but its changes can be more easily shown graphically. One can see how the applied field changes the shape of the distribution function f (k ).
1.01

1 2 0.99

vd, 10 cm s

-1

1.5

0.98

N/N
1 0.5 0 20 40 60 80

7

0.97

0.96

0.95

0.94

0.93

E, kV cm-

1

20

E, kV cm-

40

1

60

80

Fig. 5. Left panel: calculated field dependence of the drift velocity. Right panel: calculated field dependence of the relative electron concentration in the -valley. Calculations were performed using numerical-method parameters from Table 2. 821


Figure 5 shows the field dependence of the electron drift velocity and the electron concentration in the -valley up to fields of order 70 kV/cm. One can see that the prominent negative differential drift velocity is connected with the electron transfer to the A-valley where the carriers are heavier. Conclusions In this paper, the hot-electron transport in n-type InN with an electron concentration N = 9 â 1018 cm-3 at a lattice temperature T = 77 K was studied numerically. The calculated field dependence of the drift velocity is in very good agreement with the experimental data [25]. The results are obtained with the help of a new numerical method for solution of the Boltzmann transport equation. The physical grounds of the method are discussed and the optimal parameters of the numerical method are determined. The method could be a convenient tool for studying similar electron transport problems in semiconductors. The method can be used under different conditions, including the degenerate electron gas statistics; there are no fitting parameters and all relevant scattering mechanisms can be easily taken into account; the method can also be used to study charge transport in the transient regime. Furthermore, it is quite computationally effective. Acknowledgment. This work was supported by the Federal special-purpose program "Scientific and scientific-educational personnel of innovative Russia in the years 2009­2013" via contract P-2312. REFERENCES 1. A. A. Abrikosov, Fundamentals of the Theory of Metals, North-Holland, Amsterdam (1988). 2. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Holt, Rinehart and Winston, New York (1976). 3. F. G. Bass and Yu. G. Gurevich, Hot Electrons and Intense Electromagnetic Waves in Plasmas of Semiconductors and of Gas Discharge [in Russian], Nauka, Moscow (1975). 4. E. Bellotti et al. "Ensemble Monte Carlo study of electron transport in wurtzite InN," J. Appl. Phys., 85, No. 2, 916­923 (1999). 5. V. L. Bonch-Bruevich and S. G. Kalashnikov, Physics of Semiconductors [in Russian], Nauka, Moscow (1990). 6. C. Cercignani, Theory and Application of the Boltzmann Equation, Scottish Academic Press, Edinburgh (1975). 7. S. Chapman and T. G. Cowling, The Mathematical Theory of Nonuniform Gases, Cambridge Univ. Press (1952). 8. V. Yu. Davydov and A. A. Klochikhin, "The electron and vibrational states of InN and Inx Ga1-x N solid solutions," Fiz. Tekh. Poluprovodn., 38, No. 8, 897­936 (2004). 9. C. Ertler and F. Schruerrer, "A multicell matrix solution to the Boltzmann equation applied to the anisotropic electron transport in silicon," J. Phys. A, 36, 8759­8774 (2003). 10. J. H. Ferziger and H. G. Kaper, Mathematical Theory of Transport Processes in Gases, North-Holland, Amsterdam (1972). 11. B. E. Foutz, S. K. O'Leary, M. S. Shur, and L. F. Eastman, "Transient electron transport in wurtzite GaN, InN, and AlN," J. Appl. Phys., 85, No. 11, 7727­7734 (1999). 12. M. Galler and F. Schruerrer, "A deterministic solution method for the coupled system of transport equations for electrons and phonons in polar semiconductors," J. Phys. A, 37, 1479­1497 (2004). 13. V. F. Gantmakher and I. B. Levinson, Scattering of Current Carriers in Metals and Semiconductors, North-Holland, Amsterdam (1987). 14. C. J. Humphreys, "Solid-state lighting," MRS Bul l., 33, 459­470 (2008). 15. C. Jacoboni and L. Reggiani, "The Monte Carlo method for the solution of charge transport in semiconductors with applications to covalent materials," Rev. Modern Phys., 55, No. 3, 645­698 (1983). 822


16. R. Keys, "Cubic convolution interpolation for digital image processing," IEEE Trans. Acoust. Speech Signal Process, 29, 1153 (1981). 17. C. Kittel, Introduction to Solid State Physics, Wiley, New York (2004). 18. S. K. O'Leary, B. E. Foutz, M. S. Shur, U. V. Bhapkar, and L. F. Eastman, "Electron transport in wurtzite indium nitride," J. Appl. Phys., 83, No. 2, 826­829 (1998). 19. S. K. O'Leary, B. E. Foutz, M. S. Shur, and L. F. Eastman, "Steady-state and transient electron transport within bulk wurtzite indium nitride: An updated semiclassical three-valley Monte Carlo simulation anaysis," Appl. Phys. Lett., 87, 222103 (2005). 20. V. M. Polyakov and F. Schwierz, "Low-field electron mobility in wurtzite InN," Appl. Phys. Lett., 88, 032101 (2006). 21. V. M. Polyakov and F. Schwierz, "Nonparabolicity effect on bulk transport properties in wurtzite InN," J. Appl. Phys., 99, 113705 (2006). 22. C. G. Rodrigues et al. "Nonlinear transport properties of III-nitrides in an electric field," J. Appl. Phys., 98, 043702 (2005). 23. E. Starikov, et al. "Monte Carlo calculations of static and dynamic electron transport in nitrides," J. Appl. Phys., 98, 083701 (2005). 24. I. Vurgaftman, J. R. Meyer, and L. R. Ram-Mohan, "Band parameters for III­V compound semiconductors and their alloys," J. Appl. Phys., 89, No. 11, 5815­5875 (2001). 25. D. Zanato, N. Balkan, B. K. Ridley, G. Hill, and W. J. Schaff, "Hot electron cooling rates via the emission of LO-phonons in InN," Semicond. Sci. Technol., 19, 1024­1028 (2004). 26. J. M. Ziman, Principles of the Theory of Solids, Cambridge Univ. Press (1972). N. A. Masyukov and A. V. Dmitriev Moscow State University, Moscow, Russia E-mail: dmitriev@lt.phys.msu.su

823