Band structure

DAMOCLES employs a band structure obtained using empirical local or nonlocal pseudopotentials. For the local-pseudopotential calculation, see M. L. Cohen and T. K. Bergstresser, Phys. Rev. 141, 789 (1966), for the nonlocal-pseudopotential calculation, see J. R. Chelikowsky and M. L. Cohen, Phys. Rev. B 14, 556 (1976). The local approximation has been used originally. It is satisfactory in most cases and it still used almost exclusively for the conduction bands. The inclusion of the spin-orbit interaction -- necessary for the valence bands -- has required a re-working of the band structure and the nonlocal approximation has been adopted for the valence bands of those materials for which hole transport has been implemented. A transition to the nonlocal band structure will occur as necessity requires it: High-energy optical processes in Si, for example, can be better investigated with more accurate splittings at the symmetry points that nonlocal pseudopotentials provide. This is one of the reasons why the conduction bands using nonlocal pseudopotential are being recalculated and implemented in DAMOCLES (still in progress).

which is obtained by writing the Schrödinger equation using the expansion over Bloch waves

for the electron wavefunctions. The cutoff $G$_max is a numerical necessity and it corresponds to accounting only for waves such that:

The Hamiltonian matrix $H$_GG' includes four terms: a kinetic energy term, the local ionic pseudopotentials, a nonlocal contribution to the ion (pseudo)potentials, and the spin-orbit interaction:

where:

• the local form factors are expressed in terms of symmetric (V_S) and antisymmetric (V_A) form factors as follows:
  
  

  t being the vector (1,1,1)a/8

• S_s(G) is the atomic structure factor for the ion s:
  
  

  where N_a is the number of atomic species in the cell, R₀ is is location of the cell, and the sum is performed over the positions of all ions of the same species s within the cell

• the functions F_l^s are associated to the nonlocal wells (the p-well is usually ignored) and are expressed as:
  
  

  in terms of the spherical Bessel functions j_l, where d stands for radius of the well of l-symmetry, d_l^s, for the ion s

• A_l^s(E) is the energy-depth of the l-well for the ion s. Usually, the energy dependence of the d-well is ignored, while for the s-well the following empirical expression is used:
  
  

  where:

  
  

• the spin-orbit Hamiltonian, H^so_GG' is given by:
  
  

  where:

  
  

  where µ is an empirical parameter, µ₁ and µ₂ are the spin-orbit energies of the outermost p-type orbitals, and

  
  

  in terms of the radial functions of these orbitals, ß being a normalization constant fixed by the condition that B_s(k)/k approaches unity in the limit of k going to zero.

While only the band-structure information of in the 1/48 irreducible wedge of the first Brillouin Zone (BZ) is necessary, thanks to the 48-fold symmetry of the BZ, for programming simplicity and speed of execution, band-structure and carrier-phonon scattering rates are tabulated in the first octant. This is discretized into 5775 cubes with sides of length D_k=0.05, each cube identified by a k-vector k_ß (ß = 1, 5775), lying at the vertex closer to the origin. For each of these 5775 points, and for 26 additional vectors (for each of the 5775 points) displaced by 0.001, the quantities:

i.e., the energy and their gradients, and

for the 5 lowest-lying conduction bands and the 3 highest-lying valence bands are calculated and stored in look-up tables. Also stored are the mappings

which allow a fast location of any k-point in the BZ: The point is first mapped into the first octant, each of its coordinate is divided by the mesh-spacing D_k, obtaining 3 integers n_x, n_y, and n_z, after taking the integer parts and adding unity. The mapping above provides the index ß of the cube in which the point is located, since

is the k-point associated to that cube.

The energy and group velocity of a particle with crystal momentum k in band n is obtained first by mapping k into the BZ, if necessary, then into the first octant. The sign-transformations necessary to map the k point into the first octant are stored, since they are needed to convert back the group velocity to the original point outside the first octant. Finally, the following multi-linear interpolation formulae are used:

where the index ß now runs over the 8 vertices of the cube within which k lies, each weighted by:

and where a quadratic (linear) interpolation around each vertex is performed for the energy (velocity):

The 'direct' interpolation described above solves the relatively simple problem of obtaining the energy and group velocity which are needed to integrate the equations of motion during free-flights in the Monte Carlo algorithm. A much tougher problems, however, must be solved when computing scattering rates and when having to select the final state of a particle after a collision: In these cases one needs to find a particular k-vector at a given (final) energy.

To solve this problem (the most computing-intensive step in the entire DAMOCLES program), an additional look-up table is generated. The entire first octant of the BZ is discretized into two nested meshes: A coarse mesh of 108 cubes of side D_q=4D_k=0.2, and another mesh (the fine mesh), of 512 cubes within each cube of the coarse mech, each of side d_q=D_q/8. Thus, each cube in the fine mesh is labeled by a pair of indices such that the k-vector at the 'lower' vertex of the cube is:

As done before, a pair of mappings

is also defined and stored for fast location within the BZ of arbitrary k-points.

whereA_ßn is the area of the section obtained by cutting the cube with the constant-energy plane, expressed as:

where b is the half-size of the side of the cube, the coefficients a_j are the cosines of the angles formed by the gradient with respect to the cartesian axes, i.e.,

and the coefficients w_j are the dimensionless displacements of the k-point away from the center of the cube at which the equi-energy 'planes' intersect the vertices of the cube:

Thus, in this linear approximation, the maximum and minimum energies bracketing the energy-range spanned by the cube are:

The dimensionless parameter w_ß, expressing the shift from the center of the cube, k_ßc, along the direction of the gradient required to reach the plane of constant energy E is given by:

With this machinery it is now possible to perform integrations over equi-energy surfaces and to select states lying on these surfaces, as required by rates calculation and selection of final states after scattering. As an example, the density of states at energy E for all bands and both spin-states is given by:

which will be in joules^-1cm^-3, if the density-of-states S^(GR) is in 'states/(cube eV)'. In the following we shall always use these units, by measuring energies in eV and k-vectors in .

The rate at which an electron or hole of wavevector k in band n emits (lower sign) or absorb (upper sign) a phonon is computed using the adiabatic, harmonic, and weak-coupling approximations using Fermi Golden Rule:

where the carrier-phonon matrix element is given by:

for nonpolar scattering, and by the Fröhlich matrix element:

for polar scattering with longitudinal-optical phonons in the polar III-V materials. The nonpolar matrix elements are approximated as:

Scattering with ionized impurities is computed using the Born approximation with phase-shift corrections as given by J. R. Mayer and F. J. Bartoli, Phys. Rev. B 23, 5413 (1981):

where N_imp(r) is the concentrations of ionized impurities at location r, the function h(y,s) is the phase-shift correction tabulated by Mayer and Bartoli, as a function of dimensionless wavevector y (in Bohr units, taken positive for attractive and negative for repulsive interactions), and of the screening parameter s (rescaled in order to satisfy Friedel's sum rule). The Coulomb matrix element is given by:

where Ze is the charge of the impurities, I_G(k',n';kn) is the overlap factor, ß(q,0) is the static, wavevector dependent screening parameter, as given by A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems (McGraw-Hill, New York, 1971), p. 305.

During the Monte Carlo run, a 'self-scattering' rate is evaluated by computing the scattering rate ignoring overlap-factor effects. Using the usual discretization of th BZ discussed above, the 'self-scattering' rate takes the form:

Particles undergoing this fictitious scattering process are subject to another Monte Carlo step after collision: The collision is accepted or rejected with probability given by the overalp factor.

The short-range carrier-carrier interaction is treated similarly. Within the same approximations, the scattering rate for a particle of wavevector k in band n to collide with any of the 'partner particles' in the Monte carlo ensemble, each with wavevector p and in a band m, is given by:

where the Coulomb matrix element now has a slightly more complicated form, since the interaction may involve indistinguishable particles (e-e and h-h scattering). In these cases the matrix element is given by the sum of the scattering amplitudes for pairs in the triplet and singlet states of total angular momentum:

where

The long-range Coulomb interactions can be treated either as an additional Monte Carlo scattering process, introducing carrier-plasmon scattering, or by letting the self-consistent Poisson/particles procedure account for plasma effects. The second alternative is employed in DAMOCLES. The reason behind this choice is twofold:

1. Plasmons can decay via collisional damping (most efficient in materials with high low-energy scattering, such as Si, Ge, or GaP) or via Landau damping (i.e., the decay of short wavelength plasmons back into single particle excitations, most efficient in high-mobility materials, such as GaAs, InAs, etc.). The latter process is responsible for re-inserting the energy and momentum lost by the carriers to plasmons back into the carrier-gas. If carrier-plasmon coupling is treated as a Monte Carlo scattering process, it is imperative to account for the Landau damping self-consistently. This is a hard problem. On the other side, the self-consistent Poisson/Monte Carlo procedure automatically accounts for this effect.
2. Since the self-consistent Poisson/Monte Carlo procedure already couples carriers at long range, including carrier-plasmon scattering as a Monte Carlo collision process would result in double counting the long-range Coulomb interaction, unless proper measures were taken to 'dampen' fluctuations of the electrostatic potential during the Poisson-solving step. Therefore, it is more convenient to let the Poisson solution take care of the long-range interactions in a natural way.

where n is the maximum carrier density, and m_c is the conductivity mass at the band-edge. Also, in high-doping regions where plasmon effects are significant, the size of the spatial mesh must be of the order of the inverse local screening length and the number of simulated particles must be properly chosen, as explained in M. V. Fischetti and S. E. Laux, Phys. Rev. B 38, 9721 (1988) and in S. E. Laux and M. V. Fischetti, in "Monte Carlo Device Simulation: Full Band and Beyond", Karl Hess editor (Kluwer Academic, Boston, 1991), p 1.

	damoclesNO-SPAM@watson.ibm.com (last updated: January 26, 1999)