From the Schrödinger equation to a finite-difference matrix

The envelope-function equation is a differential equation; a computer solves an algebraic one. This page derives the finite-difference form of the position-dependent effective-mass operator, explains the interface averaging that keeps it Hermitian, and shows how the whole problem becomes a large sparse matrix eigenproblem \(H\Psi = E\Psi\) that the OghmaNano eigensolver handles.

1. Discretising the growth direction

The growth direction is sampled at \(N\) equally spaced nodes \(z_i = i\,\Delta z\), and the envelope is represented by its values \(\psi_i = \psi(z_i)\). Derivatives are replaced by finite differences. The subtlety is the kinetic term, which contains the mass between the derivatives,

\[ \hat{T}\,\psi = -\frac{d}{dz}\!\left[ C(z)\,\frac{d\psi}{dz} \right], \qquad C(z) = \frac{\hbar^2}{2 m^\*(z)} . \]

A naive centred difference of the outer derivative would evaluate \(C\) on the nodes, but the inner derivative naturally lives between nodes. The correct, current-conserving discretisation evaluates the coefficient at the half-nodes \(z_{i\pm 1/2}\):

\[ \hat{T}\,\psi\big|_i \approx \frac{1}{\Delta z^2}\Big[ (C_{i-1/2}+C_{i+1/2})\,\psi_i - C_{i-1/2}\,\psi_{i-1} - C_{i+1/2}\,\psi_{i+1} \Big]. \]

2. Why interface averaging is required

At a heterojunction the effective mass jumps abruptly. If the coefficient at a half-node were taken from only one side, the discrete operator would not be symmetric and probability current would not be conserved across the interface, producing spurious interface states and energies that drift with grid spacing. OghmaNano therefore averages the coefficient of adjacent nodes,

\[ C_{i+1/2} = \tfrac{1}{2}\big(C_i + C_{i+1}\big), \]

which makes the discrete kinetic matrix exactly symmetric (Hermitian) and is the discrete counterpart of the BenDaniel–Duke continuity condition. The same half-node averaging is used consistently for every position-dependent coefficient in the multiband solvers, including the Kane momentum term and the valence kinetic parameters. For the silicon Δ-valley solver the reciprocal mass \(1/m(z)\) is averaged, which is the natural quantity for that single-band operator.

3. From operator to matrix

Adding the diagonal potential and in-plane terms, the scalar problem becomes a symmetric tridiagonal matrix: each row couples a node only to itself and its two neighbours. The diagonal entry is

\[ H_{ii} = V(z_i) + \frac{\hbar^2 k_\parallel^2}{2 m^\*(z_i)} + \frac{C_{i-1/2}+C_{i+1/2}}{\Delta z^2}, \]

and the off-diagonal entries are \(H_{i,i\pm1} = -C_{i\pm1/2}/\Delta z^2\). Solving \(H\psi = E\psi\) yields the confined levels and their envelopes directly.

4. Multiband problems: block-sparse matrices

A multiband model carries several coupled envelope components at every node — eight for the zincblende and wurtzite k·p models, ten for the dilute-nitride model, six for Luttinger–Kohn. At each node the components form a small dense block set by the local Hamiltonian; between neighbouring nodes they are coupled by the discretised \(k_z = -i\,d/dz\) operators. The global matrix is therefore block tridiagonal: dense within a node, sparse between nodes. Its dimension is

\[ \dim H = N_{\text{grid}} \times N_{\text{basis}}, \]

so a 300-node grid with an 8-band model gives a \(2400 \times 2400\) complex Hermitian matrix. Because only a handful of eigenvalues near a target energy are wanted, this matrix is never formed densely; it is stored sparsely and solved with the shift-invert method described on the eigensolver page.