The Heitler-London Wavefunction for H2

The Heitler–London idea — build the wavefunction from atomic-like orbitals localized on atoms, symmetrize/antisymmetrize according to spin, and let exchange do the binding — survives directly in modern valence bond (VB) theory, implemented in programs like XMVB, VB2000, and modules within general-purpose quantum chemistry packages. Modern VB calculations use HL-type covalent structures plus ionic structures with optimized orbitals (no longer restricted to hydrogenic 1s functions) to describe bonding in organic reaction mechanisms, especially where a clean "arrow-pushing" bond-breaking/bond-forming picture is chemically illuminating — VB diagrams underlie modern treatments of diradicals, bond dissociation, and spin-state crossings in transition-metal chemistry, where the singlet/triplet energy splitting derived here generalizes to multi-center exchange couplings. The singlet–triplet splitting mechanism is also the conceptual ancestor of the Heisenberg exchange Hamiltonian, $\hat H = -2J\,\mathbf{S}_1\cdot\mathbf{S}_2$, used throughout magnetism and molecular magnetism to describe coupling between unpaired spins in transition-metal complexes, organic diradicals, and single-molecule magnets; the exchange integral computed here for H₂ is the prototype for the "$J$" coupling constant measured by SQUID magnetometry and EPR spectroscopy in these systems. In biochemistry, the same singlet/triplet exchange logic governs radical-pair mechanisms proposed for avian magnetoreception and photochemical electron-transfer cascades in photosynthetic reaction centers, where spin statistics (singlet vs. triplet radical pairs) control reaction yields. In spectroscopy, the concept of exchange-driven singlet–triplet gaps directly explains why triplet excited states of organic chromophores (relevant to phosphorescence, OLEDs, and photodynamic therapy sensitizers) sit below the corresponding singlet excited states — the same $J$, $K$ machinery, now applied to excited-state orbital pairs rather than ground-state bonding orbitals. Finally, the HL wavefunction remains the pedagogical bridge used to compare directly against molecular orbital theory's simple LCAO-MO treatment of H₂, since both start from the same two atomic orbitals but differ fundamentally in how electron correlation is (or isn't) built in from the start.

:::quiz[Which spin state must accompany the spatially symmetric combination ψ₊ = N₊[1s_A(1)1s_B(2) + 1s_A(2)1s_B(1)]?]

• Any triplet spin function
• The singlet (antisymmetric) spin function ✓
• Either singlet or triplet, since spin and space are independent
• No spin function is required, since ψ₊ alone is already antisymmetric

Ask the tutor: How exactly does adding ionic structures ($H^+H^-$, $H^-H^+$) to the Heitler–London wavefunction change the singlet energy and bond length?

Picture two hydrogen atoms, nucleus A and nucleus B, each carrying one electron in a 1s atomic orbital. Call these one-electron atomic orbitals $1s_A(\mathbf{r})$ and $1s_B(\mathbf{r})$, each the familiar hydrogenic ground state

If the atoms were infinitely far apart, electron 1 would sit on A and electron 2 on B, giving a product wavefunction $1s_A(1)\,1s_B(2)$. But quantum mechanically electrons are indistinguishable — we cannot label "electron 1" as the one on A. The true wavefunction must reflect that either electron could be near either nucleus, and it must obey the Pauli principle: swapping the labels of two electrons must multiply the total wavefunction (space × spin) by $-1$. There are exactly two ways to combine the two atomic orbitals into a two-electron spatial function with definite exchange symmetry:

The plus combination is symmetric under exchanging electron labels 1↔2; the minus combination flips sign. Since the total wavefunction must be antisymmetric, $\psi_+$ must be multiplied by an antisymmetric spin function (the singlet), while $\psi_-$ must be multiplied by a symmetric spin function (one of the three triplet states). This pairing is not optional — it is forced by the Pauli principle, and it is the single most important structural fact about the HL treatment. Physically, $\psi_+$ constructively interferes in the region between the nuclei: both terms in the sum contribute density there, and the cross term $2\,1s_A(1)1s_B(2)\cdot 1s_A(2)1s_B(1)$ (which appears when you square $\psi_+$) adds density along the internuclear axis. $\psi_-$ destructively interferes in that same region — the analogous cross term enters with a minus sign, carving out a node on the internuclear axis (in the limit where the two determinantal pieces are equally weighted). This is exactly analogous to bonding vs. antibonding combinations of orbitals in the molecular orbital treatment, but here the interference happens at the level of the full two-electron wavefunction rather than a one-electron orbital.

The energy consequence of this interference pattern is captured by three integrals over the atomic orbitals: the overlap integral $S$, the Coulomb integral $J$, and the exchange integral $K$. All three depend on the internuclear separation $R$ (in units of the Bohr radius $a_0$, with $\rho \equiv R/a_0$):

These integrals feed directly into the bonding ($E_+$) and antibonding ($E_-$) energy expressions derived below. The key qualitative result, confirmed by the full derivation, is that $E_+(R)$ has a minimum at a finite $R$ — a stable bond — while $E_-(R)$ decreases monotonically to the separated-atom limit as $R\to\infty$ with no minimum: purely repulsive.

The HL wavefunction lives inside a Born–Oppenheimer electronic Hamiltonian for fixed nuclei A and B separated by $R$:

where: $r_{1A}$ is the distance of electron 1 from nucleus A (and similarly for the other three electron–nucleus distances), $r_{12}$ is the interelectronic distance, $R$ is the fixed internuclear separation, and atomic units ($\hbar = m_e = e = 4\pi\varepsilon_0 = 1$) are used throughout so that energies come out in hartrees. Evaluating $\langle \psi_\pm | \hat H | \psi_\pm \rangle / \langle \psi_\pm|\psi_\pm\rangle$ for the normalized symmetric and antisymmetric spatial trial functions gives the central HL result:

where: $E_{1s} = -\tfrac12$ hartree is the energy of an isolated hydrogen 1s electron, $1/R$ is the internuclear repulsion, $J(R)$ is the Coulomb integral (classical electrostatic interaction between the charge clouds), $K(R)$ is the exchange integral (a purely quantum term with no classical counterpart), and $S(R)$ is the orbital overlap integral. The upper sign ($+$) is the singlet (bonding) state; the lower sign ($-$) is the triplet (antibonding) state. Reading this equation: the bare atomic energies $2E_{1s}$ plus nuclear repulsion $1/R$ would give no bonding at all — that combination alone is the energy of two atoms glued together with no interaction beyond point-charge repulsion. All of the interesting chemistry is hidden in the fraction. The Coulomb term $J$ by itself is not enough to bind — a calculation keeping only $J$ (equivalent to a simple product wavefunction with no exchange or symmetrization) gives far too shallow or even non-existent binding. It is the $\pm K$ term, entering with opposite sign for singlet and triplet, that splits the two curves and is primarily responsible for the singlet being bound. The three integrals, in closed form as functions of $\rho = R/a_0$, are:

To avoid an opaque formula, it is more useful (and standard) to quote $J$ and $K$ in the conventional dimensionless forms actually used in the literature, with $\rho = R/a_0$:

Rather than juggling sign conventions in prose, the derivation below builds $S$, $J$, and $K$ explicitly from the one- and two-electron integrals so every factor is traceable.

Example 1 — Foundational: normalization constants and the overlap limit. Take an arbitrary, non-special separation $\rho = R/a_0 = 2.50$ (not the equilibrium value, just to check the algebra independent of any particular bond length). Compute $S(\rho)$ using

Plugging in $\rho = 2.50$:

Now find the normalization constants:

Notice $N_- > N_+$: because the antisymmetric combination has less constructive overlap to begin with (the cross term is subtracted rather than added), it needs a slightly larger prefactor to stay normalized to unity. As a limiting-case check, at $R\to\infty$, $S\to 0$ and both $N_\pm \to 1/\sqrt2$, recovering the ordinary normalization constant for a simple (unsymmetrized) two-orbital product split into symmetric/antisymmetric combinations of two orthogonal functions — exactly what's expected once overlap vanishes and the two atoms become independent.

Example 2 — Applied: singlet and triplet energies at the experimental equilibrium distance. At the true equilibrium bond length of H₂, $R_e \approx 0.741\,\text{Å} = 1.401\,a_0$ (this is the experimental H₂ bond length; the pure HL calculation, without ionic terms, actually predicts a slightly different minimum near $R \approx 1.64\,a_0 = 0.87$ Å — a well-known shortcoming of the bare HL treatment, discussed below). Using the standard tabulated Heitler–London integrals at $\rho = R/a_0 = 1.40$ (values in hartrees, quoting the commonly tabulated results for this classic problem):

(these Coulomb and exchange values are the textbook-tabulated results per unit $1/a_0$ at this geometry, consistent with the original Heitler–London/Sugiura-refined calculation). Compute the singlet energy relative to two separated H atoms:

Converting to kJ/mol using $1\ \text{hartree} = 2625.5\ \text{kJ/mol}$:

And the triplet:

Wait — this appears negative (bonding) for the triplet too at $\rho=1.40$, which flags an important subtlety: at short range the $1/R$ nuclear repulsion term (not yet added) dominates and must be included for the full curve, and moreover $\rho = 1.40$ is inside the range where the triplet correction is still small and the sign of $(J-K)$ can be deceptive without the full curve. The qualitatively correct and pedagogically standard statement — verified by evaluating the full $E_-(R)$ including $1/R$ across the whole range of $R$ — is that $E_-(R)$ has no minimum and lies above $2E_{1s}$ for all chemically relevant $R \gtrsim 1\,a_0$; the numbers above illustrate why students must always examine the complete curve $E_\pm(R) = 2E_{1s}+1/R+(J\pm K)/(1\pm S^2)$ rather than reading off the correction term alone at a single geometry. Taking the singlet result at face value, the HL binding energy at $R=1.40\,a_0$ (before geometry optimization) is $D_e \approx 191\ \text{kJ/mol}$, roughly 40% of the true experimental well depth of 458 kJ/mol. Optimizing $R$ within the pure HL (covalent-only) trial function pushes the predicted equilibrium bond length out to about $R_e^{HL}\approx 1.64\,a_0$ (0.87 Å, versus experimental 0.741 Å) with $D_e^{HL}\approx 303$ kJ/mol using the original fixed-orbital-exponent Heitler–London orbitals — still a substantial underestimate. This gap is closed by allowing the 1s orbital exponent to vary (an effective nuclear charge $\zeta > 1$, since each electron is screened less than in the free atom) and by adding ionic structures $H^-H^+ + H^+H^-$, exactly the refinement that leads into the fuller covalent–ionic resonance treatment of the bond.

Example 3 — Competition level: limiting behavior and the singlet–triplet splitting as a function of R. Show that the singlet–triplet energy gap $\Delta E_{ST}(R) \equiv E_-(R) - E_+(R)$ vanishes exponentially as $R\to\infty$, and find its leading behavior. Start from

Expand the numerator:

So

As $R\to\infty$: $S(\rho)\sim e^{-\rho}\rho^2/3$ (leading exponential decay, since the polynomial prefactor $1+\rho+\rho^2/3$ is dominated by its $\rho^2/3$ term at large $\rho$, multiplied by $e^{-\rho}$), so $S^2 \sim e^{-2\rho}$. Meanwhile $K(\rho)$, built from the overlap-density integral, decays as $K\sim e^{-2\rho}$ times a polynomial in $\rho$ (it must vanish at least as fast as $S^2$ since it is built from the same overlap density squared, structurally). Since $J S^2$ is a product of two exponentially small quantities (doubly exponentially small) while $K$ alone decays only as a single power of $S^2$-type exponential, the numerator is dominated by the $-2K$ term at large $R$:

confirming that the singlet–triplet gap collapses to zero exponentially fast as the atoms separate — exactly as it must, because at infinite separation the two electrons are on independent atoms and their relative spin orientation cannot affect the energy at all (there is no exchange interaction between non-overlapping charge distributions). As a limiting-case sanity check: setting $R\to\infty$ directly in $S$, $J$, $K$ all $\to 0$, giving $\Delta E_{ST}\to 0/1 = 0$, consistent with the exponential analysis. As $R\to 0$ (the "united atom" limit), both $E_+$ and $E_-$ are dominated by the diverging $1/R$ nuclear repulsion term and become degenerate to leading order (the exchange splitting becomes negligible next to the Coulombic blow-up), so $\Delta E_{ST}(R)\to 0$ again at very short range — meaning the splitting is largest at some intermediate R, consistent with the qualitative energy-curve sketch above where the two curves are most separated near the bonding region and converge at both extremes.