MO Energy Levels and the Secular Equation

Every quantum chemistry software package — from a teaching-level extended Hückel calculation to a production-grade density functional theory code — ultimately solves a generalized matrix eigenvalue problem of the form $\mathbf{HC}=\mathbf{SCE}$, the direct linear-algebra descendant of the secular determinant derived here; the Roothaan–Hall equations of Hartree–Fock theory are exactly this equation dressed up with a self-consistent Fock matrix in place of a fixed Hamiltonian. Photoelectron spectroscopists interpret ionization energies as approximations to the negative of orbital energies (Koopmans' theorem), meaning every peak in a UV photoelectron spectrum of, say, formaldehyde or ethylene is effectively a measured root of some underlying secular equation. Materials scientists building tight-binding models of graphene, carbon nanotubes, and topological insulators use secular determinants of exactly this Hückel-style form, with the famous linear ($E\propto |k|$) dispersion of graphene's Dirac cones falling directly out of a $2\times2$ secular determinant evaluated at each crystal momentum. Drug designers reasoning about the frontier molecular orbitals (HOMO/LUMO) of a candidate inhibitor — where nucleophilic or electrophilic attack will occur, whether a chromophore will absorb in the visible — are reading off the highest and lowest occupied roots of a secular equation computed (usually) at a semi-empirical or DFT level, but conceptually identical to the two-orbital and three-orbital determinants worked out above. Spectroscopic assignment of vibronic and electronic transitions in conjugated dyes, from indigo to cyanine laser dyes, still begins with a Hückel-style secular treatment of the $\pi$ system to get qualitatively correct orbital orderings before more expensive methods refine the numbers.

Ask the tutor: Can you walk through how the secular determinant changes shape once we build in the zero-overlap approximation used in Hückel theory, and why chemists accept that trade-off?

Picture two atomic orbitals, $\phi_1$ centered on atom A and $\phi_2$ centered on atom B, each a perfectly good one-electron wavefunction on its own atom. When the atoms approach, an electron near this diatomic framework should be described by some wavefunction $\psi$ that is not exactly $\phi_1$ or $\phi_2$ but has to smoothly interpolate and superpose both. The LCAO Approximation proposes the simplest possible ansatz:

The coefficients $c_1$ and $c_2$ are not known in advance — they are the very things we are solving for. The strategy is to treat them as variational parameters: compute the energy expectation value of $\psi$ as a function of $c_1, c_2$, then find the values of $c_1, c_2$ that minimize (or, more precisely, make stationary) that energy. Quantum mechanics guarantees, via the Variational Principle, that any such stationary energy is an upper bound to a true eigenvalue of the Hamiltonian, and that the best two-parameter trial function of this form gives the two energies closest to the true bonding and antibonding orbital energies achievable within this basis. The physical content hiding inside the resulting algebra is this: the atomic orbitals $\phi_1, \phi_2$ are not orthogonal in general (they overlap in space, quantified by $S = \langle\phi_1|\phi_2\rangle$), and the Hamiltonian couples them through the resonance integral $\beta = H_{12} = \langle\phi_1|\hat H|\phi_2\rangle$, a negative number for bonding interactions between orbitals of the same phase. The diagonal Hamiltonian element $\alpha = H_{11} \approx H_{22}$ (equal for homonuclear cases) is the Coulomb integral, essentially the energy an electron would have sitting in that atomic orbital in the molecular environment. All three integrals — $\alpha$, $\beta$, $S$ — were characterized on their own terms elsewhere; here they become the raw material fed into a determinant.

Reading the diagram: the two atomic levels, degenerate at energy $\alpha$, are pulled apart into a lower bonding combination and a higher antibonding combination. The magnitude of splitting and its asymmetry — the antibonding level rises more than the bonding level falls — is a direct, unavoidable consequence of nonzero overlap $S$, derived rigorously below. The secular determinant is the formal device that finds exactly these two roots and no others, for any number of basis orbitals, not only two.

This second figure is the actual algebraic pipeline: differentiate the variational energy expression, obtain a pair of linear homogeneous equations in $c_1, c_2$, and recognize that these equations have a nonzero solution only when their coefficient determinant vanishes. That last step is not a physical assumption — it is a theorem of linear algebra, and it is the entire content of "the secular equation."

For a general $n$-orbital LCAO basis $\{\phi_1,\ldots,\phi_n\}$, define the Hamiltonian and overlap matrix elements

where: $H_{ij}$ is the resonance integral (or Coulomb integral when $i=j$) between basis orbitals $i,j$; $S_{ij}$ is their overlap integral, with $S_{ii}=1$ for normalized orbitals. Both matrices are real and symmetric here because $\hat H$ is Hermitian and the basis is chosen real, as is standard for atomic $p$, $s$ orbitals in simple MO treatments. The condition that the variationally optimized LCAO coefficients exist and are nontrivial is the secular equation:

where: $\mathbf{H}$ is the $n\times n$ matrix of $H_{ij}$; $\mathbf{S}$ is the $n\times n$ overlap matrix; $E$ is the unknown MO energy, playing the role of the generalized eigenvalue. This is a polynomial equation of degree $n$ in $E$; it has $n$ real roots (counting multiplicity) because $\mathbf{H}$ and $\mathbf{S}$ are real symmetric with $\mathbf{S}$ positive definite — these are exactly the $n$ MO energy levels built from this basis. For each root $E_k$, the coefficient vector $\mathbf{c}^{(k)}$ solving

where: $\mathbf{c}^{(k)}$ is the (unnormalized) column of LCAO coefficients for the $k$-th molecular orbital; this equation is solved after $E_k$ is known, one root at a time — gives the shape of that particular molecular orbital once normalized against $\langle\psi_k|\psi_k\rangle=1$. Specializing to the two-orbital homonuclear case ($H_{11}=H_{22}=\alpha$, $H_{12}=H_{21}=\beta$, $S_{11}=S_{22}=1$, $S_{12}=S_{21}=S$), the determinant reads

where: $\alpha$ is the atomic Coulomb integral (energy of an electron confined to one atomic orbital in the molecular field); $\beta$ is the resonance integral coupling the two orbitals; $S$ is their mutual overlap; $E$ is the sought MO energy. Expanding this $2\times2$ determinant and solving the resulting quadratic gives the two roots

where: $E_+$ (using top signs, since $\beta<0$ conventionally) is the lower, bonding orbital energy $\dfrac{\alpha+\beta}{1+S}$; $E_-$ is the higher, antibonding orbital energy $\dfrac{\alpha-\beta}{1-S}$. Reading this equation: because $S>0$ typically and $\beta<0$, the denominator $(1-S)$ on the antibonding root is smaller than $(1+S)$ on the bonding root, which is precisely why the antibonding level rises further above $\alpha$ than the bonding level falls below it — an asymmetry with real consequences for whether a four-electron interaction (both orbitals doubly filled, as in He$_2$) is net bonding or antibonding.

Example 1 — Foundational: homonuclear diatomic with explicit, non-round numbers. Take a model diatomic (parameters chosen to be numerically distinct and realistic in magnitude, not tuned to round out prettily) with

Compute the bonding and antibonding energies using $E_\pm = (\alpha\pm\beta)/(1\pm S)$. Bonding root:

Antibonding root:

So the bonding MO sits $\Delta_+ = \alpha - E_+ = -11.4-(-11.890) = 0.490$ eV below the atomic level, while the antibonding MO sits $\Delta_- = E_- - \alpha = -10.548-(-11.4) = 0.852$ eV above it. Checking the algebraic prediction from the second derivation, $\Delta_-/\Delta_+$ should equal $(1+S)/(1-S) = 1.27/0.73 = 1.740$; numerically $0.852/0.490 = 1.739$, matching (small rounding only). This single calculation demonstrates every qualitative feature of MO splitting: the average of $E_+,E_-$ is not $\alpha$ (it is pulled slightly, because of the asymmetric denominators, though for small $S$ it is close), and the antibonding destabilization visibly exceeds the bonding stabilization. Limiting case check: as $S\to 0$, both denominators collapse to 1, giving $E_\pm \to \alpha\pm\beta$ symmetric about $\alpha$ — this is exactly the zero-overlap approximation used in elementary MO cartoons and is recovered correctly from the general formula, confirming internal consistency.

Example 2 — Applied: heteronuclear bond (unequal Coulomb integrals), e.g. an HF-like $\sigma$ system. Now the two atomic orbitals are not equivalent: $H_{11} = \alpha_A = -13.6$ eV (a hydrogen-like 1s), $H_{22} = \alpha_B = -18.7$ eV (a more electronegative fluorine-like 2p$_z$), with resonance integral $\beta = H_{12} = -2.6$ eV and overlap $S=0.20$. The secular determinant is now

Expanding (careful: the two diagonal terms are no longer equal, so no simple difference-of-squares factoring applies):

Expand each piece. First, $(\alpha_A-E)(\alpha_B-E) = \alpha_A\alpha_B - E(\alpha_A+\alpha_B) + E^2$. Numerically, $\alpha_A\alpha_B = (-13.6)(-18.7) = 254.32$, and $\alpha_A+\alpha_B = -32.3$, so this piece is $254.32 + 32.3E + E^2$. Second, $(\beta-ES)^2 = \beta^2 - 2\beta S E + S^2E^2 = (-2.6)^2 - 2(-2.6)(0.20)E + (0.20)^2E^2 = 6.76 + 1.04E + 0.04E^2$. Subtracting,

Applying the quadratic formula, $E = \dfrac{-31.26 \pm \sqrt{31.26^2 - 4(0.96)(247.56)}}{2(0.96)}$. Compute the discriminant: $31.26^2 = 977.19$, and $4(0.96)(247.56) = 950.63$, so the discriminant is $977.19 - 950.63 = 26.56$, whose square root is $5.154$.

Bonding root: $E_+ = \dfrac{-31.26+5.154}{1.92} = \dfrac{-26.106}{1.92} = -13.597$ eV. Wait — check ordering: the more negative root is the lower (bonding) energy, so take the branch giving the more negative value. Antibonding root: $E_- = \dfrac{-31.26-5.154}{1.92} = \dfrac{-36.414}{1.92} = -18.966$ eV. Re-examine: $-18.966$ eV is lower than $-13.597$ eV, so in fact $E_+=-18.966$ eV (bonding, lower) and $E_-=-13.597$ eV (antibonding, higher). This makes physical sense: since $\alpha_B=-18.7$ eV is already lower than $\alpha_A=-13.6$ eV, the bonding MO — built mostly from the lower-energy fluorine-like orbital — sits close to but below $\alpha_B$, while the antibonding MO sits close to but above $\alpha_A$. The bonding MO is stabilized by only $18.966-18.7=0.266$ eV relative to $\alpha_B$, much less than the homonuclear Example 1's stabilization, because unequal atomic energies reduce mixing efficiency — a general and important trend: the more electronegativity difference between two atoms, the less the bonding MO is stabilized relative to the lower atomic level, even though the bond may still be strong in absolute terms because that level was already low. Limiting case check: if $\alpha_A\to\alpha_B$ this heteronuclear formula must reduce to the homonuclear $E_\pm=(\alpha\pm\beta)/(1\pm S)$ result; setting $\alpha_A=\alpha_B=\alpha$ in the quadratic $(\alpha-E)^2-(\beta-ES)^2=0$ reproduces exactly the Example 1 factoring, confirming the general heteronuclear treatment correctly contains the homonuclear case as a special instance.

Example 3 — Competition level: three-orbital linear chain (a minimal Hückel-style secular determinant). Consider three $p_z$ orbitals in a linear (non-cyclic) chain — a toy allyl-radical-like system — with all Coulomb integrals equal to $\alpha$, nearest-neighbor resonance integral $\beta$ between orbitals 1–2 and 2–3 (but not 1–3, since they are not adjacent), and — to make the example a genuine test of the general determinant rather than the familiar zero-overlap Hückel shortcut — retain a small nearest-neighbor overlap $S=0.15$ while setting overlap between non-adjacent orbitals 1 and 3 to exactly zero (a reasonable approximation, since they are far apart). The secular determinant is

Let $x = \alpha - E$ and $y=\beta - ES$ for brevity. The determinant expands (along the first row) as

This factors immediately into $x=0$ or $x^2 = 2y^2$, i.e. $x = \pm y\sqrt2$.

Root 1: $x=0 \Rightarrow \alpha - E = 0 \Rightarrow E = \alpha$. This is the nonbonding molecular orbital, exactly at the atomic energy — the middle rung of the allyl-type ladder, independent of both $\beta$ and $S$, a structurally important and testable prediction.

Root 2: $x = y\sqrt2 \Rightarrow \alpha - E = \sqrt2(\beta - ES) \Rightarrow \alpha - E = \sqrt2\beta - \sqrt2 ES$

Root 3: $x=-y\sqrt2 \Rightarrow \alpha - E = -\sqrt2\beta+\sqrt2 ES \Rightarrow \alpha+\sqrt2\beta = E(1+\sqrt2S) \Rightarrow E = \frac{\alpha+\sqrt2\beta}{1+\sqrt2 S}$ Now insert numbers: $\alpha = -10.5$ eV, $\beta=-2.9$ eV, $S=0.15$, so $\sqrt2\beta = -4.101$ eV and $\sqrt2 S = 0.2121$. Root 2 (using $-\sqrt2\beta = +4.101$ in numerator): $E = \dfrac{-10.5+4.101}{1-0.2121} = \dfrac{-6.399}{0.7879} = -8.122$ eV. This is the antibonding MO (highest energy of the three, since it's the least negative). Root 3: $E = \dfrac{-10.5-4.101}{1+0.2121} = \dfrac{-14.601}{1.2121} = -12.045$ eV. This is the bonding MO (most negative, lowest energy). Root 1 remains $E=\alpha=-10.5$ eV, the nonbonding MO, sitting numerically between the other two ($-12.045 < -10.5 < -8.122$), exactly as the qualitative allyl MO diagram requires: one bonding level below $\alpha$, one nonbonding level at $\alpha$, one antibonding level above $\alpha$. Limiting case check: setting $S\to 0$ recovers the standard zero-overlap Hückel allyl result $E = \alpha, \alpha\pm\sqrt2\beta$ directly, since both denominators $(1\pm\sqrt2 S)$ collapse to 1 — this is precisely the simplification that Hückel theory adopts as a standing approximation, and this worked example shows explicitly what is being neglected when that approximation is made.

The secular determinant as a chemical tool is inseparable from Erich Hückel, who in 1931 applied exactly this variational-plus-determinant machinery to $\pi$ electrons in unsaturated and aromatic hydrocarbons, reducing the electronic structure problem to a small, hand-solvable polynomial equation — a breakthrough precisely because it made quantum mechanical reasoning about aromaticity and reactivity tractable without a computer. Hückel trained originally as a physicist under Peter Debye, working on dipole moments and electrolyte theory (the Debye–Hückel theory of ionic solutions bears his name for entirely separate reasons), before turning his attention in the late 1920s to the puzzle of why benzene's six $\pi$ electrons behave so differently from three isolated double bonds. His molecular orbital treatment, built from a secular determinant with all diagonal Coulomb integrals equal and only nearest-neighbor resonance integrals nonzero, gave the now-iconic result that a closed shell of $4n+2$ $\pi$ electrons in a cyclic conjugated system is unusually stable — a rule chemists still apply by name today. The mathematical core of that theory is exactly the vanishing-determinant condition developed in this article, specialized to cyclic and acyclic conjugated $\pi$ systems with the zero-overlap simplification.