Linear Variation and the Secular Equations

Every electronic structure package in routine use — Gaussian, ORCA, Psi4, NWChem, Q-Chem — solves a secular equation of exactly this type at its core, just with $\mathbf H$ built from full electron–electron repulsion integrals (Hartree–Fock or Kohn–Sham operators) rather than empirical Hückel parameters, and with basis sets of hundreds or thousands of Gaussian functions rather than three $p_z$ orbitals. The Overlap and Hamiltonian Matrix Elements that populate $\mathbf S$ and $\mathbf H$ in these calculations are computed analytically for Gaussian basis functions, but the linear-algebra step that follows — Solving the Secular Determinant via matrix diagonalization — is the identical generalized eigenvalue problem derived here. Simple Hückel theory itself, despite its crudeness, remains in active pedagogical and even research use for rapid, back-of-envelope predictions about π-conjugated systems: estimating relative stabilities of polyene and polycyclic aromatic isomers, rationalizing regiochemistry in electrophilic aromatic substitution through frontier-orbital coefficients, and teaching the origin of aromaticity and antiaromaticity through the pattern of bonding/nonbonding/antibonding roots. Extended Hückel theory, a step up in sophistication that includes all valence orbitals (not just π electrons), is still used in solid-state and materials chemistry for qualitative band-structure reasoning in extended systems, again via the same secular-determinant construction scaled up to a much larger basis. In drug design and materials science, configuration interaction and multireference methods build their many-electron wavefunctions as linear combinations of Slater determinants — a direct many-electron analog of the linear combination of one-electron basis functions treated here — and diagonalize the resulting CI Hamiltonian matrix by the same secular-equation logic to obtain ground and excited electronic states relevant to photochemistry, chromophore design, and reaction mechanism studies. Even machine-learning-accelerated quantum chemistry, which increasingly predicts $\mathbf H$ or the electronic density directly, still typically returns to a secular-equation solve as the final step connecting a learned Hamiltonian to physically interpretable orbital energies and coefficients.

Ask the tutor: How do you actually evaluate the Hamiltonian and overlap matrix elements for a real Gaussian basis set before plugging them into the secular determinant?

Start from the variational principle in its most familiar form: for any normalizable trial function $\Psi_{\text{trial}}$,

where $E_0$ is the true ground-state energy. In the simple variational method you might pick a single functional form with one continuous parameter — an exponential $e^{-\alpha r}$, a Gaussian $e^{-\alpha r^2}$ — and minimize $E_{\text{trial}}(\alpha)$ over $\alpha$. The linear variation method generalizes this in a specific, extremely fruitful direction: instead of one function with a shape parameter, use a fixed set of $N$ basis functions $\phi_1, \phi_2, \dots, \phi_N$ (real, for simplicity, though the complex case is a direct generalization) and write

The $\phi_i$ are chosen for chemical or mathematical convenience — atomic orbitals centered on each atom in a π system, a set of Gaussian primitives, particle-in-a-box eigenfunctions, whatever spans a subspace believed to contain a good approximation to the true wavefunction. They do not need to be orthogonal or normalized in general, though many practical choices (Hückel theory's assumption of orthonormal 2p orbitals) simplify things by assuming so. The coefficients $c_i$ are now the variational parameters. There are $N$ of them, and the variational theorem still applies: whatever values of $c_1,\dots,c_N$ minimize the trial energy give the best possible energy and wavefunction within the space spanned by the $\phi_i$. This is the crucial restriction to keep in mind — the linear variation method never gives you the exact energy unless the true eigenfunction happens to lie exactly in your chosen subspace; it gives you the best answer achievable from the functions you decided to trust. A second, deeper feature emerges once you carry out the minimization (done in full in the derivations below): because $E_{\text{trial}}$ is a ratio of quadratic forms in the $c_i$, setting $\partial E/\partial c_k = 0$ for every $k$ produces not one equation but a whole system of $N$ coupled linear equations, one per coefficient, all sharing the same unknown energy $E$. This system is homogeneous — it has the trivial all-zero solution, which is physically useless (no wavefunction) — so a nontrivial solution exists only for special values of $E$ that make the determinant of the coefficient matrix vanish. Those special values are exactly the variational estimates of the energy levels, and there are $N$ of them for an $N$-function basis, ordered so that the lowest is the tightest upper bound on the true ground state, the next lowest bounds the true first excited state (subject to the basis containing enough structure to represent it), and so on — this ordering result is sometimes called MacDonald's theorem.

The bookkeeping that makes this concrete uses two families of integrals over the basis functions. The overlap integrals

measure how much basis functions $i$ and $j$ resemble each other in the sense of wavefunction overlap; $S_{ii} = 1$ if each $\phi_i$ is normalized, and $S_{ij} = 0$ for $i \neq j$ if the basis is orthogonal. The Hamiltonian matrix elements

encode the energetics: $H_{ii}$ is close to (though not exactly) the energy of an electron sitting purely in $\phi_i$, while $H_{ij}$ for $i \ne j$ is a coupling or "resonance" term describing how strongly $\phi_i$ and $\phi_j$ interact through the Hamiltonian. The detailed physical and computational meaning of these two matrices — how to evaluate them for real basis sets, what makes $H_{ij}$ large or small chemically — is developed at length elsewhere; here the point is only that once $S_{ij}$ and $H_{ij}$ are known numbers, the secular equations are pure linear algebra.

The trial function is

where:

• $\phi_i$ are $N$ fixed, linearly independent basis functions (not necessarily orthonormal),
• $c_i$ are real variational coefficients to be determined,
• $N$ is the size of the basis, fixed in advance by the modeler's choice.

Reading: this is not a single function to optimize but a whole family of functions parametrized by the point $(c_1,\dots,c_N)$ in an $N$-dimensional coefficient space; the variational method searches that space. The trial energy is the Rayleigh quotient

where:

• $H_{ij} = \int \phi_i \hat H \phi_j\, d\tau$ is the Hamiltonian matrix element,
• $S_{ij} = \int \phi_i \phi_j\, d\tau$ is the overlap matrix element,
• the sums run over all $i,j$ from 1 to $N$.

Reading: both numerator and denominator are quadratic forms in the coefficients; $E(\mathbf c)$ is bounded below by the true ground-state energy no matter what $\mathbf c$ is chosen, by the ordinary variational theorem. Minimizing $E(\mathbf c)$ with respect to every $c_k$ simultaneously (worked out below) yields the secular equations,

where:

• $E$ is now the specific trial energy at the stationary point, shared by all $N$ equations,
• each value of $k$ gives one linear equation in the unknowns $c_1,\dots,c_N$.

Reading: this is a homogeneous linear system; the only way it has a solution other than $c_1 = c_2 = \cdots = c_N = 0$ is if the determinant of its coefficient matrix is zero. That determinant condition is the secular determinant equation,

where:

• $\mathbf H$ is the $N\times N$ matrix with entries $H_{ij}$,
• $\mathbf S$ is the $N\times N$ overlap matrix with entries $S_{ij}$,
• the determinant, expanded, is a polynomial of degree $N$ in $E$.

Reading: solving this polynomial for $E$ gives $N$ roots $E_1 \le E_2 \le \cdots \le E_N$; each root, substituted back into the secular equations, yields a coefficient vector $\mathbf c^{(m)}$ and hence an approximate wavefunction $\Psi_m$. ##

derivation[From the secular equations to the secular determinant, and back to physical coefficients]

Goal. Show why the secular equations force $\det(\mathbf H - E\mathbf S) = 0$, and how, once a root $E_m$ is found, the corresponding coefficients $c_j^{(m)}$ are recovered and normalized.

Step 1 — Write the system in matrix form. The $N$ secular equations $\sum_j (H_{kj} - ES_{kj})c_j = 0$ for $k = 1,\dots,N$ are exactly the matrix equation

where $\mathbf c = (c_1, c_2, \dots, c_N)^{\mathsf T}$ is a column vector and $\mathbf 0$ is the zero vector. This is a homogeneous linear system with coefficient matrix $\mathbf M(E) \equiv \mathbf H - E\mathbf S$.

Step 2 — Invoke the theorem on homogeneous linear systems. A basic result of linear algebra states that a homogeneous system $\mathbf{M}\mathbf{c} = \mathbf 0$ has a solution other than $\mathbf c = \mathbf 0$ if and only if $\mathbf M$ is singular, i.e.

Physically we require a nonzero solution, since $\mathbf c = \mathbf 0$ means $\Psi \equiv 0$, not a wavefunction at all. So the energy $E$ appearing in a physically meaningful stationary point must satisfy

Step 3 — Expand the determinant and count roots. Written out, $\det(\mathbf H - E\mathbf S)$ is a sum of products of $N$ matrix entries (by the Leibniz expansion of an $N\times N$ determinant), each entry linear in $E$. The highest power of $E$ that can appear comes from the product of all $N$ diagonal-type entries $(H_{ii} - E S_{ii})$, giving an $E^N$ term. Hence

is a genuine degree-$N$ polynomial in $E$ (assuming $\mathbf S$ nonsingular, which holds whenever the basis functions are linearly independent). A degree-$N$ polynomial has exactly $N$ roots (real, in this Hermitian case, and in fact it can be shown they are all real because $\mathbf H$ and $\mathbf S$ are Hermitian and $\mathbf S$ is positive definite). Label them in increasing order $E_1 \le E_2 \le \cdots \le E_N$.

Step 4 — Recover coefficients for each root. For each root $E_m$, substitute $E = E_m$ back into the secular equations,

Because $\det(\mathbf H - E_m \mathbf S) = 0$ by construction, this system has rank at most $N-1$, so it determines the ratios $c_1^{(m)} : c_2^{(m)} : \cdots : c_N^{(m)}$ uniquely (up to an overall scale) but not their absolute magnitudes — exactly as many independent equations as needed to fix $N-1$ ratios, with one equation redundant.

Step 5 — Fix the scale by normalization. The overall scale is fixed by requiring $\langle \Psi_m|\Psi_m\rangle = 1$, i.e.

This is one additional algebraic condition (not a new physical equation, just bookkeeping) that pins down the magnitude of $\mathbf c^{(m)}$, completing the determination of the approximate eigenfunction $\Psi_m = \sum_i c_i^{(m)}\phi_i$ belonging to variational energy $E_m$.

Step 6 — Interpret multiplicity of roots. Each of the $N$ roots $E_1,\dots,E_N$, together with its eigenvector, is a distinct approximate stationary state of the system restricted to the chosen basis — not just the ground state. MacDonald's theorem guarantees that if the exact Hamiltonian has true eigenvalues $\varepsilon_0 \le \varepsilon_1 \le \varepsilon_2 \le \cdots$, then the linear variation roots satisfy $E_m \ge \varepsilon_{m-1}$ for every $m = 1, \dots, N$ (using 1-based indexing for $E_m$ and 0-based for $\varepsilon$), so every root is a genuine upper bound to the corresponding true level, not merely the lowest one.

Example 1 — Foundational: two non-orthogonal basis functions, done by hand

Consider a schematic two-basis-function problem (the numbers are chosen to be asymmetric and dimensionless, standing in for energies in units of some reference, e.g. eV) with

Notice the basis is not orthogonal ($S_{12}\ne0$) and the two diagonal energies are unequal — this is deliberately asymmetric so that no accidental symmetry disguises an error.

Step 1: write the secular determinant.

i.e.

Step 2: expand. First term:

Second term:

Subtracting,

Step 3: solve the quadratic. Using the quadratic formula,

Compute the discriminant:

So

The two roots are

Check ordering: $E_1 = -8.003 < H_{11}=-8.0 < H_{22} = -5.0 < E_2 = -4.634$. This is exactly the expected qualitative behavior of a two-level mixing problem — the lower root is pushed below the lower diagonal element and the upper root is pushed above the higher diagonal element, because off-diagonal coupling always splits levels further apart (level repulsion).

Step 4: get coefficients for $E_1 = -8.003$. Use the first secular equation, $(H_{11}-E_1)c_1 + (H_{12}-E_1S_{12})c_2 = 0$:

(The near-zero coefficient of $c_1$ reflects rounding sensitivity right at the numerically nearly-degenerate root; carrying one more decimal in $E_1$, e.g. $E_1=-8.0033$, gives a stable ratio $c_1/c_2 \approx 3.1$, consistent with the lower root being dominated by, but not purely equal to, $\phi_1$ since $H_{11}$ is the more stable diagonal energy.) Normalizing with $S_{12}=0.30$ via $c_1^2 + c_2^2 + 2c_1c_2S_{12}=1$ fixes the overall scale, completing $\Psi_1$.

Limiting case check. If the off-diagonal coupling is switched off, $H_{12}\to 0$ and $S_{12}\to 0$, the secular determinant becomes simply $(H_{11}-E)(H_{22}-E)=0$, whose roots are $E = H_{11}$ and $E = H_{22}$ exactly — the two basis functions decouple into two independent, unmixed one-function problems, as they must.

Example 2 — Applied: allyl radical π system by simple Hückel theory

The allyl system (three contiguous $p_z$ orbitals, as in the allyl radical or allyl cation/anion π framework) is the smallest chemically realistic case with $N=3$, and it is solved using the two standard Hückel Molecular Orbital Theory approximations: all $S_{ii}=1$, all $S_{ij}=0$ for $i\ne j$ (orthonormal basis of $2p_z$ AOs), all $H_{ii} = \alpha$ (same Coulomb integral on every carbon), and $H_{ij} = \beta$ only between bonded neighbors (here $\beta <0$), zero otherwise. With this orthonormal basis, $\mathbf S = \mathbf 1$ and the secular equations reduce to the ordinary eigenvalue problem $\mathbf{Hc} = E\mathbf c$. The Hamiltonian matrix is

Step 1: substitute $x = (\alpha-E)/\beta$, a standard Hückel change of variable (dividing every entry by $\beta$, a negative number, and defining $x$ this way keeps the algebra in clean integers). The determinant condition $\det(\mathbf H - E\mathbf 1)=0$ becomes

Step 2: expand along the first row.

Step 3: roots.

Converting back via $E = \alpha - x\beta$ (remember $\beta<0$, so a negative $x$ gives an energy below $\alpha$, i.e. more bonding):

Numerically, using the standard reference values $\alpha = -11.2$ eV and $\beta = -2.7$ eV (typical carbon 2p Hückel parameters),

Wait — since $\beta<0$, check signs carefully: $E_1=\alpha-\sqrt2\beta$ with $\beta$ negative means $-\sqrt2\beta$ is positive, so $E_1 = -11.2+3.818=-7.38$ eV is the least stable (highest) of the bonding-labeled root, and $E_3=\alpha+\sqrt2\beta=-11.2-3.818=-15.02$ eV is in fact the lowest, most stable orbital. Relabeling by actual energy order: the three π molecular orbital energies are $-15.02$ eV (bonding, fully in-phase), $-11.2$ eV (nonbonding, node through the central carbon), and $-7.38$ eV (antibonding), which is exactly the qualitative allyl MO diagram taught in every organic mechanism course: one bonding, one nonbonding, one antibonding π MO, populated by 3 electrons in the radical (2 in bonding, 1 in nonbonding) or 2 electrons in the cation (2 in bonding only) or 4 electrons in the anion (2 bonding + 2 nonbonding).

Step 4: coefficients for the lowest root, $x=-\sqrt2$ i.e. $E_3=\alpha+\sqrt2\beta$. Substituting into row 1 of $(\mathbf H - E\mathbf 1)\mathbf c=0$ with $x=-\sqrt2$:

Row 2: $c_1 - \sqrt2 c_2 + c_3 = 0 \implies c_1 + c_3 = \sqrt2 c_2 = 2c_1 \implies c_3 = c_1$. So the (unnormalized) lowest MO has coefficients in ratio $c_1:c_2:c_3 = 1:\sqrt2:1$. Normalizing with $\sum c_i^2=1$ (orthonormal basis):

This matches the textbook allyl bonding MO $\psi_1 = 0.500\,\phi_1 + 0.707\,\phi_2 + 0.500\,\phi_3$ — largest amplitude on the central carbon, equal and same-sign amplitude on the two termini, with no nodes, exactly the fully bonding combination expected as the lowest-energy π MO.

Example 3 — Exam/competition level: three-function basis with unequal overlaps, extracting the ground-state energy bound

A common "trap" style question gives a basis that is neither fully orthonormal nor uses simplified Hückel parameters, forcing full use of the generalized eigenvalue machinery. Take $N=2$ again (to keep the arithmetic checkable exactly by hand) but with deliberately unequal, non-Hückel numbers, in eV:

Step 1: secular determinant.

Step 2: expand the first product.

Step 3: expand the square.

Step 4: subtract and collect.

Divide through by 0.96 for convenience:

Step 5: quadratic formula.

Step 6: sanity-check against the diagonal elements. $E_1 = -14.323$ eV is (correctly) more negative than $H_{11}=-14.30$ eV, and $E_2 = -11.677$ eV is (correctly) less negative than $H_{22}=-11.90$ eV — the coupling has pushed the lower root down and the upper root up, exactly the expected level-repulsion behavior seen already in Example 1.

Step 7: coefficients for the ground-state root $E_1=-14.323$ eV. From row 1: $(H_{11}-E_1S_{11})c_1 + (H_{12}-E_1S_{12})c_2 = 0$:

So the ground-state trial function is heavily dominated by $\phi_1$ (consistent with $\phi_1$'s much lower diagonal energy $-14.30$ eV compared to $\phi_2$'s $-11.90$ eV), with only a small admixture of $\phi_2$. Normalizing via $c_1^2 + c_2^2 + 2c_1c_2S_{12} = 1$ with $c_1 = 10.24\,c_2$:

Limiting case check. If instead $H_{11}=H_{22}$ (equal diagonal energies, the symmetric special case deliberately avoided above so as not to mask an error), the two roots collapse to $E_{1,2} = (H_{11}\pm H_{12})/(1\pm S_{12})$ and the coefficients become exactly $c_1=\pm c_2$ (in-phase and out-of-phase combinations) — the textbook symmetric two-orbital splitting pattern, recoverable as a special case of the general asymmetric formula used throughout this example.

The linear variation method as used throughout quantum chemistry owes its structural form to the general variational principle developed across the early quantum-mechanical literature (Rayleigh's classical work on vibrating systems, later adapted quantum mechanically), but its decisive chemical application — turning a determinant of Coulomb and resonance integrals into predictions about aromatic stability, bond orders, and reactivity — is due to Erich Hückel, who in the early 1930s applied exactly this secular-determinant machinery to the π-electron systems of unsaturated and aromatic hydrocarbons.