Chapter 3
Physical Modeling

The SHE equations (2.34) and (2.35) incorporate material-specific properties by the velocity term v , the modulus |k| of the wave vector as a function of energy, the generalized density of states Z  , and the scattering operator Q  . However, the velocity term v and the density of states Z  are not independent and depend on the dispersion relation ε(k)  by

       v = ℏ ∇kε ,                             (3.1)
Z(ε,θ,φ) = ---3---- .                          (2.3)
           4π   ∂ε
Similarly, the modulus of the wave vector is obtained from inverting the dispersion relation. Consequently, the choice of the dispersion relation plays a central role for the accuracy of the SHE method. Spherically symmetric approximations are presented in Sec. 3.1, but only poor approximations are obtained at higher energies. Sec. 3.2 deals with the state-of-the-art on the inclusion of full-band effects up to high energies.

The second half of this chapter is devoted to the various scattering effects. Since scattering balances the energy gain of carriers due to the electric field, accurate expressions for the scattering operators are mandatory. Scattering mechanisms leading to a linear operator are discussed in Sec. 3.3, while the case of nonlinear scattering operators is investigated in Sec. 3.4.

3.1 Spherically Symmetric Energy Band Models

Due to the nonuniformity of the crystal lattice in silicon with respect to a change in direction, the dispersion relation linking the particle momentum with the particle energy cannot be accurately described by the idealized setting of an infinitely deep quantum well, for which the solution of Schrödinger’s equation yields a parabolic dependence of the energy on the wave vector:

        2   2
ε(k) = ℏ-|k|-,                              (3.2)
        2m ∗
where ℏ  is the scaled Planck constant and m ∗ is the effective mass. Still, this quadratic relationship termed parabolic band approximation is a good approximation near the minimum of the energy valley. Due to its simple analytical form, the parabolic dispersion relation is often used up to high energies, for which it fails to provide an accurate description of the material.

A more accurate approximation to the band structure in silicon can be obtained by a slight modification of the form

        2 2
γ(ε) = ℏ-k--,                              (3.3)
       2m ∗
where the typical choice
γ(ε) = ε(1+ α ε)                             (3.4)
is known as Kane’s model  [55]. The parameter α  is called nonparabolicity factor; in the case α = 0  one obtains again (3.2). More complicated analytical dependencies of the energy on the wave vector are certainly possible, but (3.4) already provides very good results in the low-energy regime. The choice α = 0.5  provides a good approximation of the dispersion relation for electrons in relaxed silicon. However, for kinetic energies above 1.75  eV the nonparabolic approximation fails to describe the nonmonotonicity of the density of states in silicon.


Figure 3.1: Multi-band structure of silicon. Bands 1 and 3 are electron-like (increasing density of states), bands 2 and 4 are hole-like (decreasing density of states).

The deficiencies of the nonparabolic dispersion relation can be mitigated by a combination of four energy bands as proposed by Brunetti et. al.  [11]:

           2  2
ε+ α ε2 = -ℏ-k-- ,                       (Band  1)            (3.5)
         2m ∗(1)
          (2)   -ℏ2k2-
     ε = εmax − 2m ∗(2) ,                (Band  2)            (3.6)
     ε = ε(m3i)n + ---∗(3) ,                (Band  3)            (3.7)
     ε = ε(4) − -ℏ2k2- .                (Band  4)            (3.8)
          max   2m ∗(4)
The bands are valid up to energy limits  (2)
εmax   ,  (3)
εmin   and  (4)
εmax   respectively. The specific form of each band is given by
γ(1)(ε) = ε + αε2 ,                  (Band  1)              (3.9)
 (2)      (2)
γ  (ε) = εmax − ε ,                 (Band  2)             (3.10)
 (3)          (3)
γ  (ε) = ε − εmin ,                 (Band  3)             (3.11)
γ(4)(ε) = ε(m4a)x − ε .                 (Band  4)             (3.12)
The relationships between the energy and the wave vector can be written in compact form as
γ (ν)(ε) = --∗(ν) ,   ν = 1,2,3,4 ,                   (3.13)
and are depicted in Fig. 3.1. The total density of states Z(ε)  is computed from the individual densities of states by a weighted sum:
Z(ε) = β (1)Z (1)(ε) + β(2)Z (2)(ε)+ β(3)Z(3)(ε)+ β (4)Z (4)(ε) ,       (3.14)
where the band multiplicities β (ν)   , ν = 1,2,3,4  account for the number of equivalent symmetrical bands of the ν  -th band can be found together with the other parameters of the multi-band model in Tab. 3.1.

Band  ∗
m ∕m0   εmin   εmax   β  α

1 0.320 0 1.75 6 0.35
2 0.712 1.75 3.02 6 0
3 0.750 2.60 3.00 12 0
4 0.750 3.00 3.40 12 0

Table 3.1: Parameters for the analytical band structure proposed by Brunetti et. al.  [11]. m
  0   denotes the electron mass, β  the multiplicity of the respective band and α  is the nonparabolicity factor. Each band extends from εmin   to εmax   .

The BTE has in principle to be solved in each energy band with index ν  for a distribution function   (ν)
f   . However, since the individual distribution functions for each energy band are not of particular interest, it is preferred to have a single dispersion relation describing the total distribution function f  .

3.2 Incorporation of Full-Band Effects

In order to account for the angular dependence of the energy on the wave vector in a semiconductor, the inverse dispersion relation is expanded as

           ∑L  ∑l
k(ε,θ,φ) =         kl,m(ε)Y l,m ,                    (3.15)
           l=0 m= −l
where due to symmetries of the band structure only even values of l  and nonnegative m  -values, which are a multiple of four, lead to nonzero expansion coefficients. The expansion coefficients kl,m  can either be determined by an integration over spheres in k -space, cf. (2.2), or by minimization of certain target functionals such as the quadratic error of the scaled moments
 V0(ε) :=  (2π)3 B δ(ε− ε(k))dk ,                          (3.16)
             1     ∫
Vn (ε) :=  ----3-----   |v(k )|nδ(ε− ε(k))dk ,               (3.17)
         (2π) V0(ε) B
where B  denotes the Brillouin zone.

A SHE of the valence band has been carried out by Kosina et al.  [60] up to 1.27  eV. Pham et al.  [73] proposed a refined method for an expansion also including higher energies. A fitted band structure based on SHE for the conduction band has been presented by Matz et al.  [69]. Due to the bijective mapping between energy and wave vector, the velocity and the density of states are not in perfect agreement with full-band data. Nevertheless, good results compared to the Monte Carlo method are obtained  [45].

(a) Density of States.
(b) Carrier Velocity.

Figure 3.2: Comparison of the density of states and the carrier velocity for the band models (3.2), (3.3), (3.5)-(3.8), and (3.15) with full-band Monte Carlo values from  [69]. The nonparabolic band (3.3) using (3.4) and α =  0.5  is labelled Modena model, as it is also used in e.g.  [42].

A comparison of the presented energy band models is given in Fig. 3.2. Slight deviations of the many-band model and the full-band density of states are due to different Monte Carlo data used in  [11]. While the many-band model provides a good fit for the density of states, it fails to approximate the carrier velocity and is worse than the nonparabolic model (3.3). As expected, the fitted band model provides the best accuracy, even though approximations are less accurate at energies above 2  eV.

Vecchi et al.  [107] used full-band Monte Carlo data for the velocity v and the generalized density of states Z  for a first-order SHE. This is possible because in this case all terms with an explicit representation of the band structure, i.e. Γ l′,m ′
  0,0   , vanish. However, higher-order SHE does not allow for a similar procedure, because the angular coupling term (2.33) does not vanish any longer and an explicit expression for the modulus of the wave vector is required.

Recently, Jin et al.  [48] suggested a reformulation as follows: Consider

-Z--   -|k|-∂-|k|
ℏ|k | = 4π3ℏ ∂ε  .                           (3.18)
With the relation                  2
2|k|∂|k |∕∂ε = ∂|k| ∕∂ε  , the expression can be further rearranged to
   (       )
∂--   1-∂ε-
∂ε  Z ℏ∂ |k |  .                            (3.19)
With (3.1), one thus obtains
-Z--≈ 1-∂|v|Z-,                            (3.20)
ℏ|k |  2  ∂ ε
where equality holds true only in a homogeneous material with a spherically symmetric dispersion relation. Nevertheless, the results presented in  [48] demonstrate that the use of full-band data for the velocity and the density of states using the approximation (3.20) leads to an accuracy comparable with the anisotropic fitted band (3.15). Consequently, the coupling term (2.33) is then approximated by
               ∫  (                       )
 l′,m′   1∂-|v|Z-     ∂Y-l,m-    --1- ∂Yl,m--     l′,m ′
Γl,m  ≈  2 ∂H    Ω    ∂θ  eθ + sinθ  ∂φ  eφ  Y     dΩ          (3.21)
and similarly for the adjoint term (2.36). The integral can be tabulated easily, since it does not contain any band structure information any longer. As will become clear in Chap. 4, the use of band structure parameters without angular dependencies brings a number of additional benefits from the computational point of view compared to the fitted band structure (3.15).

3.3 Linear Scattering Operators

While the band structure links the particle energy with the particle momentum, it does not fully describe the propagation of carriers. In the presence of an electrostatic force, carriers would be accelerated and thus gain energy indefinitely unless scattering with the crystal lattice or with other carriers is included in the model. The most important scattering mechanisms are discussed in the following. The scattering operator is assumed to be given in the form

        --1--         ′        ′            ′           ′
Q {f} = (2π)3  Bs(x,k ,k)f(x, k,t)− s(x, k,k )f (x, k,t)dk ,       (2.13)
where it has to emphasized that scaling factors in front of the scattering integral, here       3
1∕(2π )   , need to be taken into account when comparing scattering rates from different sources. Note that the numerator of the prefactor differs from the numerator used for the spherical projection (2.6) due to the assumption that scattering does not change spin. Moreover, it should be noted again that the commonly written small sample volume Vs   as prefactor for the scattering integral is not written explicitly in the following.

3.3.1 Acoustic Phonon Scattering

Atoms in the crystal lattice vibrate around their fixed equilibrium locations at nonzero temperature. These vibrations are quantized by phonons with energy ℏωphon   . Acoustic vibrations refer to a coherent movement of the lattice atoms out of their equilibrium positions. Depending on the displacements with respect to the direction of propagation of the lattice wave, transversal (TA) and longitudinal (LA) acoustic modes are distinguished.

Since the change in particle energy due to acoustic phonon scattering is very small, the process is typically modelled as an elastic process  [47], which does not couple different energy levels. The scattering rate can thus be written as

sac(x,k, k′) = σacδ(ε(k )− ε(k′)) ,

where the coefficient σac   is given by

σac = 2πkBT-E--,                            (3.23)
where E is the deformation potential, ρ  is the density of mass, and ul   is the longitudinal sound velocity, cf. Tab. 3.2.

Si Ge

ρ  2.33  g/cm3   5.32 g/cm3
ul           5
9.05 × 10   cm/s         5
5.40× 10   cm/s
ϵ  11.7 × ϵ
       0   16.0× ϵ
Ee   8.90  eV 8.79  eV
Eh   5.12  eV 7.40  eV

Table 3.2: Material parameters for silicon and germanium, cf.  [52]. The subscripts e  and h  are used to distinguish between electrons and holes.

3.3.2 Optical Phonon Scattering

Optical phonon scattering refers to an out-of-phase movement of lattice atoms. In ionic crystals, these vibrations can be excited by infrared radition, which explains the name. Similar to acoustic phonon scattering, transversal (TO) and longitudinal (LO) modes are distinguised.

Since the involved phonon energies are rather high, cf. Tab. 3.3, optical phonon scattering is typically modeled as an inelastic process leading to a change of the particle energy. With the phonon occupation number Nphon   given by the Bose-Einstein statistics

N     = -------1-------,                        (3.24)
  phon   exp (ℏωphon)− 1
the scattering rate for the initial state k and the final state  ′
k can be written as
sop(x,k,k′) = σopNphon  δ(ε(k) − ε(k′) + ℏωop)
           + (1+ N     )δ(ε(k) − ε(k′) − ℏω  )] ,
                   phon                  op

where           ′
σop(x,k, k)  is symmetric in k and   ′
k and given by

σop =    ρων    ,                           (3.26)
with coupling constant D  K
  t ν  , mass density ρ  , and phonon frequency ω
 ν  . Values for the individual modes can be found in Tab. 3.3. It should be noted that optical phonon scattering couples the energy levels H  − ℏωop   , H  and H + ℏωop   in an asymmetric manner, because scattering from higher energy to lower energy is more likely than vice versa.


ν Mode DtK ν  ℏω  DtK ν  ℏω

1  TA 0.470× 108   eV/cm 12.1  meV 0.479 × 108   eV/cm 5.60  meV
2  LA 0.740× 108   eV/cm 18.5  meV 0.772 × 108   eV/cm 8.60  meV
3  LO 10.23× 108   eV/cm 62.0  meV 9.280 × 108   eV/cm 37.0  meV
4  TA 0.280× 108   eV/cm 19.0  meV 0.283 × 108   eV/cm 9.90  meV
5  LA 1.860× 108   eV/cm 47.4  meV 1.940 × 108   eV/cm 28.0  meV
6  TO 1.860× 108   eV/cm 58.6  meV 1.690 × 108   eV/cm 32.5  meV

Table 3.3: Modes, coupling constants and energies for inelastic phonon scattering with electrons, cf.  [52]. (TA: transveral acoustic, LA: longitudinal acoustic, TO: transversal optical, LO: longitudinal optical)

3.3.3 Ionized Impurity Scattering

Dopants in a semiconductor are fixed charges inside the crystal lattice. Since carriers are charged particles as well, their trajectories are influenced by these fixed charges, leading to a change of their momentum. The model by Brooks and Herring  [1047] suggests an elastic scattering process with scattering coefficient

           ′             ′            ′
simp(x,k,k ) = σimp(x,k,k )δ(ε(x )− ε(x )) ,              (3.27)
where            ′
σimp(x,k, k)  is symmetric in k and   ′
k and given by
           ′   2π-NIq4---------1---------
σimp(x,k,k ) = ℏ   ϵ2 ((k ′ − k)2 + 1∕λ2 )2 .            (3.28)
The ordinality of the impurity charge is assumed to be one, and Na   and Nd   denote the acceptor and donator concentrations respectively. The Debye length λD   under assumption of local equilibrium is given by
λD =  q2(n+ p) .                            (3.29)
Similar to elastic acoustical phonon scattering, ionized impurity scattering does not couple different energy levels. However, a considerable complication stems from the angular dependence of the coefficient σimp   . This complication can be circumvented by approximating the anisotropic coefficient (3.27) by an elastic-isotropic process with the same momentum relaxation time τm;ii    [52]. The momentum relaxation time is computed for an isotropic dispersion relation by an integration over the whole Brillouin zone and by weighting the change of direction of the momentum  [67]:
τm;ii(k ) =  1
(2π ) Bsimp(x,k,k)(1 cos(θ))dk3
Here, the z  -axis for the integration in the Brillouin zone is chosen such that it is aligned with k , hence the angle between k and k′ is given by the inclination θ  . Transformation to spherical coordinates leads to
τm;ii(k ) =   3   4
  ℏϵ2 0 0π(------1−--cos(θ)-----)-
 4|k|2 sin2(θ∕2)+ 1 ∕λ2D 2 sinθdθZdε
where the density of states Z  is independent of the angles because of the assumption of an isotropic dispersion relation. The integral over the inclination θ  can be computed analytically as
 (k ′ − k)2 + 1∕λ2D sinθdθ = --1-4
4|k|[                     2   2  ]
  ln(1 + 4λ2D|k|2)− --4λD-|k2|---
                  1 + 4λD|k|2 .
Therefore, the isotropic scattering coefficient
                      4      [                    2   2  ]
σimp;iso(x,k,k ′) = π-NIq---1--- ln(1 + 4λ2|k|2)− --4λD-|k-|--        (3.30)
                 ℏ  ϵ2  4|k |4          D       1 + 4λ2D |k|2
has the same momentum relaxation time as the anisotropic coefficient (3.27). Note that |k| should be evaluated consistently with the approximated band, therefore the transformation (3.20) needs to be employed for the full-band case. Moreover, since the Brooks-Herring model fails to correctly describe the carrier mobility at high doping concentrations, an empirical fit factor depicted in Fig. 3.3 is usually employed additionally  [52] in order to reproduce the Caughey-Thomas expression for the mobility  [15].


Figure 3.3: Fit factor for ionized impurity scattering in order to reproduce the Caughey-Thomas expression for the mobility  [5215].

3.4 Nonlinear Scattering Operators

The scattering operator in low-density approximation for single carrier processes is linear, which is very attractive from a computational point of view. Consequently, scattering processes leading to a nonlinear scattering operator are often neglected in order to avoid nonlinear iteration schemes. In the following, two such types of scattering mechanisms are considered.

3.4.1 Pauli Exclusion Principle

A high population of the conduction band can lead to the case that the distribution function takes large values near the band edge, thus the term f(1− f)  cannot be approximated with f  any longer:

          1   ∫       ′        ′
Q {f} = (2π)3   s(x,k ,k)f(x,k ,t)(1− f(x, k,t))
               B        ′                  ′     ′
              − s(x,k,k )f(x,k, t)(1− f (x,k ,t))dk

Repeating the steps from Sec. 2.2, one finally obtains for the projected in-scattering operator assuming velocity randomization and a scattering rate independent of the angles

                          [              ][f]  (x,ε±  ℏω ,t)
Qiηn,l,m{f } = σ η(x, ε± ℏω η,ε) Zl,m(ε)− [f]l,m---0,0---------η---.        (3.32)
Similarly, the projected out-scattering operator is found as
                          [                                ]
Qout {f } = σ (x,ε,ε± ℏω  )Z   (ε± ℏω  )− [f]  (x, ε± ℏω  ,t) [f]l,m .     (3.33)
  η,l,m        η           η   0,0       η     0,0         η    Y0,0
The quadratic nonlinearity needs to be handled by a suitable nonlinear iteration scheme. In  [43] it was reported that no complications arose within a Newton scheme.

Numerical results in  [43] confirm that the low-density approximation does not have a high impact on macroscopic quantities such as the electron density or carrier velocities, but notable differences in the distribution function are obtained near the band edge. The carrier population is then shifted towards higher energies, because all states at lower energies are already populated.

3.4.2 Carrier-Carrier Scattering

The linear scattering operators in Sec. 3.3 stem from the scattering of carriers with noncarriers. Very important for particularly the high energy tail of the distribution function is carrier-carrier scattering  [79]. A carrier-carrier scattering mechanism requires that the two source states are occupied, and the two final states after scattering are empty. This leads to a scattering operator of the form

              ∫ ∫  ∫
        --1--              ′    ′         ′                     ′
Q {f} = (2π)3  B B  Bs(x,k ,k,k 2,k2)f(x,k ,t)(1 − f(x,k,t))f(x,k2,t)(1− f(x,k2, t))
                ′     ′                   ′                      ′     ′      ′
      − s(x,k, k,k2, k2)f(x,k,t)(1− f(x,k ,t))f(x,k2,t)(1 − f(x,k 2,t))dk dk2dk 2 ,

where the scattering coefficient s(⋅,⋅,⋅,⋅,⋅)  now depends on the spatial location and on two pairs of initial and final states. With a low-density approximation, the nonlinearity of degree four of the carrier-carrier scattering operator reduces to second order:

               ∫  ∫ ∫
Qcc{f } = --1--        s(x,k′,k,k′,k2)f(x, k′,t)f(x,k ′,t)
          (2π)3 B  B  B          2                   2
                    − s(x,k,k′,k2,k′)f(x,k, t)f(x,k2, t)dk ′dk2dk ′
                                   2                           2

The scattering coefficient can be derived to be of the form  [9699]

s(x,k′,k,k′2,k2 ) = σcc(x,k,k′,k2,k′2)δ(k + k′ − k2 − k′2)δ(ε+ ε′ − ε2 − ε′2) , (3.36)
where the two delta distributions in (3.36) refer to conservation of momentum and energy respectively, and
                         2π q4n(x )        1
σcc = σcc(x,k,k ′,k2,k ′2) =-----2------′----2------2-2 .        (3.37)
                          ℏ   ϵ    ((k  − k) + 1∕λ D)
The set of parameters is similar to that of ionized impurity scattering, with the impurity density Na + Nd   replaced by the carrier density n  . Since σcc   only depends via β  on the difference of the initial and the final state of one of the two carriers, the shorthand notation σcc(x, k,k′)  is used in the following.

Carrier-carrier scattering, particularly electron-electron scattering, has so far been discussed for first-order SHE only  [108110]. In a joint work with Peter Willibald Lagger  [62], the author has recently extended the method to arbitrary-order SHE, and the derivation is given in the following.

The scattering operator (3.35) is again split into an in-scattering term Qicnc   and an out-scattering term Qocuct   as in Sec. 2.2. Inserting (3.36) into (3.35), one integration can be carried out due to the momentum conservation. An integration with respect to k2   yields

           ∫  ∫
Qin=  --1--     σ  (x,k, k′)f(k ′)f(k′)δ(ε+ ε∗ − ε′ − ε′)dk ′dk ′ ,      (3.38)
 cc   (2π)3 B  B  cc                2                2      2
with the initial energy ε∗ := ε(k′ + k′ − k)
            2  of the second particle involved in the scattering process. A transformation of the two integrals to spherical coordinates leads to
      ∫   ∫   ∫  ∫
  in     ∞   ∞                ′   ′    ′     ′    ′
Q cc =  0   0   Ω  Ωσcc(x,k,k )f(k )f(k2)Z(k )Z(k 2)
                            ∗    ′   ′   ′ ′   ′  ′
                    × δ(ε+ ε  − ε − ε2)d εdε2dΩ dΩ 2 ,

where k = k (ε,θ,φ )  and similarly for   ′
k and  ′
k2   . A projection onto the spherical harmonic Y l,m  yields

                   ∫ ∞ ∫ ∞ ∫  ∫  ∫
Qincc,l,m{f}(x, ε,t) =                 σcc(x,k,k ′)f(k′)f(k′2)Z(k ′)Z (k ′2)
                    0   0   Ω  Ω  Ω
                           × δ(ε + ε∗ − ε′ − ε′2)Y l,mZd ΩdΩ ′d Ω′2dε′dε′2 .

Up to now, no approximations have been applied. However, a direct evaluation of these nested integrals at each node in the simulation domain is certainly prohibitive for the use within a simulator due to excessive execution times. Consequently, the further derivation is based on the following two assumptions:

With these assumptions, one can split the integrands to

            ∫ ∞ ∫ ∞
Qicnc,l,m {f} =         σcc(x, ε,ε′)δ(ε+ ε∗ − ε′ − ε′2)
             0   0    ∫  ∫  ∫
                                  ′    ′    ′    ′           ′   ′  ′ ′
                    ×  Ω  Ω  Ω f(k)f (k 2)Z (k )Z (k2)Yl,mZd Ωd Ω dΩ 2d εdε2 .

An expansion of f  into spherical harmonics, the use of the delta distribution for an elimination of the integral over ε2   , and the assumption of spherical energy bands leads to

              ∑   ∑   ∫ ∞
Qincc,l,m {f} = Z             ccc(x, ε,ε′)fl,m(ε′)fl′,m′(ε+ ε∗ − ε′)Z (ε′)Z(ε + ε∗ − ε′)dε′
              l′,m ′l′,m ′ 0
                  2  2  ∫         ∫          ∫
                      ×    Yl,md Ω   Yl′,m ′dΩ ′  Yl′,m ′dΩ′2 .
                         Ω         Ω          Ω  2  2

Summing up, the projected in-scattering operator is given by

             Z  ∫  ∞
Qincc,l,m{f} = --3-     σcc(x,ε,ε′)f0,0(ε′)f0,0(ε+ ε∗ − ε′)Z (ε)Z (ε+ ε∗ − ε′)dε′ .
            Y 0,0  0

For the more general case of an expansion of the scattering coefficient (3.41), one obtains again with the assumption of spherical energy bands

                  ∫ ∞
Qin   {f} = -1--Z    σ      ′ ′(x,ε,ε′)f   (ε′)f′  ′(ε + ε∗ − ε′)
 cc,l,m       Y0,0   0   cc;l,m,l,m         l,m     l,m
                      × Z (ε)Z (ε + ε∗ − ε′)dε′ .

The projection of the out-scattering operator starts with the same steps as for the in-scattering operator in order to arrive at

               ∑       ∑   ∫ ∞ ∫ ∞
Qouctc,l,m{f } = Z    fl′,m ′            σcc(x,ε,ε′)Z (ε′)Z (ε′2)
              l′,m′     l2,m2 0   0
                    ∫  ∫  ∫         ∗       ∗   ′   ′
                  ×          fl2,m2(ε )δ(ε + ε − ε − ε2)
                      Ω Ω  Ω                ∗   ∗      ′  ′  ′ ′
                           × Yl,mYl′,m ′Yl2,m2(θ ,φ )dΩd Ω dΩ 2dεdε2 .

The angles θ∗ and φ∗ depend in a complicated way on the other angles, therefore rather crude simplifications are applied. Scattering with another carrier is to a first approximation determined by the density of carriers at the particular location inside the device. Consequently, the distribution function of the second particle, which is described by  ∗
ε and  ∗  ∗
θ ,φ , is approximated by the isotropic part of the distribution function only, since it fully describes the carrier density. With    ∗       ′   ′             ∗
f(k ) = f(k + k2 − k) ≈ f0,0(ε )Y0,0   and considering  ∗
ε again as the average energy, (3.46) simplifies to

                             ∫ ∞
Qout  {f} = -Z--f  (ε∗)f   (ε)    σ  (x,ε,ε′)Z (ε′)Z(ε + ε∗ − ε′)d ε′ .
 cc,l,m       Y03,0 0,0    l,m     0   cc

The more general case of an expansion of the scattering coefficient (3.41) results in

Qout  {f } = Z--f0,0(ε∗)    fl′,m ′
  cc,l,m       Y20,0       l′,m′
               ∑    ∫ ∞
            ×           σcc;lcc,mcc,0,0(x, ε,ε′)Z(ε′)Z(ε+ ε∗ − ε′)dε′
              lcc,mcc  0
            ×    Yl,mYl′,m ′Ylcc,mccdΩ .

For the full projected scattering operator, one thus obtains

             Z  ∫ ∞           [
Qcc,l,m {f} = --3-    σcc(x,ε,ε′)f0,0(ε′)f0,0(ε+ ε∗ − ε′)
            Y0,0 0                             ]
                                − f0,0(ε∗)fl,m (ε)Z (ε)Z(ε+  ε∗ − ε′)dε′ .

With the use of an expansion of the scattering coefficient (3.41), one obtains

                 ∑     ∑   ∫ ∞
Qcc,l,m {f} = -Z--               σcc;l ,m  ,l′,m′ (x, ε,ε′)
            Y02,0l ,m   l′,m′  0     cc  cc cc  cc
                 cc  cccc cc    [
                                Y0,0flcc,mcc(ε′)fl′cc,m′cc(ε + ε∗ − ε′)
                                         ∑       ∫                   ]
                               − f0,0(ε∗)    fl′,m′   Yl,mYl′,m ′Ylcc,mccdΩ
                                        l′,m′      Ω
                                            ∗    ′  ′
                               × Z (ε)Z(ε+ ε  − ε)dε  .

One can immediately see that the full scatter operator vanishes for the equilibrium case, where f  is given by a Maxwell distribution. Therefore, even though simplifications have been used for the separate derivation of the projected equations for the in- and the out-scattering operators, the resulting expressions are consistent.