Sparse matrices, shift-invert and Rayleigh–Ritz

Discretising any of the quantum-well models produces the same mathematical object: a large, sparse, complex Hermitian matrix whose eigenvalues are the confined energies and whose eigenvectors are the envelopes. Only a handful of states near the band edge are physically relevant, so a full diagonalisation would be wasteful. OghmaNano instead uses a shift-and-invert block iteration that converges rapidly onto exactly the states of interest. This page explains the mathematics; it is deliberately not a description of the linear-algebra library calls.

1. The eigenproblem

The task is to find eigenpairs \((E_n, \Psi_n)\) of

\[ H\,\Psi_n = E_n\,\Psi_n, \]

where \(H\) has dimension \(N_{\text{grid}}\times N_{\text{basis}}\) (see the discretisation page). \(H\) is Hermitian, so its eigenvalues are real and its eigenvectors can be chosen orthonormal. It is also sparse — block tridiagonal — so matrix–vector products and factorisations are far cheaper than for a dense matrix of the same size.

2. Why only near-edge states matter

A \(2400\times2400\) Hamiltonian has 2400 eigenvalues, but a laser or detector uses only the lowest few conduction subbands and the highest few valence subbands. Computing all 2400 would cost orders of magnitude more effort for no benefit. The challenge is that the wanted states are in the interior of the spectrum (the valence-band maximum, the conduction-band minimum), not at its extremes, where simple iterative methods converge. Shift-and-invert solves precisely this interior-eigenvalue problem.

3. The shift-and-invert transformation

Choose a target energy \(\sigma\) near the states of interest and consider the transformed operator \((H-\sigma I)^{-1}\). If \(\Psi\) is an eigenvector of \(H\) with eigenvalue \(E\), then

\[ (H - \sigma I)^{-1}\,\Psi = \frac{1}{E - \sigma}\,\Psi. \]

The eigenvectors are unchanged, but the eigenvalues are mapped to \(1/(E-\sigma)\). States with \(E\) close to \(\sigma\) are mapped to very large magnitudes, while distant states are mapped to small ones. An iteration that amplifies the largest-magnitude eigenvalues — ordinary power or inverse iteration — therefore converges onto exactly the states nearest \(\sigma\). This is the heart of the method.

4. Inverse iteration in practice

Applying \((H-\sigma I)^{-1}\) does not require inverting the matrix. OghmaNano forms the shifted matrix \(A = H - \sigma I\) once and computes its sparse LU factorisation. Each application of the transformed operator is then a forward/back substitution — solving \(A\,y = q\) for \(y\) — which is fast because the factors inherit the sparsity of \(A\). A single factorisation is reused for every iteration and every vector.

5. Block iteration and orthonormalisation

Iterating a single vector finds only one state and struggles when states are nearly degenerate (as quantum-well subbands often are). OghmaNano therefore iterates a whole block of vectors at once — twice as many as the number of states requested — so that clustered levels are resolved together. After each application of \((H-\sigma I)^{-1}\) the block is re-orthonormalised by modified Gram–Schmidt (with two passes for numerical robustness), which prevents all the vectors from collapsing onto the single dominant direction and keeps the block spanning the wanted subspace.

6. Rayleigh–Ritz projection

The orthonormalised block \(Q\) spans a good approximation to the wanted subspace, but its columns are not yet eigenvectors. The Rayleigh–Ritz procedure extracts the best eigenpairs from the subspace. It projects \(H\) onto the block,

\[ T = Q^\dagger H\,Q, \]

a small dense matrix of the block size, and diagonalises it. Its eigenvalues \(\theta_i\) are the Ritz values (approximate energies) and the corresponding combinations \(x_i = Q\,c_i\) are the Ritz vectors (approximate envelopes). Because \(T\) is small, this dense diagonalisation is cheap, and it optimally combines the block vectors into the best available eigenvector estimates.

7. Residual and convergence

The quality of each Ritz pair is measured by its residual norm,

\[ r_i = \big\lVert H\,x_i - \theta_i\,x_i \big\rVert, \]

which is zero for an exact eigenpair. OghmaNano selects the requested number of Ritz values nearest \(\sigma\) and declares convergence when all their residuals fall below a tolerance (typically \(10^{-6}\)–\(10^{-8}\) eV). If not converged, the Ritz vectors reseed the block for another inverse iteration, and the process repeats. This combination — shift-invert to target the interior, block iteration for clustered states, and Rayleigh–Ritz to extract the eigenpairs — is a robust and efficient interior eigensolver used identically by all the multiband models.