Band Structure Calculations in Semiconductor Manufacturing

Mathematical Framework

1. The Fundamental Problem

We need to solve the many-body Schrödinger equation for electrons in a crystal:

$$ \hat{H}\Psi = E\Psi $$

The full Hamiltonian includes kinetic energy, ion-electron interaction, and electron-electron repulsion:

$$ \hat{H} = -\sum_i \frac{\hbar^2}{2m} abla_i^2 + \sum_i V_{\text{ion}}(\mathbf{r}_i) + \frac{1}{2}\sum_{i eq j} \frac{e^2}{|\mathbf{r}_i - \mathbf{r}_j|} $$

Key challenges:

• The system contains ~$10^{23}$ electrons
• Electron-electron interactions couple all particles
• Analytical solution is impossible for real materials
• Requires a hierarchy of approximations

2. Density Functional Theory (DFT)

The workhorse of modern band structure calculations rests on the Hohenberg-Kohn theorems:

1. Ground-state properties are uniquely determined by electron density $n(\mathbf{r})$ 2. The true ground-state density minimizes the energy functional

2.1 Kohn-Sham Equations

The many-body problem is mapped to non-interacting electrons in an effective potential:

$$ \left[-\frac{\hbar^2}{2m} abla^2 + V_{\text{eff}}(\mathbf{r})\right]\psi_i(\mathbf{r}) = \epsilon_i\psi_i(\mathbf{r}) $$

where the effective potential is:

$$ V_{\text{eff}}(\mathbf{r}) = V_{\text{ion}}(\mathbf{r}) + V_H(\mathbf{r}) + V_{xc}[n] $$

Components of $V_{\text{eff}}$:

• Ionic potential: $V_{\text{ion}}(\mathbf{r})$ — interaction with nuclei
• Hartree potential: $V_H(\mathbf{r}) = \int \frac{n(\mathbf{r}')}{|\mathbf{r}-\mathbf{r}'|}d\mathbf{r}'$ — classical electrostatic repulsion
• Exchange-correlation: $V_{xc}[n] = \frac{\delta E_{xc}[n]}{\delta n(\mathbf{r})}$ — quantum many-body effects

The density is reconstructed self-consistently:

$$ n(\mathbf{r}) = \sum_i^{\text{occupied}} |\psi_i(\mathbf{r})|^2 $$

2.2 Exchange-Correlation Functionals

The unknown piece requiring approximation:

• Local Density Approximation (LDA):

$$ E_{xc}^{\text{LDA}}[n] = \int n(\mathbf{r})\,\epsilon_{xc}^{\text{homog}}(n(\mathbf{r}))\,d\mathbf{r} $$

• Generalized Gradient Approximation (GGA):

$$ E_{xc}^{\text{GGA}}[n] = \int f\left(n(\mathbf{r}), abla n(\mathbf{r})\right)\,d\mathbf{r} $$

• Hybrid Functionals (HSE06):

$$ E_{xc}^{\text{HSE}} = \frac{1}{4}E_x^{\text{HF,SR}}(\mu) + \frac{3}{4}E_x^{\text{PBE,SR}}(\mu) + E_x^{\text{PBE,LR}}(\mu) + E_c^{\text{PBE}} $$

• Mixing parameter: $\alpha = 0.25$
• Screening parameter: $\mu \approx 0.2\,\text{Å}^{-1}$

3. Bloch's Theorem and Reciprocal Space

For a periodic crystal with lattice vectors $\mathbf{R}$, the fundamental symmetry relation:

$$ \psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}}\,u_{n\mathbf{k}}(\mathbf{r}) $$

where:

• $u_{n\mathbf{k}}(\mathbf{r})$ has lattice periodicity: $u_{n\mathbf{k}}(\mathbf{r} + \mathbf{R}) = u_{n\mathbf{k}}(\mathbf{r})$
• $\mathbf{k}$ is the crystal momentum (wavevector)
• $n$ is the band index

3.1 Reciprocal Lattice

Reciprocal lattice vectors $\mathbf{G}$ satisfy:

$$ \mathbf{G} \cdot \mathbf{R} = 2\pi m \quad (m \in \mathbb{Z}) $$

For a cubic lattice with parameter $a$:

$$ \mathbf{G} = \frac{2\pi}{a}(h\hat{\mathbf{x}} + k\hat{\mathbf{y}} + l\hat{\mathbf{z}}) $$

The band structure $E_n(\mathbf{k})$ emerges as eigenvalues indexed by:

• Band number $n$
• Wavevector $\mathbf{k}$ within the first Brillouin zone

4. Basis Set Expansions

4.1 Plane Wave Basis

Expand the periodic part in Fourier series:

$$ u_{n\mathbf{k}}(\mathbf{r}) = \sum_{\mathbf{G}} c_{n,\mathbf{k}+\mathbf{G}}\,e^{i\mathbf{G}\cdot\mathbf{r}} $$

The Schrödinger equation becomes a matrix eigenvalue problem:

$$ \sum_{\mathbf{G}'} H_{\mathbf{G},\mathbf{G}'}(\mathbf{k})\,c_{\mathbf{G}'} = E_{n\mathbf{k}}\,c_{\mathbf{G}} $$

Matrix elements:

$$ H_{\mathbf{G},\mathbf{G}'} = \frac{\hbar^2|\mathbf{k}+\mathbf{G}|^2}{2m}\delta_{\mathbf{G},\mathbf{G}'} + V(\mathbf{G}-\mathbf{G}') $$

Basis truncation via kinetic energy cutoff:

$$ \frac{\hbar^2|\mathbf{k}+\mathbf{G}|^2}{2m} < E_{\text{cut}} $$

Typical values: $E_{\text{cut}} \sim 30\text{--}80\,\text{Ry}$ (400–1000 eV)

4.2 Localized Basis (LCAO/Tight-Binding)

Linear Combination of Atomic Orbitals:

$$ \psi_{n\mathbf{k}}(\mathbf{r}) = \sum_{\alpha} c_{n\alpha\mathbf{k}} \sum_{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}}\phi_\alpha(\mathbf{r} - \mathbf{R} - \mathbf{d}_\alpha) $$

This yields a generalized eigenvalue problem:

$$ H(\mathbf{k})\,\mathbf{c} = E(\mathbf{k})\,S(\mathbf{k})\,\mathbf{c} $$

where:

• $H_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}}\langle\phi_i(\mathbf{r})|\hat{H}|\phi_j(\mathbf{r}-\mathbf{R})\rangle$ — Hamiltonian matrix
• $S_{ij}(\mathbf{k}) = \sum_{\mathbf{R}} e^{i\mathbf{k}\cdot\mathbf{R}}\langle\phi_i(\mathbf{r})|\phi_j(\mathbf{r}-\mathbf{R})\rangle$ — Overlap matrix

4.3 Slater-Koster Parameters

For empirical tight-binding with direction cosines $(l, m, n)$:

$$ \begin{aligned} E_{s,s} &= V_{ss\sigma} \\ E_{s,x} &= l \cdot V_{sp\sigma} \\ E_{x,x} &= l^2 V_{pp\sigma} + (1-l^2) V_{pp\pi} \\ E_{x,y} &= lm(V_{pp\sigma} - V_{pp\pi}) \end{aligned} $$

Harrison's universal parameters:

5. Pseudopotential Theory

Core electrons are chemically inert but computationally expensive. Replace true potential with smooth pseudopotential.

5.1 Norm-Conserving Conditions

(Hamann, Schlüter, Chiang):

1. Matching: $\psi^{\text{PS}}(r) = \psi^{\text{AE}}(r)$ for $r > r_c$ 2. Norm conservation: $$ \int_0^{r_c}|\psi^{\text{PS}}(r)|^2 r^2 dr = \int_0^{r_c}|\psi^{\text{AE}}(r)|^2 r^2 dr $$ 3. Eigenvalue matching: $\epsilon^{\text{PS}} = \epsilon^{\text{AE}}$ 4. Log-derivative matching: $$ \left.\frac{d}{dr}\ln\psi^{\text{PS}}\right|_{r_c} = \left.\frac{d}{dr}\ln\psi^{\text{AE}}\right|_{r_c} $$

5.2 Ultrasoft Pseudopotentials (Vanderbilt)

Relaxes norm conservation for smoother potentials:

$$ \hat{H}|\psi_i\rangle = \epsilon_i\hat{S}|\psi_i\rangle $$

where:

$$ \hat{S} = 1 + \sum_{ij}q_{ij}|\beta_i\rangle\langle\beta_j| $$

5.3 Projector Augmented Wave (PAW) Method

Linear transformation connecting pseudo and all-electron wavefunctions:

$$

$$

Components:

• $|\tilde{\psi}\rangle$ — smooth pseudo-wavefunction
• $|\phi_i\rangle$ — all-electron partial waves
• $|\tilde{\phi}_i\rangle$ — pseudo partial waves
• $|\tilde{p}_i\rangle$ — projector functions

6. Brillouin Zone Integration

Physical observables require integration over $\mathbf{k}$-space:

$$ \langle A \rangle = \frac{1}{\Omega_{BZ}}\int_{BZ} A(\mathbf{k})\,d\mathbf{k} $$

6.1 Monkhorst-Pack Grid

Systematic $\mathbf{k}$-point sampling:

$$ \mathbf{k}_{n_1,n_2,n_3} = \sum_{i=1}^{3} \frac{2n_i - N_i - 1}{2N_i}\mathbf{b}_i $$

where:

• $n_i = 1, 2, \ldots, N_i$
• $\mathbf{b}_i$ are reciprocal lattice vectors
• Grid specified as $N_1 \times N_2 \times N_3$

6.2 Density of States

The tetrahedron method improves integration accuracy:

$$ g(E) = \frac{1}{\Omega_{BZ}}\int_{BZ}\delta(E - E_{n\mathbf{k}})\,d\mathbf{k} $$

Practical evaluation:

• Divide Brillouin zone into tetrahedra
• Linear interpolation of $E_n(\mathbf{k})$ within each tetrahedron
• Analytical integration of $\delta$-function

7. Self-Consistent Field (SCF) Iteration

7.1 Algorithm

1. Initialize density $n^{(0)}(\mathbf{r})$ 2. Construct $V_{\text{eff}}[n]$ 3. Diagonalize Kohn-Sham equations → obtain $\{\psi_i, \epsilon_i\}$ 4. Compute new density: $$ n^{\text{new}}(\mathbf{r}) = \sum_i^{\text{occ}}|\psi_i(\mathbf{r})|^2 $$ 5. Mix densities: $$ n^{\text{in}} = (1-\alpha)n^{\text{old}} + \alpha n^{\text{new}} $$ 6. Repeat until $\|n^{\text{new}} - n^{\text{old}}\| < \epsilon$

7.2 Mixing Schemes

• Linear mixing: Simple but slow convergence

$$ n^{(i+1)} = (1-\alpha)n^{(i)} + \alpha n^{\text{out},[i]} $$

• Pulay mixing (DIIS): Minimizes residual over history

$$ n^{\text{in}} = \sum_j c_j n^{(j)}, \quad \text{where } \{c_j\} \text{ minimize } \left\|\sum_j c_j R^{(j)}\right\| $$

• Broyden mixing: Quasi-Newton approach

$$ n^{(i+1)} = n^{(i)} - \alpha B^{(i)} R^{(i)} $$

8. Beyond DFT: The Band Gap Problem

DFT-LDA/GGA systematically underestimates band gaps.

Typical underestimation:

8.1 GW Approximation

The self-energy captures many-body corrections:

$$ \Sigma(\mathbf{r}, \mathbf{r}'; \omega) = \frac{i}{2\pi}\int G(\mathbf{r}, \mathbf{r}'; \omega+\omega')\,W(\mathbf{r}, \mathbf{r}'; \omega')\,d\omega' $$

Components:

• $G$ — single-particle Green's function
• $W$ — screened Coulomb interaction:

$$ W = \epsilon^{-1}v $$

Dielectric function (RPA):

$$ \epsilon(\mathbf{r}, \mathbf{r}'; \omega) = \delta(\mathbf{r} - \mathbf{r}') - \int v(\mathbf{r} - \mathbf{r}'')P^0(\mathbf{r}'', \mathbf{r}'; \omega)\,d\mathbf{r}'' $$

Quasiparticle correction:

$$ E_{n\mathbf{k}}^{\text{QP}} = E_{n\mathbf{k}}^{\text{DFT}} + \langle\psi_{n\mathbf{k}}|\Sigma(E^{\text{QP}}) - V_{xc}|\psi_{n\mathbf{k}}\rangle $$

This typically adds 0.5–2 eV to band gaps.

9. Effective Mass and k·p Theory

Near band extrema, expand energy to quadratic order:

$$ E_n(\mathbf{k}) \approx E_n(\mathbf{k}_0) + \frac{\hbar^2}{2}\sum_{ij}k_i\left(\frac{1}{m^*}\right)_{ij}k_j $$

9.1 Effective Mass Tensor

From second-order perturbation theory:

$$ \left(\frac{1}{m^*}\right)_{ij} = \frac{1}{m}\delta_{ij} + \frac{2}{m^2}\sum_{n' eq n}\frac{\langle n|\hat{p}_i|n'\rangle\langle n'|\hat{p}_j|n\rangle}{E_n - E_{n'}} $$

Alternate form using band curvature:

$$ \left(\frac{1}{m^*}\right)_{ij} = \frac{1}{\hbar^2}\frac{\partial^2 E_n}{\partial k_i \partial k_j} $$

9.2 8-Band Kane Model

For zincblende semiconductors (GaAs, InP, etc.):

$$ H_{\text{Kane}} = \begin{pmatrix} E_c + \frac{\hbar^2k^2}{2m_0} & \frac{P}{\sqrt{2}}k_+ & -\sqrt{\frac{2}{3}}Pk_z & \cdots \\ \frac{P}{\sqrt{2}}k_- & E_v - \frac{\hbar^2k^2}{2m_0} & \cdots & \cdots \\ \vdots & \vdots & \ddots & \vdots \end{pmatrix} $$

where:

• $k_\pm = k_x \pm ik_y$
• $P = \langle S|\hat{p}_x|X\rangle$ is the Kane momentum matrix element
• Includes: conduction band, heavy hole, light hole, split-off bands

10. Spin-Orbit Coupling

For heavier elements (Ge, GaAs, InSb):

$$ H_{\text{SO}} = \frac{\hbar}{4m^2c^2}( abla V \times \mathbf{p})\cdot\boldsymbol{\sigma} $$

10.1 Effects

• Lifts degeneracies: Valence band splitting ~0.34 eV in GaAs
• Essential for:
• Topological insulators
• Spintronics
• Optical selection rules

10.2 Matrix Form

The Hamiltonian becomes a $2 \times 2$ spinor structure:

$$ H = \begin{pmatrix} H_0 + H_{\text{SO}}^{zz} & H_{\text{SO}}^{+-} \\ H_{\text{SO}}^{-+} & H_0 - H_{\text{SO}}^{zz} \end{pmatrix} $$

where:

• $H_{\text{SO}}^{zz} = \lambda L_z S_z$
• $H_{\text{SO}}^{+-} = \lambda L_+ S_-$

11. Semiconductor Manufacturing Applications

11.1 Strain Engineering

Biaxial strain modifies band structure via deformation potentials:

$$ \Delta E_c = \Xi_d \cdot \text{Tr}(\boldsymbol{\epsilon}) + \Xi_u \cdot \epsilon_{zz} $$

Strain tensor components:

$$ \boldsymbol{\epsilon} = \begin{pmatrix} \epsilon_{xx} & \epsilon_{xy} & \epsilon_{xz} \\ \epsilon_{yx} & \epsilon_{yy} & \epsilon_{yz} \\ \epsilon_{zx} & \epsilon_{zy} & \epsilon_{zz} \end{pmatrix} $$

Valence band (Bir-Pikus Hamiltonian):

$$ H_{\epsilon} = a(\epsilon_{xx} + \epsilon_{yy} + \epsilon_{zz}) + 3b\left[(L_x^2 - \frac{1}{3}L^2)\epsilon_{xx} + \text{c.p.}\right] $$

Manufacturing application:

• Strained Si channels: ~30–50% mobility enhancement
• SiGe virtual substrates for strain control

11.2 Heterostructures and Quantum Wells

At interfaces, the envelope function approximation:

$$ \left[-\frac{\hbar^2}{2} abla\cdot\frac{1}{m^*(\mathbf{r})} abla + V(\mathbf{r})\right]F(\mathbf{r}) = EF(\mathbf{r}) $$

Ben Daniel-Duke boundary conditions:

$$ \begin{aligned} F_A(z_0) &= F_B(z_0) \\ \frac{1}{m_A^}\left.\frac{\partial F}{\partial z}\right|_A &= \frac{1}{m_B^}\left.\frac{\partial F}{\partial z}\right|_B \end{aligned} $$

Band alignment types:

• Type I (straddling): Both carriers confined in same layer (e.g., GaAs/AlGaAs)
• Type II (staggered): Electrons and holes in different layers (e.g., InAs/GaSb)
• Type III (broken gap): Conduction and valence bands overlap

11.3 Defects and Dopants

Supercell approach — create periodic array of defects.

Formation energy:

$$ E_f[D^q] = E_{\text{tot}}[D^q] - E_{\text{tot}}[\text{bulk}] - \sum_i n_i\mu_i + q(E_F + E_V + \Delta V) $$

where:

• $D^q$ — defect in charge state $q$
• $n_i$ — number of atoms of species $i$ added/removed
• $\mu_i$ — chemical potential of species $i$
• $E_F$ — Fermi level referenced to valence band maximum $E_V$
• $\Delta V$ — potential alignment correction

Charge transition levels:

$$ \epsilon(q/q') = \frac{E_f[D^q; E_F=0] - E_f[D^{q'}; E_F=0]}{q' - q} $$

Classification:

• Shallow donors/acceptors: $\epsilon$ near band edges
• Deep levels: $\epsilon$ in mid-gap (recombination centers)

11.4 Alloy Effects

Virtual Crystal Approximation (VCA):

$$ V_{\text{VCA}} = xV_A + (1-x)V_B $$

Bowing parameter:

$$ E_g(x) = xE_g^A + (1-x)E_g^B - bx(1-x) $$

Advanced methods:

• Coherent Potential Approximation (CPA) for disorder
• Special Quasirandom Structures (SQS) for explicit alloy supercells

12. Computational Complexity

12.1 Parallelization Strategies

• k-point parallelism: Different k-points on different processors
• Band parallelism: Different bands distributed across processors
• Real-space decomposition: Domain decomposition for large systems
• FFT parallelism: Distributed 3D FFTs for plane-wave methods

12.2 Key Software Packages

13. Workflow

┌─────────────────────────────────────────────────────────────┐
│                    INPUT: Crystal Structure                 │
│            (atomic positions, lattice vectors)              │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              SELECT METHOD                                  │
│   • DFT (LDA/GGA/Hybrid) for accuracy                       │
│   • Tight-binding for speed                                 │
│   • GW for accurate band gaps                               │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              COMPUTATIONAL SETUP                            │
│   • Choose k-point grid (Monkhorst-Pack)                    │
│   • Set energy cutoff (plane waves)                         │
│   • Select pseudopotentials                                 │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              SELF-CONSISTENT CALCULATION                    │
│   • Iterate until density converges                         │
│   • Obtain ground-state energy                              │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              POST-PROCESSING                                │
│   • Band structure along high-symmetry paths                │
│   • Density of states                                       │
│   • Effective masses                                        │
│   • Optical properties                                      │
└─────────────────────────────────────────────────────────────┘
                              │
                              ▼
┌─────────────────────────────────────────────────────────────┐
│              VALIDATION & APPLICATION                       │
│   • Compare with ARPES, optical data                        │
│   • Extract parameters for device simulation (TCAD)         │
└─────────────────────────────────────────────────────────────┘

14. Key Equations Reference Card

Schrödinger Equation $$ \hat{H}\psi = E\psi $$

Bloch Theorem $$ \psi_{n\mathbf{k}}(\mathbf{r}) = e^{i\mathbf{k}\cdot\mathbf{r}}u_{n\mathbf{k}}(\mathbf{r}) $$

Kohn-Sham Equation $$ \left[-\frac{\hbar^2}{2m} abla^2 + V_{\text{eff}}[n]\right]\psi_i = \epsilon_i\psi_i $$

Effective Mass $$ \frac{1}{m^*_{ij}} = \frac{1}{\hbar^2}\frac{\partial^2 E}{\partial k_i \partial k_j} $$

GW Self-Energy $$ \Sigma = iGW $$

Formation Energy $$ E_f = E_{\text{tot}}[\text{defect}] - E_{\text{tot}}[\text{bulk}] - \sum_i n_i\mu_i + qE_F $$