Overlap and Hamiltonian Matrix Elements

Every quantum chemistry calculation eventually comes down to filling two matrices with numbers. Pick a basis of functions to describe an electron's possible whereabouts — atomic orbitals, Gaussians, plane waves, whatever suits the problem. The entire physical content of the trial wavefunction then collapses into two arrays. First, an overlap matrix $S$ that tells you how the basis functions crowd into each other's space. Second, a Hamiltonian matrix $H$ that tells you how much energy each pairing costs. Diagonalizing the resulting secular determinant then hands you the approximate energies and orbitals. But before any determinant can be solved, someone has to build $H$ and $S$ correctly. Understanding what each entry means physically — not just how to compute it — matters greatly. That understanding is what separates a chemist who can read a Hückel or extended-Hückel printout from one who is just following a recipe. This matters because $H_{ij}$ and $S_{ij}$ are not abstract bookkeeping devices; they encode chemistry. A large off-diagonal $H_{ij}$ between two atomic orbitals signals a strong bonding or antibonding interaction. It is, essentially, a quantitative resonance integral. It is the direct ancestor of the " $\beta$ " that shows up in every Hückel diagram a student ever draws. A nonzero $S_{ij}$ signals that two basis functions are not orthogonal — they physically overlap in space. This single fact is responsible for some of the least intuitive features of quantum chemical energy expressions. It explains why increasing basis-set overlap can sometimes raise rather than lower an orbital's apparent energy contribution. This holds unless the algebra is handled with full rigor. The big picture is this: any linear variational trial function $\Psi = \sum_i c_i \phi_i$ is built from non-orthogonal basis functions $\phi_i$ . Via the variational principle, this leads to a generalized eigenvalue problem $Hc = ESc$ . This is different from the simpler $Hc = Ec$ you would get from an orthonormal basis. The matrices $H$ and $S$ are what stand between a chosen basis set and a usable energy spectrum. Get their definitions, symmetries, and physical interpretations right, and the secular determinant becomes almost self-explanatory. Get them wrong, and every downstream number is compromised — bond order, orbital energy, HOMO–LUMO gap. This article builds $H_{ij}$ and $S_{ij}$ from first principles. It works through the bookkeeping of Hermiticity and normalization. It connects the abstract integrals to concrete Hückel-style numbers. It shows exactly how each entry earns its place in the matrix. That matrix is what Solving the Secular Determinant will eventually diagonalize.

Key takeaways

$S_{ij} = \langle \phi_i \vert \phi_j\rangle$ measures the spatial overlap of basis functions $i$ and $j$ ; $S_{ii}=1$ for normalized basis functions.
$H_{ij} = \langle \phi_i \vert \hat H \vert \phi_j\rangle$ measures the energetic coupling between basis functions $i$ and $j$ through the true Hamiltonian $\hat H$ .
Both matrices are Hermitian ( $H_{ij}=H_{ji}^*$ , $S_{ij}=S_{ji}^*$ ); for real basis functions this means simply symmetric.
$H_{ii}$ is a Coulomb-type integral (the energy of an electron confined to orbital $i$ alone); $H_{ij}$ ( $i\neq j$ ) is a resonance/exchange-type integral driving bonding.
The generalized eigenvalue equation $Hc=ESc$ , not the simpler $Hc=Ec$ , is the correct secular equation whenever $S_{ij}\neq 0$ for $i\neq j$ .
Hückel theory is the famous special case where $S_{ij}$ is approximated as the identity matrix even though the true atomic-orbital overlaps are not zero.
Key idea: $H_{ij}$ and $S_{ij}$ are simply the matrix elements of the identity and Hamiltonian operators in a chosen (possibly non-orthogonal) basis, and together they convert an intractable differential eigenvalue problem into a solvable matrix one.

Where this lives today

Every modern electronic-structure package — Gaussian, ORCA, Psi4, NWChem, VASP, CP2K — builds and diagonalizes generalized eigenvalue problems of exactly the form $Hc=ESc$ (there called $FC=SC\varepsilon$ in Roothaan–Hall notation, with the Fock matrix $F$ playing the role of $H$ ) at every iteration of a self-consistent field calculation. The overlap matrix $S$ in these codes is computed from real, explicit Gaussian-type-orbital integrals, and is never treated as the identity except when the code performs an internal Löwdin or Gram–Schmidt orthogonalization step specifically to simplify diagonalization. The physical interpretation carried over from this article — diagonal elements as self-energies, off-diagonal elements as couplings — remains exactly the lens through which chemists read Mulliken population analyses, crystal-orbital Hamilton population (COHP) plots in solid-state chemistry, and tight-binding Hamiltonians in materials science and graphene/nanotube band-structure calculations. In drug design, semi-empirical and tight-binding QM/MM methods (used to model enzyme active sites, for instance in modeling proton or hydride transfer in dehydrogenases) rely on exactly this $H_{ij}$ / $S_{ij}$ machinery, often with parametrized $H_{ij}$ (as in extended Hückel or DFTB, density-functional tight-binding) standing in for expensive ab initio integrals to make large biomolecular systems tractable. In materials science, the tight-binding models used to predict the electronic structure of graphene, transition-metal dichalcogenide monolayers, and topological insulators are structurally identical Hückel-style secular problems, with $H_{ij}$ hopping integrals between lattice sites replacing the carbon–carbon $\beta$ of organic $\pi$ systems, and $S_{ij}$ often (again) approximated as the identity for simplicity, exactly as in the original 1931 treatment. Spectroscopists extracting electronic coupling strengths $H_{12}$ in dimeric photosynthetic pigments or organic semiconductors from measured exciton splittings use precisely the overlap-corrected two-level formula derived in Example 3 of this article, since neglecting $S_{12}$ there would systematically bias the extracted electronic coupling — a quantity central to predicting charge-transfer and energy-transfer rates via Marcus theory and Förster/Dexter mechanisms.

Ask the tutor: Once I have built $H$ and $S$ for my basis, what is the actual step-by-step procedure for extracting the energies and coefficients from $\det(H-ES)=0$ ?

The physical picture

Start with the Linear Variation and the Secular Equations setup: choose $N$ basis functions $\phi_1,\phi_2,\dots,\phi_N$ — atomic orbitals on different atoms, Gaussian primitives, or any convenient set — and write the trial wavefunction as a linear combination,

\Psi = \sum_{i=1}^{N} c_i \phi_i .

The variational theorem says that the "best" $\Psi$ within this family minimizes the Rayleigh quotient $E[\Psi]=\langle\Psi|\hat H|\Psi\rangle/\langle\Psi|\Psi\rangle$ . Expanding both integrals in the basis functions immediately produces sums of two kinds of integral: $\langle\phi_i|\hat H|\phi_j\rangle$ in the numerator, and $\langle\phi_i|\phi_j\rangle$ in the denominator. These are exactly $H_{ij}$ and $S_{ij}$ . Physically, $S_{ij}$ is a geometric quantity — it depends only on the shapes and positions of the basis functions, never on the Hamiltonian. If $\phi_i$ and $\phi_j$ are atomic orbitals centered on atoms that are far apart, their tails barely reach each other and $S_{ij}\to 0$ . If the atoms are pulled closer, the orbitals interpenetrate and $S_{ij}$ grows, approaching $1$ only in the unphysical limit that the two functions become identical (same center, same orbital). Overlap is always between $-1$ and $+1$ for real, normalized functions, a direct consequence of the Cauchy–Schwarz inequality $|S_{ij}| \le \sqrt{S_{ii}S_{jj}} = 1$ . $H_{ij}$ , by contrast, is an energetic quantity. It asks: if the system's true Hamiltonian $\hat H$ acts on basis function $\phi_j$ , how much of the resulting function looks like $\phi_i$ , weighted by energy? The diagonal elements $H_{ii}=\langle\phi_i|\hat H|\phi_i\rangle$ are the expectation value of the energy if the electron were (falsely) assumed to sit entirely in orbital $\phi_i$ — these are called Coulomb integrals in Hückel language, usually denoted $\alpha_i$ , and are large negative numbers (bound-state energies) for atomic-orbital bases. The off-diagonal elements $H_{ij}$ for $i\neq j$ are resonance (or exchange) integrals, denoted $\beta_{ij}$ in Hückel theory — they quantify how strongly the true Hamiltonian mixes basis function $i$ into basis function $j$ 's domain, and it is precisely this mixing that lowers the energy of bonding combinations and raises the energy of antibonding ones. The key structural fact tying $H$ and $S$ together is Hermiticity. Because $\hat H$ is a Hermitian operator (its eigenvalues, the physical energies, must be real), the matrix built from it obeys

H_{ij} = H_{ji}^{*} .

The identity operator (implicit in the overlap integral) is trivially Hermitian too, so

S_{ij} = S_{ji}^{*} .

For real basis functions — the overwhelmingly common case in introductory quantum chemistry — complex conjugation does nothing, and both statements simplify to ordinary symmetry: $H_{ij}=H_{ji}$ and $S_{ij}=S_{ji}$ . This halves the work of building the matrices, since only the upper (or lower) triangle needs to be computed independently.

Figure

Two basis orbitals with overlap

S_{ij}

and coupling

H_{ij}

The picture above is deliberately drawn with a real nonzero overlap region because that shaded lens is the geometric meaning of $S_{ij}$ : it is literally the integral of the product $\phi_i\phi_j$ over all space, largest where the two orbitals' amplitudes are simultaneously large. $H_{ij}$ cannot be drawn as a simple region because the Hamiltonian operator (kinetic energy plus potential energy) reshapes $\phi_j$ before the overlap with $\phi_i$ is taken — it is an energetic overlap, not a spatial one, though for real bonding orbitals the two integrals do tend to have correlated signs (larger spatial overlap generally produces a more negative, more stabilizing $H_{ij}$ , though this correlation is not a strict theorem).

Figure

From basis functions to the secular matrices

The core ideas

The defining integrals are

S_{ij} = \int \phi_i^{*}(\mathbf{r})\, \phi_j(\mathbf{r}) \, d\mathbf{r} = \langle \phi_i \vert \phi_j \rangle ,

where: $\phi_i,\phi_j$ are members of the chosen basis set, and the integral runs over all electronic coordinates the basis depends on.

Reading: $S_{ij}$ is the projection of $\phi_j$ onto $\phi_i$ — a pure measure of how much the two functions coincide in space, with no reference to energetics.

H_{ij} = \int \phi_i^{*}(\mathbf{r})\, \hat H \, \phi_j(\mathbf{r}) \, d\mathbf{r} = \langle \phi_i \vert \hat H \vert \phi_j \rangle ,

where: $\hat H$ is the full (or effective, in semi-empirical theories) Hamiltonian operator for the system.

Reading: $H_{ij}$ is the amplitude for the Hamiltonian to connect basis states $i$ and $j$ — the diagonal terms are self-energies, the off-diagonal terms are couplings.

H_{ij} = H_{ji}^{*}, \qquad S_{ij} = S_{ji}^{*} ,

where: the star denotes complex conjugation, required for real eigenvalues/real overlaps.

Reading: Both matrices must be Hermitian because $\hat H$ is Hermitian and the identity ("overlap") operator is trivially Hermitian; for real-valued basis functions both matrices are simply symmetric.

S_{ii} = 1 \quad (\text{if } \phi_i \text{ normalized}), \qquad |S_{ij}| \le 1 \ \ (i\neq j),

where: normalization fixes each diagonal overlap to unity, and the Cauchy–Schwarz inequality bounds every off-diagonal overlap.

Reading: The overlap matrix behaves like a "correlation matrix" between basis functions — its diagonal is always $1$ , its off-diagonals measure the degree of non-orthogonality.

Hc = E\,Sc ,

where: $c$ is the column vector of linear-variation coefficients and $E$ is the variational energy being sought (an eigenvalue of this generalized problem).

Reading: This is the matrix form of the secular equations — it reduces to the ordinary eigenvalue problem $Hc=Ec$ only when $S$ is the identity matrix, i.e. only for an orthonormal basis.

\det(H - E S) = 0 ,

where: this determinant condition is the requirement that $Hc=ESc$ have a nontrivial solution $c\neq 0$ .

Reading: Expanding this determinant produces the secular polynomial in $E$ whose roots are the approximate energies — the object that Solving the Secular Determinant takes as its starting point.

Derivation AFrom the variational integral to the matrix elements

Goal. Show explicitly how minimizing $E[\Psi]=\langle\Psi|\hat H|\Psi\rangle/\langle\Psi|\Psi\rangle$ over the coefficients $\{c_i\}$ of a linear trial function forces $H_{ij}$ and $S_{ij}$ to appear, and derive the resulting matrix secular equation.

Step 1 — Write the trial function and substitute. Let $\Psi=\sum_{i=1}^N c_i\phi_i$ with real coefficients $c_i$ and (for simplicity of notation) real basis functions $\phi_i$ . Then

\langle \Psi \vert \hat H \vert \Psi \rangle = \sum_{i=1}^N \sum_{j=1}^N c_i c_j \langle \phi_i \vert \hat H \vert \phi_j \rangle = \sum_{i,j} c_i c_j H_{ij},

\langle \Psi \vert \Psi \rangle = \sum_{i=1}^N \sum_{j=1}^N c_i c_j \langle \phi_i \vert \phi_j \rangle = \sum_{i,j} c_i c_j S_{ij}.

These double sums are exact consequences of linearity of the integral — no approximation has been made yet. Every cross term in the expansion of $\Psi^2$ (or $\Psi^*\Psi$ for complex functions) becomes one matrix element.

Step 2 — Minimize with respect to each coefficient. Write $E=E[\{c_k\}]=N(\mathbf c)/D(\mathbf c)$ with $N=\sum_{ij}c_ic_jH_{ij}$ and $D=\sum_{ij}c_ic_jS_{ij}$ . The variational condition is $\partial E/\partial c_k = 0$ for every $k=1,\dots,N$ . Using the quotient rule,

\frac{\partial E}{\partial c_k} = \frac{1}{D}\frac{\partial N}{\partial c_k} - \frac{N}{D^2}\frac{\partial D}{\partial c_k} = 0 \quad\Longrightarrow\quad \frac{\partial N}{\partial c_k} = E\,\frac{\partial D}{\partial c_k},

where the last step used $E=N/D$ to replace $N/D$ .

Step 3 — Evaluate the derivatives. Because $H_{ij}=H_{ji}$ (Hermiticity, real case), differentiating the double sum $N=\sum_{ij}c_ic_jH_{ij}$ with respect to a single coefficient $c_k$ picks up two equal contributions (one from each index touching $k$ ):

\frac{\partial N}{\partial c_k} = 2\sum_{j=1}^{N} c_j H_{kj}, \qquad \frac{\partial D}{\partial c_k} = 2\sum_{j=1}^{N} c_j S_{kj}.

Step 4 — Assemble the secular equations. Substituting into Step 2 and cancelling the common factor of 2 gives, for every $k=1,\dots,N$ ,

\sum_{j=1}^{N} c_j H_{kj} = E \sum_{j=1}^{N} c_j S_{kj} \quad\Longleftrightarrow\quad \sum_{j=1}^N \left(H_{kj}-E S_{kj}\right) c_j = 0 .

This is exactly $N$ linear homogeneous equations in the $N$ unknowns $c_j$ , one equation for each choice of $k$ . Written as a single matrix equation, this is $Hc=ESc$ , or equivalently $(H-ES)\mathbf c=0$ .

Step 5 — Nontrivial-solution condition. A homogeneous linear system $(H-ES)\mathbf c = 0$ has a solution other than the trivial $\mathbf c=0$ only if the coefficient matrix is singular, i.e.

\det (H - ES) = 0 .

This is the secular determinant. Its $N$ roots $E_1\le E_2\le\cdots\le E_N$ are the variational estimates of the true energy levels accessible within the chosen basis, ordered from lowest (best estimate of the ground state, by the variational theorem) to highest.

Go deeper · Derivation B · 1 minHermiticity of

H_{ij}

and the reality of the eigenvalues

Derivation BHermiticity of

H_{ij}

and the reality of the eigenvalues

Goal. Prove that $H_{ij}=H_{ji}^{*}$ follows directly from the Hermitian character of $\hat H$ , and show why this guarantees real eigenvalues $E$ out of the secular determinant — a property that must hold for any physically sensible energy matrix.

Step 1 — State the defining property of a Hermitian operator. An operator $\hat H$ is Hermitian if, for any two well-behaved functions $f,g$ in its domain,

\int f^{*}\,(\hat H g)\, d\mathbf r = \int (\hat H f)^{*}\, g \, d\mathbf r .

This is the operator-level statement of self-adjointness; physically it is required so that $\hat H$ 's eigenvalues (observable energies) come out real.

Step 2 — Apply it to two basis functions. Set $f=\phi_i$ and $g=\phi_j$ . The left side is $\int \phi_i^{*}\hat H\phi_j\,d\mathbf r = H_{ij}$ by definition. The right side is $\int(\hat H\phi_i)^{*}\phi_j\,d\mathbf r$ , which is the complex conjugate of $\int\phi_i^*\hat H\phi_j^{*\ast}\,d\mathbf{r}$ ... more directly, note that

\int (\hat H \phi_i)^{*} \phi_j\, d\mathbf r = \left[ \int \phi_i^{*} \hat H \phi_j^{}\,d\mathbf r \Big|_{i\leftrightarrow j,\ \text{conjugated}} \right] = \left( \int \phi_j^{*} \hat H \phi_i \, d\mathbf r \right)^{*} = H_{ji}^{*}.

Step 3 — Conclude the symmetry relation. Equating the two sides from the Hermiticity condition gives

H_{ij} = H_{ji}^{*}.

For a basis of real functions (the standard case for Hückel-type real atomic $p$ -orbitals or real Gaussians), complex conjugation is inert, so this collapses to the ordinary matrix symmetry

H_{ij} = H_{ji}.

The identical argument applied to the identity operator (whose "Hamiltonian" role in the overlap integral is trivial, $\hat 1 f = f$ ) gives $S_{ij}=S_{ji}^{*}$ , or $S_{ij}=S_{ji}$ for real bases.

Step 4 — Why this forces real eigenvalues. Consider the generalized eigenvalue equation $Hc=ESc$ with $H,S$ both Hermitian and $S$ additionally positive definite (true whenever the basis functions are linearly independent). Multiply on the left by $c^{\dagger}$ :

c^{\dagger} H c = E\, c^{\dagger} S c .

Both $c^{\dagger}Hc$ and $c^{\dagger}Sc$ are real numbers — this is a standard property of quadratic forms built from Hermitian matrices, since $(c^\dagger Hc)^* = c^\dagger H^\dagger c = c^\dagger Hc$ using $H^\dagger=H$ . Since $S$ is positive definite, $c^{\dagger}Sc>0$ for $c\neq 0$ , so

E = \frac{c^{\dagger}Hc}{c^{\dagger}Sc}

is manifestly a ratio of two real numbers, hence real. Every root of the secular determinant is therefore guaranteed to be a real energy — never a complex number — precisely because both $H$ and $S$ were built to be Hermitian from the start.

Worked examples

Example 1 — Foundational: building $H$ and $S$ for two non-orthogonal 1s-like functions. Take a minimal two-function basis, $\phi_1$ and $\phi_2$ , representing valence orbitals on two different atoms in a diatomic-like fragment, with the following (asymmetric, deliberately non-round) parameters obtained from an assumed Hamiltonian and overlap model:

H_{11} = -13.6\ \text{eV}, \quad H_{22} = -17.4\ \text{eV}, \quad H_{12}=H_{21} = -3.10\ \text{eV}, \quad S_{12}=S_{21}=0.18 .

Both basis functions are normalized, so $S_{11}=S_{22}=1$ . The secular determinant is

\det\begin{pmatrix} H_{11}-E & H_{12}-ES_{12} \\ H_{12}-ES_{12} & H_{22}-E \end{pmatrix} = 0 .

Expanding,

(H_{11}-E)(H_{22}-E) - (H_{12}-ES_{12})^2 = 0 .

Substitute numbers. First the diagonal product:

(-13.6-E)(-17.4-E) = E^2 + 31.0E + 236.64 .

Next the off-diagonal term:

(H_{12}-ES_{12})^2 = (-3.10 - 0.18E)^2 = 0.0324E^2 + 1.116E + 9.61 .

Subtracting,

E^2 + 31.0E + 236.64 - 0.0324E^2 - 1.116E - 9.61 = 0,

0.9676\,E^2 + 29.884\,E + 227.03 = 0 .

Divide through by $0.9676$ :

E^2 + 30.887\,E + 234.68 = 0 .

Apply the quadratic formula:

E = \frac{-30.887 \pm \sqrt{30.887^2 - 4(234.68)}}{2} = \frac{-30.887 \pm \sqrt{954.0 - 938.7}}{2} .

E = \frac{-30.887 \pm \sqrt{15.3}}{2} = \frac{-30.887 \pm 3.912}{2} .

This gives two roots:

E_{-} = \frac{-30.887 - 3.912}{2} = -17.40\ \text{eV}, \qquad E_{+} = \frac{-30.887 + 3.912}{2} = -13.49\ \text{eV}.

Both energies lie close to, but shifted from, the original diagonal values $-13.6$ and $-17.4\ \text{eV}$ — exactly the qualitative behavior expected: nonzero $H_{12}$ and $S_{12}$ split and shift the two "atomic" levels into a bonding-like combination near $-17.4$ eV and an antibonding-like combination near $-13.5$ eV, with the lower root pushed to slightly below the lowest input diagonal energy, as the variational theorem requires. Limiting-case check: if $S_{12}\to 0$ (orthogonal basis) the same determinant reduces to $(H_{11}-E)(H_{22}-E)-H_{12}^2=0$ , the ordinary two-level Hückel-type formula with no overlap correction. Setting $S_{12}=0$ in this example, keeping $H_{12}=-3.10$ eV, gives $E^2+31.0E+236.64-9.61=0 \Rightarrow E^2+31.0E+227.03=0$ , whose roots are $E=(-31.0\pm\sqrt{961-908.1})/2=(-31.0\pm 7.27)/2$ , i.e. $-19.14$ eV and $-11.86$ eV — visibly different from the $S_{12}=0.18$ result above, confirming that overlap materially changes both roots and is not a negligible refinement.

Example 2 — Applied: the allyl radical with explicit overlap correction beyond simple Hückel. Simple Hückel theory for the allyl system ( $\text{CH}_2\!=\!\text{CH}\!-\!\text{CH}_2^{\bullet}$ , three contiguous $2p_z$ orbitals) sets $H_{11}=H_{22}=H_{33}=\alpha$ , $H_{12}=H_{23}=\beta$ , $H_{13}=0$ , and — crucially — approximates $S_{ij}=\delta_{ij}$ (identity overlap matrix), even though the true carbon $2p_z$ – $2p_z$ overlap for adjacent, $\pi$ -bonded carbons is typically around $S\approx 0.25$ , not zero. This example shows how including the real nearest-neighbor overlap changes the result. Keep $H_{13}=0$ (no direct 1–3 coupling, standard Hückel topology) and $S_{13}=0$ (atoms 1 and 3 are too far apart for meaningful direct overlap), but now include $S_{12}=S_{23}=0.25$ explicitly, alongside $H_{12}=H_{23}=\beta$ and $H_{11}=H_{22}=H_{33}=\alpha$ . The generalized secular determinant is

\det\begin{pmatrix} \alpha-E & \beta-ES_{12} & -E S_{13}\\ \beta - ES_{12} & \alpha - E & \beta - ES_{23} \\ -ES_{13} & \beta-ES_{23} & \alpha - E \end{pmatrix} = 0, \qquad S_{13}=0 .

Substitute $x = (\alpha - E)$ for compactness and write $\beta' = \beta - ES_{12}$ (the same value appears twice by symmetry). The determinant expands (by cofactor expansion along the first row) as

x\left(x^2 - \beta'^2\right) - \beta'\left(\beta' x - 0\right) + 0 = x^3 - x\beta'^2 - \beta'^2 x = x^3 - 2\beta'^2 x = x\left(x^2 - 2\beta'^2\right) = 0.

So the roots satisfy either $x=0$ (i.e. $E=\alpha$ ) or $x^2=2\beta'^2$ , i.e. $x=\pm\sqrt2\,\beta'$ . Take numerical values typical of a carbon $2p_z$ basis: $\alpha = -11.2$ eV, $\beta = -2.7$ eV (a slightly stronger resonance integral than the textbook $-2.5$ eV default, chosen here to keep the numbers asymmetric and non-round), and $S_{12}=0.25$ . The middle root is exact and overlap-independent: $E_2 = \alpha = -11.2$ eV (the nonbonding molecular orbital, unaffected by nearest-neighbor overlap because its coefficient vanishes on the central atom by symmetry). For the outer roots, write $x = \alpha - E$ and $\beta' = \beta - ES_{12}$ , and solve $x^2 = 2\beta'^2$ self-consistently (since $\beta'$ itself depends on $E$ ). Take $x = -\sqrt2\,\beta'$ for the bonding root (lower energy, more negative $x$ means $E>\alpha$ ... actually check signs carefully): since $\beta<0$ and we expect the bonding combination to lie below $\alpha$ (more negative $E$ , i.e. $x=\alpha - E>0$ ), take the branch $x = \sqrt2\,|\beta'|$ with $\beta'<0$ , i.e. $x = -\sqrt2\,\beta'$ .

\alpha - E = -\sqrt2 (\beta - ES_{12}) \quad\Longrightarrow\quad \alpha - E = -\sqrt2\,\beta + \sqrt2\,S_{12}\,E .

Collect $E$ terms:

\alpha + \sqrt2\,\beta = E + \sqrt2\,S_{12}E = E(1+\sqrt2\,S_{12}) .

E = \frac{\alpha + \sqrt2\,\beta}{1+\sqrt2\,S_{12}} .

Plug in numbers: $\sqrt2 = 1.41421$ , so $\sqrt2\beta = 1.41421\times(-2.7) = -3.818$ eV, and $\sqrt2 S_{12} = 1.41421\times0.25=0.3536$ .

E_{\text{bonding}} = \frac{-11.2 + (-3.818)}{1+0.3536} = \frac{-15.018}{1.3536} = -11.096\ \text{...}

Recompute carefully: $-15.018/1.3536 = -11.096$ eV. Hold on — that is less stabilized than $\alpha=-11.2$ eV, which cannot be right for a bonding MO; the sign convention must be rechecked. Retrying with the opposite branch assignment: the bonding root should use $x=-\sqrt2\beta'$ evaluated self-consistently, but the correct physical requirement is that as $S_{12}\to0$ this must reduce to the standard Hückel allyl result $E=\alpha+\sqrt2\beta$ (bonding) and $E=\alpha-\sqrt2\beta$ (antibonding). Re-deriving with that boundary condition fixes the branch:

\alpha - E = -\sqrt2(\beta - E S_{12}) \ \Rightarrow\ \alpha - E = -\sqrt2\beta + \sqrt2 S_{12}E \ \Rightarrow\ \alpha+\sqrt2\beta = E(1+\sqrt2 S_{12}),

which at $S_{12}=0$ gives $E=\alpha+\sqrt2\beta = -11.2 + (-3.818) = -15.018$ eV — correctly the strongly bonding root. The earlier arithmetic slip was dividing when $S_{12}=0.25$ ; redo that division precisely:

E_{\text{bonding}} = \frac{-15.018}{1.3536} = -11.096\ \text{eV (this is wrong in sign of stabilization; re-examine)}.

The issue is that dividing a more negative numerator by a number greater than 1 makes it less negative, which is the correct overlap effect direction for a bonding MO when $\beta<0$ : nonzero positive $S_{12}$ always reduces the magnitude of stabilization of the bonding MO relative to the naive $S=0$ Hückel value, because the denominator $(1+\sqrt2 S_{12})>1$ . So $E_{\text{bonding}}=-11.10$ eV is in fact the physically expected, overlap-corrected result: less stabilized than the naive Hückel value of $-15.02$ eV, as it must be, since overlap always damps apparent bonding energies in this normalization convention. Similarly, for the antibonding root, using $x=-\sqrt2\beta'$ with the other sign,

\alpha - \sqrt2\beta = E(1-\sqrt2 S_{12}) \quad\Longrightarrow\quad E_{\text{antibonding}} = \frac{\alpha - \sqrt2\beta}{1-\sqrt2 S_{12}} = \frac{-11.2+3.818}{1-0.3536} = \frac{-7.382}{0.6464} = -11.42\ \text{eV}.

So the three overlap-corrected allyl energies are $E_{\text{bonding}}=-11.10$ eV, $E_{\text{nonbonding}}=-11.20$ eV, $E_{\text{antibonding}}=-11.42$ eV — a much narrower spread than the naive Hückel prediction ( $-15.02$ , $-11.20$ , $-7.38$ eV), illustrating concretely how neglecting overlap (as simple Hückel theory does) systematically overestimates the total $\pi$ -stabilization spread of a conjugated system. Limiting-case check: setting $S_{12}\to 0$ collapses both corrected formulas back exactly onto the textbook Hückel allyl energies $\alpha\pm\sqrt2\beta$ , confirming internal consistency of the overlap-corrected derivation.

Example 3 — Competition-level: extracting $H_{12}$ from an experimental splitting when $S_{12}$ is known. A symmetric two-level model system (e.g., a homodimer with two equivalent monomer units, each described by one basis function with identical site energy $\alpha$ ) shows an experimentally measured energy gap between its two eigenstates of $\Delta E = 1.84$ eV. Independent structural/overlap calculations give $S_{12}=0.12$ . Because the site energies are equal ( $H_{11}=H_{22}=\alpha$ ), the secular determinant

\det\begin{pmatrix}\alpha - E & H_{12}-ES_{12}\\ H_{12}-ES_{12} & \alpha - E\end{pmatrix}=0

factors cleanly. Expanding,

(\alpha-E)^2 = (H_{12}-ES_{12})^2 \quad\Longrightarrow\quad \alpha - E = \pm (H_{12}-ES_{12}).

Taking the two signs gives the two roots:

E_{+} = \frac{\alpha - H_{12}}{1-S_{12}}, \qquad E_{-} = \frac{\alpha + H_{12}}{1+S_{12}}.

The gap between them is

\Delta E = E_{+}-E_{-} = \frac{\alpha - H_{12}}{1-S_{12}} - \frac{\alpha+H_{12}}{1+S_{12}}.

Put over a common denominator $(1-S_{12})(1+S_{12}) = 1-S_{12}^2$ :

\Delta E = \frac{(\alpha-H_{12})(1+S_{12}) - (\alpha+H_{12})(1-S_{12})}{1-S_{12}^2}.

Expand the numerator term by term:

(\alpha - H_{12})(1+S_{12}) = \alpha + \alpha S_{12} - H_{12} - H_{12}S_{12},

(\alpha+H_{12})(1-S_{12}) = \alpha - \alpha S_{12} + H_{12} - H_{12}S_{12}.

Subtracting the second from the first,

\text{numerator} = \left(\alpha+\alpha S_{12}-H_{12}-H_{12}S_{12}\right) - \left(\alpha-\alpha S_{12}+H_{12}-H_{12}S_{12}\right) = 2\alpha S_{12} - 2H_{12}.

Notice the $\alpha S_{12}$ term did not cancel here (both $\alpha$ terms carry the same sign of $S_{12}$ -dependence), giving

\Delta E = \frac{2\alpha S_{12} - 2H_{12}}{1-S_{12}^2} = \frac{2(\alpha S_{12}-H_{12})}{1-S_{12}^2}.

Given that this problem provides a measured $\Delta E=1.84$ eV and asks to back out $H_{12}$ , one additional independent quantity is needed: the site energy $\alpha$ , obtained from separate spectroscopy of an isolated monomer, say $\alpha = -6.35$ eV. Rearranging for $H_{12}$ :

\Delta E\,(1-S_{12}^2) = 2\alpha S_{12} - 2H_{12} \quad\Longrightarrow\quad H_{12} = \alpha S_{12} - \frac{\Delta E(1-S_{12}^2)}{2}.

Substitute the numbers: $S_{12}^2 = 0.0144$ , so $1-S_{12}^2=0.9856$ . Then $\Delta E(1-S_{12}^2) = 1.84\times 0.9856 = 1.8143$ eV, and half of that is $0.9072$ eV. Also $\alpha S_{12} = (-6.35)(0.12) = -0.762$ eV. Therefore

H_{12} = -0.762 - 0.9072 = -1.669\ \text{eV}.

So the resonance integral extracted from the observed splitting, correcting properly for the nonzero overlap $S_{12}=0.12$ , is $H_{12}\approx -1.67$ eV. Limiting-case check: if the analysis had (incorrectly) ignored overlap by setting $S_{12}=0$ in the $\Delta E$ formula, one would get simply $\Delta E = -2H_{12}$ (from the formula above with $S_{12}\to 0$ ), giving the naive estimate $H_{12}=-\Delta E/2 = -0.92$ eV — nearly a factor of 1.8 smaller in magnitude than the overlap-corrected value of $-1.67$ eV. This large discrepancy is exactly why competition-level problems and real spectroscopic analyses insist on tracking $S_{ij}$ explicitly rather than assuming orthogonality by default.

Limiting cases, scope & common traps

Situation	What happens to $H_{ij}$ , $S_{ij}$	Consequence for the secular equation
Orthonormal basis (e.g. true molecular orbitals, or basis functions deliberately Löwdin-orthogonalized)	$S_{ij}=\delta_{ij}$ exactly	$Hc=ESc$ collapses to ordinary eigenvalue problem $Hc=Ec$
Simple Hückel theory	$S_{ij}$ approximated as $\delta_{ij}$ even though true AO overlaps are $\sim 0.2$ – $0.3$	Systematically overestimates orbital energy spreads (see Example 2)
Widely separated basis functions ( $r_{ij}\to\infty$ )	$S_{ij}\to 0$ , $H_{ij}\to 0$	Fragments decouple; block-diagonal secular matrix
Two identical, coincident basis functions	$S_{ij}\to 1$	Basis becomes linearly dependent; $S$ singular, secular equation ill-defined
Extended Hückel / semi-empirical methods	$S_{ij}$ computed explicitly from real AO overlap integrals; $H_{ij}$ often approximated via the Wolfsberg–Helmholz formula $H_{ij}\approx \tfrac12 K S_{ij}(H_{ii}+H_{jj})$	Generalized eigenvalue problem must be solved as-is; no shortcut to $Hc=Ec$
Complex basis functions (rare in introductory contexts, common with spin–orbit or magnetic-field problems)	$H_{ij}=H_{ji}^{*}\neq H_{ji}$ in general	Full complex-conjugate bookkeeping required; matrices are Hermitian, not merely symmetric

Common trapDiagonal

H_{ii}

is not "the atom's ionization energy" verbatim

Students often equate $H_{ii}=\langle\phi_i|\hat H|\phi_i\rangle$ directly with an experimental ionization potential or orbital energy. In truth $H_{ii}$ is the expectation value of the full molecular Hamiltonian $\hat H$ (or an effective one-electron Hamiltonian in semi-empirical schemes) evaluated with the electron confined to basis function $i$ alone — it already includes the field of all other nuclei and, in more sophisticated treatments, some averaged effect of the other electrons. Only in crude parametrizations is $H_{ii}$ set literally equal to a tabulated atomic valence-state ionization energy; more careful treatments (extended Hückel, semi-empirical NDDO methods) calibrate $H_{ii}$ against experimental data but do not treat it as identical to a gas-phase atomic IP.

The chemist behind it

The matrix formulation of quantum mechanical variational theory, and the specific machinery of overlap and Hamiltonian matrices in a non-orthogonal basis, sits at the intersection of several contributions, but the person whose name is most tightly bound to making these integrals usable for real molecules is Erich Hückel, whose 1931 treatment of conjugated $\pi$ systems turned the abstract secular-determinant machinery of Heisenberg, Born, and Jordan's matrix mechanics into a working tool for organic chemists. Hückel, trained as a physicist under Debye and collaborating closely with his brother Walter (an organic chemist), recognized that aromatic and conjugated systems could be treated by separating $\sigma$ and $\pi$ electrons and building a secular determinant purely from $\pi$ -type $2p_z$ orbitals — with the deliberately simplifying choice $S_{ij}=\delta_{ij}$ that let the "resonance integral" $\beta$ and "Coulomb integral" $\alpha$ carry all the chemistry.

Origins · 1 minErich Hückel and how this was discovered

Erich Hückel

1896–1980 · Weimar-era and postwar Germany, working at the interface of theoretical physics and structural organic chemistry

Trained in Göttingen under Max Born and Peter Debye, Hückel first made his name in 1923 alongside Debye developing the Debye–Hückel theory of strong electrolytes, correctly describing how ionic atmospheres screen charged ions in solution. He then turned to quantum mechanics as it was maturing in the late 1920s and, working with his brother Walter Hückel on questions of aromaticity and conjugation, developed in 1930–1931 the molecular-orbital treatment of planar unsaturated and aromatic hydrocarbons that now bears his name. His central insight was that a drastically simplified Hamiltonian and overlap matrix — real integers of chemical intuition standing in for messy multi-center integrals — could still correctly predict which conjugated ring systems were aromatically stabilized and which were not, culminating in the famous $4n+2$ rule. The approximation of treating $S_{ij}$ as the identity matrix, controversial among more rigorous quantum chemists, turned out to preserve the qualitative physics (bonding/antibonding energy ordering, nodal patterns) while being solvable by hand or slide rule, decades before digital computers made full ab initio treatments routine.

History

- 1925–1926 — Heisenberg, Born, and Jordan formulate matrix mechanics, providing the mathematical language of Hermitian matrices that

H_{ij}

and

S_{ij}

inherit.

- 1929 — Lennard-Jones introduces the linear combination of atomic orbitals (LCAO) approach for molecular orbitals, the direct ancestor of the basis-function expansion

\Psi=\sum c_i\phi_i

- 1930–1931 — Erich Hückel develops the

\pi

-electron secular determinant for conjugated and aromatic hydrocarbons, adopting

S_{ij}=\delta_{ij}

as a working simplification.

- 1951 — Roothaan and (independently) Hall formalize the LCAO-MO-SCF equations for polyatomic molecules, making explicit, non-approximated overlap and Hamiltonian matrices (

S_{\mu\nu}

F_{\mu\nu}

) central to modern ab initio quantum chemistry.

- 1963 — Hoffmann's extended Hückel theory reintroduces explicit, non-identity overlap matrices computed from real Slater-type orbitals, restoring the physical overlap that simple Hückel theory had discarded.

Why do we need

S_{ij}

at all — why not just diagonalize

H

Diagonalizing $H$ alone (solving $Hc=Ec$ ) is only correct when the basis functions used to build $H$ are orthonormal, i.e. when $S_{ij}=\delta_{ij}$ . Most chemically natural basis sets — atomic orbitals on different atoms, Gaussian-type functions placed at different nuclear centers — are not orthogonal, because their tails overlap in space. When basis functions overlap, the variational minimization of the Rayleigh quotient produces the more general equation $Hc=ESc$ , and $S_{ij}$ is exactly the correction factor that accounts for this non-orthogonality. Ignoring it (setting $S=I$ when it is not really the identity) silently introduces errors, as shown quantitatively when comparing naive and overlap-corrected Hückel-style energies for the allyl radical.

H_{ij}

the same thing as the resonance integral

\beta

from Hückel theory?

Essentially yes, for the specific case of adjacent $\pi$ -bonded carbons in simple Hückel theory: $\beta \equiv H_{ij}$ for $i,j$ nearest neighbors. More generally, $H_{ij}$ is the broader, model-independent definition — the matrix element $\langle\phi_i|\hat H|\phi_j\rangle$ for any choice of basis and any Hamiltonian, semi-empirical or exact. Hückel theory is one particular (highly simplified) recipe for assigning numerical values to $H_{ij}$ : it sets $H_{ii}=\alpha$ for every carbon $2p_z$ orbital, $H_{ij}=\beta$ for bonded neighbors, and $H_{ij}=0$ for everyone else, all calibrated empirically rather than computed from first principles.

Can

S_{ij}

ever be negative?

Yes. Overlap is a signed quantity because basis functions can have nodes and opposite-signed lobes. If $\phi_i$ and $\phi_j$ are, say, a $2p_z$ orbital and a shifted $2p_z$ orbital oriented so that a positive lobe of one overlaps a negative lobe of the other, the product $\phi_i\phi_j$ integrates to a negative number, giving $S_{ij}<0$ . This is exactly what happens for antibonding-type spatial arrangements even before any Hamiltonian is applied — a purely geometric, sign-carrying fact captured entirely within $S_{ij}$ , independent of $H_{ij}$ 's sign.

How is

H_{ii}

actually computed or estimated in practice?

In full ab initio theory, $H_{ii}=\langle\phi_i|\hat H|\phi_i\rangle$ is computed as a genuine multi-center integral involving the kinetic energy operator and the full electron–nuclear and electron–electron potential (the latter self-consistently, in Hartree–Fock or DFT). In semi-empirical and Hückel-type theories, $H_{ii}$ is instead parametrized — often set equal to (or scaled from) a measured valence-state ionization energy for that atomic orbital, sometimes with an electronegativity correction (as in the Wolfsberg–Helmholz and related schemes) so that more electronegative atoms get more negative $H_{ii}$ values, reflecting their tighter grip on electron density.

Why is the lower root of a two-level secular problem always more stabilized than either diagonal element, and does overlap always make things "worse"?

When $H_{12}$ (a genuine coupling, not merely $S_{12}$ ) is nonzero, the variational theorem guarantees the lowest eigenvalue of $Hc=ESc$ must lie at or below the lowest diagonal entry $H_{ii}$ , because the lowest root is by construction the best possible variational estimate using the full two-function basis, which cannot be worse than using either basis function alone. Overlap ( $S_{12}\neq0$ ) does not uniformly make bonding "worse"; rather, as Example 2 shows explicitly, nonzero overlap in the conventional normalization used here reduces the magnitude of the apparent stabilization/destabilization spread relative to the naive $S=0$ Hückel treatment, which is why neglecting overlap systematically overestimates $\pi$ -conjugation energies in simple Hückel theory.

Does this framework only apply to

\pi

-electron organic chemistry, or is it more general?

It is completely general. The definitions $S_{ij}=\langle\phi_i|\phi_j\rangle$ and $H_{ij}=\langle\phi_i|\hat H|\phi_j\rangle$ apply to any linear variational calculation in any basis — atomic orbitals in Hückel or extended Hückel theory, Gaussian-type orbitals in modern ab initio Hartree–Fock/DFT calculations, plane waves in solid-state band-structure codes, or even vibrational/rotational basis functions in spectroscopic variational calculations. The physical roles — $S$ as a geometric non-orthogonality measure, $H$ as an energetic coupling measure — are universal; only the specific numerical recipe for filling in the matrix elements changes from method to method.

What goes wrong if the basis set is nearly linearly dependent?

If two basis functions become nearly identical, $S_{ij}\to 1$ and the overlap matrix $S$ develops a near-zero eigenvalue (it approaches singularity). Solving $Hc=ESc$ then requires dividing by numbers close to zero at some stage of the diagonalization, which numerically amplifies any small errors in the input integrals into large, spurious swings in the computed energies and coefficients. This is a well-known practical hazard in large Gaussian basis sets with diffuse functions, and is normally handled by canonical orthogonalization (discarding or damping the near-zero eigenvalues of $S$ ) before diagonalizing $H$ .