Solving the Secular Determinant

Every electronic-structure calculation performed today — density functional theory, Hartree–Fock, configuration interaction, coupled cluster — ultimately reduces, at some stage, to diagonalizing a Hamiltonian-like matrix in a finite basis, which is precisely the secular-determinant problem scaled up to hundreds or thousands of basis functions and solved numerically rather than by hand. Drug-design software computing frontier molecular orbital energies to predict reactivity, materials scientists computing band structures of periodic solids (where the secular determinant becomes a function of crystal momentum $\mathbf{k}$), and spectroscopists fitting vibronic coupling models to absorption spectra are all, underneath the software, solving $\det(\mathbf{H}(\mathbf{k})-E\mathbf{S})=0$ or its close relatives. In inorganic and organometallic chemistry, ligand-field and crystal-field splitting diagrams — the $t_{2g}/e_g$ gap in octahedral complexes, the basis of color and magnetism in transition-metal compounds — are secular-determinant results for a five-orbital d-basis interacting with a ligand symmetry field. In photochemistry and spectroscopy, the "avoided crossing" of two electronic potential energy surfaces near a conical intersection is a $2\times2$ secular problem at every nuclear geometry, with the off-diagonal coupling element controlling how sharply the surfaces repel — directly the algebra of the asymmetric $2\times2$ derivation above. Even simple qualitative reasoning tools such as Walsh diagrams and perturbational molecular orbital theory in organic mechanism drawing are secular-determinant intuitions applied without writing down the matrix explicitly. In materials science, tight-binding models of graphene, carbon nanotubes, and topological insulators are Hückel-style secular determinants with a periodic basis, whose roots reproduce the famous linear (Dirac-cone) dispersion near the Fermi level.

Ask the tutor: How would you set up and solve the secular determinant for cyclobutadiene's four π orbitals, and what does the degeneracy of its roots tell you about its instability?

Picture two atomic orbitals, $\phi_1$ on atom A and $\phi_2$ on atom B, each an eigenfunction of its own isolated atomic Hamiltonian with energy roughly $H_{11}$ and $H_{22}$. Bring the atoms together. The true molecular wavefunction is not $\phi_1$ or $\phi_2$ alone but some blend $\psi = c_1\phi_1 + c_2\phi_2$. The variational method says: try this blend, compute the energy expectation value as a function of $c_1,c_2$, and minimize. The minimization conditions are two linear, homogeneous equations in $c_1$ and $c_2$. Homogeneous linear systems only have nonzero solutions when their determinant vanishes — otherwise the only solution is the trivial $c_1=c_2=0$, which describes no electron at all. Demanding a nontrivial mixture is exactly what forces the energy to take one of two special values, $E_+$ and $E_-$. Physically, $E_-$ is the bonding combination — the coefficients $c_1,c_2$ have the same sign, constructive overlap builds up electron density between the nuclei, and the energy drops below both $H_{11}$ and $H_{22}$. $E_+$ is antibonding — opposite signs, a node between the nuclei, energy pushed above both atomic levels. The secular determinant is the algebraic engine that produces this splitting quantitatively, including the unequal case where atom A and atom B have different electronegativities ($H_{11} \neq H_{22}$), which is the realistic situation for any heteronuclear bond, C–N π bonds in imines, or a metal d-orbital interacting with a ligand lone pair.

The same mathematics scales directly to $n$ basis functions: an $n\times n$ secular determinant produces $n$ roots. For conjugated π-systems this is the Hückel method — each carbon 2p orbital is a basis function, the determinant is built from Coulomb integrals $\alpha$ on the diagonal and resonance integrals $\beta$ (only between bonded neighbors) off the diagonal, and the roots are the π molecular orbital energies.

The starting point is the variational expectation value for a trial function that is a linear combination of $n$ real basis functions,

where: $c_i$ are (real, for this treatment) mixing coefficients to be determined, and $\phi_i$ are fixed basis functions (atomic orbitals, Gaussians, etc.) that are not necessarily orthogonal or normalized. Reading: this is simply "build the unknown wavefunction as a weighted sum of known functions" — the coefficients carry all the physics still to be determined. The variational energy functional is

where: $H_{ij}$ are Hamiltonian matrix elements and $S_{ij}$ are overlap matrix elements, both computable once the basis is chosen — this is exactly the content of Overlap and Hamiltonian Matrix Elements. Reading: $E(\mathbf{c})$ is a ratio of quadratic forms in the coefficients — the Rayleigh quotient — and the variational theorem guarantees $E(\mathbf{c}) \ge E_0$ (the true ground-state energy) for any choice of coefficients. Minimizing $E$ with respect to each $c_k$ (treating $\mathbf{c}$ as a free real vector, using the quotient rule) and setting $\partial E/\partial c_k = 0$ for every $k$ gives the matrix eigenvalue equation

where: this is a generalized eigenvalue problem; $E$ here has become the stationary value of the Rayleigh quotient, not an independent variable. Reading: at a stationary point, the vector $\mathbf{Hc}$ is parallel to $\mathbf{Sc}$ with proportionality constant $E$ — the same structure as an ordinary eigenvalue equation, generalized to a non-orthogonal basis. Rearranging as a homogeneous linear system $(\mathbf{H}-E\mathbf{S})\mathbf{c} = \mathbf{0}$, linear algebra says a nonzero $\mathbf{c}$ exists only if the coefficient matrix is singular:

where: this is the secular equation (equivalently, the secular determinant set to zero); it is a polynomial equation of degree $n$ in $E$. Reading: this single condition is the gateway from "an infinite continuum of possible $E$" to "exactly $n$ physically allowed energies" — everything downstream (Hückel levels, ligand-field splittings, CI state energies) is this equation solved for a specific $\mathbf{H}$ and $\mathbf{S}$. Once a root $E_k$ is found, the corresponding coefficient vector $\mathbf{c}^{(k)}$ is obtained from

where: the second equation is the normalization condition fixing the overall scale of $\mathbf{c}^{(k)}$ (and sign, up to an arbitrary overall $\pm$). Reading: the determinant gives you where the energy levels are; this final linear solve (plus normalization) gives you what the orbital looks like at each level — the physically meaningful mixing ratios.

Example 1 — Foundational: a genuinely asymmetric heteronuclear two-orbital problem. Take two atomic orbitals with unequal Coulomb integrals, deliberately non-round numbers and nonzero overlap, mimicking a polarized σ bond:

Using the quadratic-formula result from the derivation, first compute the pieces:

Next the constant term:

And $1-S^2 = 1-0.0324 = 0.9676$. The discriminant is

So the two roots are

Explicitly,

Check against Step 3 of the derivation: the lower root $-18.32$ eV lies below both $H_{11}=-14.7$ and $H_{22}=-18.3$ eV (just barely below the more stable atomic level, as expected for weak-to-moderate coupling), and the upper root $-14.63$ eV lies above both — the bracketing theorem holds. For coefficients, take $E_-=-18.32$ eV:

A ratio this large means the bonding MO here is overwhelmingly $\phi_2$ (the lower-energy, more electronegative orbital) with only a small admixture of $\phi_1$ — exactly the qualitative expectation for a polarized bond: the bonding MO resembles the more stable AO.

Example 2 — Applied: Hückel allyl system (three contiguous 2p orbitals, C1–C2–C3). The Hückel secular determinant, using $\alpha$ for every Coulomb integral (identical carbons), $\beta$ for adjacent resonance integrals, zero for non-adjacent pairs, and $S=0$ off-diagonal (Hückel's usual orthogonality assumption), in the reduced variable $x=(\alpha-E)/\beta$:

Expanding along the first row:

Factor: $x(x^2-2)=0$, giving $x=0,\ \pm\sqrt{2}$. Converting back via $E=\alpha - x\beta$ (note $\beta<0$ by convention, so larger $x$ means lower energy):

With typical carbon values $\alpha=-11.2$ eV, $\beta=-2.7$ eV (asymmetric, chemically realistic numbers, not the bare Hückel symbols):

The allyl radical has 3 π electrons: 2 fill $E_1$, 1 occupies $E_2$ (the nonbonding, purely radical-character orbital split evenly across C1 and C3 with a node at C2 — the classic allylic delocalization result). Total π energy is $2E_1+E_2$, and the delocalization stabilization relative to a hypothetical localized single double bond ($2(\alpha+\beta)$ for one π bond plus one nonbonding electron at $\alpha$) is

a genuine, quantitatively meaningful resonance stabilization, negative (stabilizing) because $\beta<0$.

Example 3 — Exam/competition level: perturbing a symmetric 2×2 case toward asymmetry and testing the weak-coupling limit. Consider two orbitals with a small energy mismatch $\Delta = H_{22}-H_{11}$, zero overlap ($S=0$, the Hückel-style simplification), and coupling $H_{12}=\beta$. Take specific asymmetric numbers: $H_{11}=-10.0$ eV, $H_{22}=-10.9$ eV (so $\Delta=-0.9$ eV), $\beta = -0.35$ eV — a case where the mismatch $|\Delta|$ is comparable to, but larger than, $|\beta|$, testing whether the exact formula and a perturbative estimate agree. With $S=0$ the quadratic simplifies (set $S=0$ in the boxed formula from the derivation):

Numerically, the average is $(-10.0-10.9)/2=-10.45$ eV, and the half-difference is $(-10.0-(-10.9))/2=0.45$ eV, so

Now compare to the weak-coupling (perturbative) limit, where $|\beta|\ll|\Delta|$ predicts the shift is approximately $\beta^2/\Delta$ added/subtracted from the unperturbed levels — foreshadowing exactly the structure that reappears in First-Order Perturbation Theory as second-order energy correction. The unperturbed levels are $-10.0$ and $-10.9$ eV; the perturbative correction magnitude is

Perturbative estimate: lower level (originally $-10.9$) pushed down to $-11.036$ eV; upper level (originally $-10.0$) pushed up to $-9.864$ eV. Compare to the exact secular roots, $-11.020$ and $-9.880$ eV: the perturbative estimate is close (within $0.02$ eV) but not identical, because $|\beta/\Delta| = 0.39$ is not truly small — this is exactly the regime where the exact secular-determinant solution and the perturbative shortcut start to visibly diverge, illustrating precisely when the full quadratic must be trusted over the leading-order approximation.

The secular determinant is the meeting point of two traditions: celestial mechanics, where "secular" perturbations (slow, long-term orbital drifts) were first analyzed via eigenvalue problems in the 18th and 19th centuries by Lagrange and Laplace, and quantum chemistry, where Erich Hückel adapted exactly this linear-algebra machinery to π-electron systems in 1931. Hückel, a German physicist trained under Peter Debye, was trying to explain the peculiar chemical stability of benzene and other aromatic compounds — a problem organic chemists had puzzled over since Kekulé's structural proposals. Hückel's radical simplification was to treat the π-electrons of a conjugated system as moving in an effective field set by only nearest-neighbor interactions, reducing the many-electron problem to a set of one-electron secular equations with just two empirical parameters, $\alpha$ and $\beta$. The Hückel method's astonishing longevity — it remains a teaching staple and a first qualitative check nearly a century later — owes everything to the transparency of the secular-determinant approach: any organic chemist with pencil and paper (or, later, a slide rule, and now a laptop) could set up the determinant for a new conjugated system, expand the polynomial, and read off bonding, nonbonding, and antibonding levels, along with aromatic stabilization energies that correlated remarkably well with reactivity and spectroscopy. The same determinant-and-eigenvector logic, generalized to non-orthogonal bases and many more basis functions, became the backbone of the Roothaan–Hall equations for Hartree–Fock method-based molecular orbital theory decades later.