Hydrogen Molecule

The hydrogen molecule, H₂, is the smallest neutral molecule and the smallest molecular system with electron–electron repulsion. Adding one electron to H₂⁺ changes the problem qualitatively: the electronic equation is no longer separable, the total state must be antisymmetric under electron exchange, and no single orbital product gives both an accurate equilibrium bond and the correct neutral-atom dissociation limit.

That combination makes H₂ the canonical laboratory for covalent bonding and electron correlation. Its ground state is simple enough that localized and delocalized descriptions can be related algebraically, yet rich enough to expose distinctions that remain important in large molecules:

• spin symmetry versus spatial symmetry;
• covalent and ionic structures versus molecular-orbital configurations;
• exchange required by antisymmetry versus correlation beyond a determinant;
• dynamic correlation near equilibrium versus static correlation during dissociation;
• a qualitatively useful orbital picture versus a quantitatively converged wavefunction.

Required background. H₂⁺ provides the two-center Hamiltonian and dissociation convention; molecular orbitals and valence-bond theory supply the two representations. Spin and spatial wavefunctions explains their exchange symmetry.

Two Electrons and Two Nuclei

Place protons $A$ and $B$ at $\pm R\mathbf e_z/2$ and define

$$r_{iA}=|\mathbf r_i-\mathbf R_A|, \qquad r_{iB}=|\mathbf r_i-\mathbf R_B|,$$

for electrons $i=1,2$. In Hartree atomic units, the clamped-nuclei electronic Hamiltonian excluding nuclear repulsion is

$$\hat H_e(R) ={}-\frac12\nabla_1^2-\frac12\nabla_2^2 -\sum_{i=1}^{2} \left(\frac1{r_{iA}}+\frac1{r_{iB}}\right) +\frac1{r_{12}}.$$

The Born–Oppenheimer potential for a normalized electronic state is

$$U(R) = \langle\Psi|\hat H_e(R)|\Psi\rangle +\frac1R.$$

The term $1/r_{12}$ couples the electron coordinates and prevents the separation that works for H₂⁺. Even with clamped nuclei, the exact state is a function of six spatial coordinates and two spin variables,

$$\Psi(x_1,x_2;R), \qquad x_i=(\mathbf r_i,\sigma_i).$$

In the infinite-nuclear-mass convention, two separated ground-state hydrogen atoms have energy

$$U(\infty)=2\left(-\frac12E_h\right)=-1E_h.$$

That threshold is the reference for the Born–Oppenheimer well depth $D_e$. It is not the finite-mass dissociation energy from a specific rovibrational state.

Symmetries of the ground state

The nonrelativistic electrostatic Hamiltonian commutes with total electronic spin, inversion through the midpoint, rotations about the internuclear axis, and electron exchange. The ground electronic term is

$$X\,{}^1\Sigma_g^+.$$

Thus it has total spin $S=0$, zero orbital projection on the molecular axis, even inversion parity, and positive reflection parity. The total electronic wavefunction must also satisfy

$$\hat P_{12}\Psi(x_1,x_2)=-\Psi(x_1,x_2).$$

Nuclear exchange and ortho/para H₂ enter when proton spin and nuclear rotation are restored. They do not alter the electronic antisymmetry condition.

Singlet and Triplet Spatial Functions

The normalized antisymmetric spin singlet is

$$\chi_{00}=\frac{\alpha(1)\beta(2)-\beta(1)\alpha(2)}{\sqrt2}.$$

It multiplies a symmetric spatial function $\Phi_S$. A symmetric triplet spin function $\chi_{1M}$ instead multiplies an antisymmetric spatial function $\Phi_T$, which vanishes at $\mathbf r_1=\mathbf r_2$. The total state is antisymmetric in both cases. A triplet exchange node exists even without Coulomb repulsion; opposite-spin Coulomb avoidance in a singlet requires correlated spatial dependence.

Heitler–London Picture

Let $a(\mathbf r)$ and $b(\mathbf r)$ be normalized real hydrogenic $1s$ orbitals centered on $A$ and $B$, with overlap

$$S=\langle a|b\rangle.$$

Use the shorthand

$$|AB\rangle=a(1)b(2), \qquad |BA\rangle=b(1)a(2).$$

The normalized covalent spatial combinations are

$$\Phi_{\mathrm{cov}}^{+} = \frac{|AB\rangle+|BA\rangle} {\sqrt{2(1+S^2)}},$$ $$\Phi_{\mathrm{cov}}^{-} = \frac{|AB\rangle-|BA\rangle} {\sqrt{2(1-S^2)}}.$$

The plus spatial function pairs with the spin singlet, while the minus spatial function pairs with a triplet spin function:

$$\Psi_{\mathrm{HL}}^{S} = \Phi_{\mathrm{cov}}^{+}\chi_{00},$$ $$\Psi_{\mathrm{HL}}^{T,M} = \Phi_{\mathrm{cov}}^{-}\chi_{1M}.$$

These are the Heitler–London states. Each product term places one electron in an orbital associated with each nucleus, but the two assignments occur coherently. The sum is not a classical mixture of “electron 1 on $A$” and “electron 1 on $B$.”

Direct and exchange matrix elements

Define

$$H_D=\langle AB|\hat H_e|AB\rangle,$$ $$H_X=\langle AB|\hat H_e|BA\rangle.$$

The fixed-nuclei Heitler–London energies are

$$U_S^{\mathrm{HL}} = \frac{H_D+H_X}{1+S^2}+\frac1R,$$ $$U_T^{\mathrm{HL}} = \frac{H_D-H_X}{1-S^2}+\frac1R.$$

Their splitting is

$$U_T^{\mathrm{HL}}-U_S^{\mathrm{HL}} = \frac{2(S^2H_D-H_X)}{1-S^4}.$$

Near the physical bond length, the symmetric spatial state is lower. Its spin partner is therefore a singlet. The splitting vanishes exponentially as the atoms separate because overlap and exchange matrix elements vanish.

Calling this an “exchange force” can be misleading. No additional force has been added to the Coulomb Hamiltonian. The energy difference arises because antisymmetry restricts singlet and triplet states to different spatial subspaces, whose kinetic and Coulomb expectation values differ.

Strength and limitation

Heitler–London theory has the correct neutral dissociation structure:

$$\Psi_{\mathrm{HL}}^{S} \xrightarrow{R\to\infty} \frac{|AB\rangle+|BA\rangle}{\sqrt2}\chi_{00}.$$

It omits ionic charge-fluctuation structures,

$$|AA\rangle=a(1)a(2), \qquad |BB\rangle=b(1)b(2),$$

as well as orbital contraction, polarization, and explicit dependence on $r_{12}$. Ionic configurations are energetically important near equilibrium even though their amplitude must disappear appropriately at neutral dissociation. A pure covalent ansatz therefore gets the asymptote right but is not quantitatively complete in the bonding region.

Molecular-Orbital Picture

Form normalized bonding and antibonding orbitals

$$g = \frac{a+b}{\sqrt{2(1+S)}}, \qquad u = \frac{a-b}{\sqrt{2(1-S)}}.$$

For H₂ these are the minimal $1s\sigma_g$ and $1s\sigma_u$ orbitals. The simplest closed-shell molecular-orbital state doubly occupies $g$:

$$\Psi_{\mathrm{RHF}} = g(1)g(2)\chi_{00}.$$

In a self-consistent calculation this is the restricted Hartree–Fock ansatz: the $\alpha$ and $\beta$ electrons share one optimized spatial orbital. Its minimal-basis spatial part expands as

$$g(1)g(2) =\frac1{2(1+S)}\Big[{} |AA\rangle+|BB\rangle +|AB\rangle+|BA\rangle\Big].$$

The state mixes covalent and ionic structures automatically. Near equilibrium, that flexibility and orbital optimization make the restricted MO description compact and useful. Its simple orbital bond order is

$$\frac{N_g-N_u}{2} = \frac{2-0}{2}=1.$$

The bond-order count describes the configuration; it is not a derivation of the binding energy and does not guarantee correct dissociation.

Exact minimal-basis identities

Define unnormalized gerade singlet structure sums

$$C=|AB\rangle+|BA\rangle,$$ $$I=|AA\rangle+|BB\rangle.$$

Then

$$\begin{aligned} C&=(1+S)|g^2\rangle-(1-S)|u^2\rangle,\\ I&=(1+S)|g^2\rangle+(1-S)|u^2\rangle. \end{aligned}$$

At finite $S$, $C$ and $I$ are nonorthogonal:

$$\frac{\langle C|I\rangle} {\sqrt{\langle C|C\rangle\langle I|I\rangle}} = \frac{2S}{1+S^2}.$$

Consequently, squared coefficients of finite-$R$ covalent and ionic structures are not unique probabilities. Orthogonalization or another population prescription must be declared before numerical “weights” are interpreted.

Heitler–London and MO are basis languages

The two pictures are not competing exact theories. Within a complete one-electron basis and a complete many-electron expansion, localized and delocalized orbitals are related by invertible transformations and represent the same Hilbert space. Their practical difference comes from truncation:

Truncated description	Compact strength	Principal weakness
covalent Heitler–London	neutral dissociation and spin coupling	insufficient ionic, polarization, and dynamic-correlation flexibility
restricted $g^2$ determinant	equilibrium orbital picture and efficient mean field	fixed ionic contamination at dissociation
$g^2/u^2$ configuration interaction	correct static-correlation channel in the minimal space	still limited by the one-electron basis

Exchange and Covalent Bonding

H₂ binds when its total energy, including proton–proton repulsion, falls below the neutral-fragment threshold. Orbital contraction, polarization, ionic fluctuations, and electron–electron avoidance all contribute. The total energy is independent of orbital representation; its division into chemical contributions requires a specified scheme.

Effective spin coupling

When one localized electron is associated with each center, the singlet and triplet multiplets can be represented by

$$\hat H_{\mathrm{spin}} = E_0(R)\mathbf1 +\frac{J(R)}{\hbar^2}\mathbf S_A\cdot\mathbf S_B.$$

For angular momenta $\mathbf S_A,\mathbf S_B$ and an energy-valued $J$,

$$J(R)=E_T(R)-E_S(R).$$

The ordinary H₂ bond has $J(R)>0$ in its bonding region, so the singlet lies lower. This effective Hamiltonian reproduces the singlet–triplet splitting but not the charge fluctuations, orbital relaxation, or excited electronic spectrum from which that splitting arises. As $R\to\infty$, $J(R)\to0$ and the isolated atomic spins become degenerate.

Density alone is not the whole bond

The one-electron density is

$$\rho(\mathbf r) = 2\int |\Phi_S(\mathbf r,\mathbf r_2)|^2 \,d^3r_2.$$

It shows where electronic charge is distributed but integrates out the electron–electron geometry. The pair density

$$P(\mathbf r_1,\mathbf r_2) = 2|\Phi_S(\mathbf r_1,\mathbf r_2)|^2$$

retains that information. Two approximations can have similar $\rho$ and substantially different pair densities, correlation energies, and dissociation behavior. A picture of charge accumulation between nuclei is therefore evidence to interpret, not a complete microscopic definition of a covalent bond.

The Dissociation Problem

The adiabatic continuation of the nondegenerate singlet ground state to large separation produces two neutral hydrogen atoms while retaining total spin and inversion parity. At strict infinite separation the triplet becomes degenerate; a singlet is no longer selected by energy alone. In the $S\to0$ limit, define normalized covalent and ionic structures

$$|C_N\rangle = \frac{|AB\rangle+|BA\rangle}{\sqrt2},$$ $$|I_N\rangle = \frac{|AA\rangle+|BB\rangle}{\sqrt2}.$$

The normalized MO configurations obey

$$|g^2\rangle = \frac{|C_N\rangle+|I_N\rangle}{\sqrt2},$$ $$|u^2\rangle = \frac{-|C_N\rangle+|I_N\rangle}{\sqrt2}.$$

Restricted Hartree–Fock retains only $|g^2\rangle$. It therefore assigns one-half probability to the ionic sector at strict dissociation, even though the physical ground channel is H $+$ H.

Quantitative fixed-orbital failure

For isolated hydrogenic $1s$ orbitals,

$$(1s,1s|1s,1s)=\frac58E_h.$$

In $|g^2\rangle$ at $R\to\infty$, the probability that both electrons occupy the same center is $1/2$. The electron–electron contribution is consequently

$$\langle r_{12}^{-1}\rangle = \frac12\left(\frac58E_h\right) = \frac5{16}E_h.$$

The two one-electron energies sum to $-1E_h$, so the fixed-$1s$ restricted result tends to

$$U_{\mathrm{RHF}}(\infty) = -1+\frac5{16} = -\frac{11}{16}E_h.$$

The error is not a small basis-set imperfection: the ansatz has the wrong charge-fluctuation structure. Orbital optimization changes the numerical plateau but cannot remove the restricted determinant’s qualitative dissociation defect.

Three repairs with different tradeoffs

Method	Large-$R$ behavior
Heitler–London or spin-projected localized pair	preserves singlet symmetry and yields neutral fragments
unrestricted Hartree–Fock	reaches the neutral-fragment energy by localizing opposite spins, but breaks spin symmetry
two-configuration CI or multiconfigurational SCF	mixes $g^2$ and $u^2$ and preserves both neutral dissociation and total spin

In unrestricted Hartree–Fock, a determinant such as $|a\alpha,b\beta|$ becomes energetically favorable beyond a basis-dependent Coulson–Fischer point. At infinite separation it has the correct energy and charge density, but it is not an eigenstate of $\hat S^2$:

$$\frac{\langle\hat S^2\rangle}{\hbar^2} \xrightarrow{R\to\infty}1.$$

For two electrons this corresponds to equal singlet and triplet contamination. Broken symmetry can be a useful mean-field device, but its determinant is not the exact spin state.

The exact dissociation state is entangled

The spin-adapted neutral limit can be written schematically as

$$|\Psi_{S=0}\rangle = \frac{ |A\alpha;B\beta\rangle -|A\beta;B\alpha\rangle }{\sqrt2}.$$

No single spin-orbital determinant with one localized electron on each center represents this pure singlet. In the symmetry-adapted MO basis, the same fact appears through equal-magnitude $g^2$ and $u^2$ amplitudes. The spatial natural-orbital occupations approach

$$n_g\to1, \qquad n_u\to1,$$

rather than the restricted-determinant values $n_g=2$, $n_u=0$. The multiconfigurational character is therefore basis-invariant at the level of the one-body density-matrix spectrum, even though particular “covalent” and “ionic” weights are not.

Correlation and Configuration Interaction

At a fixed geometry, the conventional correlation energy is

$$E_c(R) ={}E_{\mathrm{exact}}^{\mathrm{BO}}(R) -E_{\mathrm{HF}}(R),$$

where both energies use the same nonrelativistic clamped-nuclei Hamiltonian and each is converged in its respective complete one-electron basis. Specify whether the HF reference is restricted to a singlet determinant or permits unrestricted spin breaking; the two reference energies differ beyond the Coulson–Fischer point. Because Hartree–Fock is variational,

$$E_c(R)\leq0.$$

For H₂, two regimes are especially clear.

Dynamic correlation near equilibrium

Near $R_e$, the $g^2$ determinant is dominant, but it treats opposite-spin repulsion only through an average field. The exact wavefunction suppresses configurations in which the electrons approach one another and enhances configurations in which they avoid each other. Away from a nucleus, spherical averaging over the electron-pair relative direction at fixed pair center gives $\overline\Phi_S(r_{12})$. It satisfies the electron–electron cusp condition

$$\left. \frac{\partial\overline\Phi_S}{\partial r_{12}} \right|_{r_{12}=0} = \frac12 \Phi_S(r_{12}=0).$$

Products of smooth one-electron orbitals do not reproduce this linear cusp efficiently. Orbital configuration interaction converges toward it through many angular momenta, while explicitly correlated $R12/F12$, Hylleraas, James–Coolidge, or Jastrow factors include $r_{12}$ directly.

Static correlation during dissociation

As $R$ increases, the $g$ and $u$ orbital energies approach degeneracy. The $g^2$ and $u^2$ configurations must then enter with comparable amplitudes. This is static, nondynamic, or strong correlation: more than one configuration is required before short-range dynamical avoidance is considered.

The distinction is diagnostic rather than absolute. Near equilibrium, many small excitations improve the Coulomb hole; near dissociation, one particular double excitation becomes order one. A method designed only for small corrections to one determinant can therefore work near $R_e$ and fail catastrophically as the bond is stretched.

Minimal-basis full CI

In the gerade singlet sector of the two-orbital minimal basis, use the orthonormal configurations

$$|G\rangle=|g^2\rangle, \qquad |U\rangle=|u^2\rangle.$$

The full-CI Hamiltonian in this symmetry sector is

$$\mathbf H_{\mathrm{CI}} = \begin{pmatrix} H_{GG}&H_{GU}\\ H_{GU}&H_{UU} \end{pmatrix}.$$

Define the mean diagonal energy and half-difference by

$$\begin{aligned} \bar H&=\frac{H_{GG}+H_{UU}}2,\\ \delta_H&=\frac{H_{GG}-H_{UU}}2. \end{aligned}$$

The lower eigenvalue is then

$$E_-=\bar H-\sqrt{\delta_H^2+H_{GU}^2}.$$

The eigenvector

$$|\Psi_-\rangle = c_G|G\rangle+c_U|U\rangle$$

is exact only within the chosen two-orbital one-particle space. Near equilibrium $|c_G|$ dominates. At dissociation,

$$c_G\to\frac1{\sqrt2}, \qquad c_U\to-\frac1{\sqrt2},$$

which cancels the ionic structure and leaves $|C_N\rangle$.

Full CI means all symmetry-allowed determinants or configuration-state functions in a specified finite orbital basis. It removes determinant truncation error in that basis but not basis-set error. For two electrons, the space is small enough to make that distinction unusually transparent.

Explicitly correlated H₂ wavefunctions

The high-accuracy route does not stop at a two-orbital CI. James–Coolidge and related bases use the two-center coordinates

$$\xi_i=\frac{r_{iA}+r_{iB}}R, \qquad \eta_i=\frac{r_{iA}-r_{iB}}R,$$

and include $r_{12}$ explicitly. A schematic basis function is

$$\Phi_k = e^{-\alpha(\xi_1+\xi_2)} \xi_1^{p_k}\xi_2^{q_k} \eta_1^{m_k}\eta_2^{n_k} r_{12}^{s_k},$$

followed by the required exchange and inversion symmetrization. Increasing the polynomial space gives variationally convergent Born–Oppenheimer energies. Explicit $r_{12}$ powers resolve the Coulomb cusp far more efficiently than a determinant expansion built only from smooth orbitals.

Quantitative Benchmarks

Accurate nonrelativistic clamped-nuclei calculations locate the minimum near

$$R_e\simeq1.4011a_0 \simeq0.7414\,\text{\AA},$$ $$U^{\mathrm{BO}}(1.4011a_0) \simeq -1.1744759314002167E_h.$$

The displayed high-precision value is at $R=1.4011a_0$, not a claim that this rounded geometry is the exact minimizer. Relative to two infinite-mass hydrogen atoms, the well depth to the shown precision is

$$\begin{aligned} D_e^{\mathrm{BO}} &=(-1E_h)-U_e^{\mathrm{BO}}\\ &=0.1744759314E_h\\ &\simeq4.74773\,\mathrm{eV}. \end{aligned}$$

At $R=1.4a_0$, representative complete-basis values are

$$E_{\mathrm{HF}}\simeq-1.13362957E_h,$$ $$E_{\mathrm{exact}}^{\mathrm{BO}} = -1.174475714220E_h,$$

so

$$E_c\simeq-0.04084614E_h \simeq-1.1115\,\mathrm{eV}.$$

This correlation energy is evaluated at one geometry. It should not be silently combined with an independently optimized Hartree–Fock minimum or a finite-mass experimental threshold.

For comparison, Liu and colleagues (2009) obtained the following hybrid experimental–theoretical determination of the dissociation energy from $v=0,J=0$ in ordinary H₂:

$$D_0=36118.06962(37)\,\mathrm{cm}^{-1} \simeq4.47807\,\mathrm{eV}.$$

The difference from $D_e^{\mathrm{BO}}$ includes nuclear zero-point motion and finite-mass, adiabatic, nonadiabatic, relativistic, and radiative effects under the chosen convention. Ortho-H₂ begins at $J=1$ and therefore has a different quoted $D_0$; isotope and rotational state are part of the datum.

Exercises

1. Derive the singlet–triplet splitting

Use the direct and exchange matrix elements to derive $U_T^{\mathrm{HL}}-U_S^{\mathrm{HL}}$. Why does nuclear repulsion not appear in the final splitting?

Solution

Subtracting the two variational energies gives

$$\begin{aligned} U_T-U_S ={}& \frac{H_D-H_X}{1-S^2} -\frac{H_D+H_X}{1+S^2}\\ ={}& \frac{2(S^2H_D-H_X)}{1-S^4}. \end{aligned}$$

The same scalar $1/R$ is added to both electronic symmetry sectors, so it cancels in their difference.

2. Optimize the restricted dissociation plateau

Use normalized $1s$ functions with a common exponent $\zeta$ on the far-separated centers. The one-electron energy of each function is $\zeta^2/2-\zeta$ and its same-center Coulomb integral is $5\zeta/8$. Find the optimum restricted energy. Does this orbital relaxation recover neutral dissociation?

Solution

The restricted state retains ionic probability $1/2$, so

$$E(\zeta)=\zeta^2-2\zeta+\frac{5\zeta}{16} =\zeta^2-\frac{27}{16}\zeta.$$

Its minimum occurs at $\zeta=27/32$, with $E=-729/1024\,E_h\simeq-0.711914\,E_h$. This improves the fixed-exponent $-11/16\,E_h$ result, but remains far above $-1\,E_h$. Relaxing this orbital family cannot eliminate the determinant’s ionic sector.

3. Natural occupations as a correlation diagnostic

Compare the limiting spatial natural-orbital occupations of restricted $g^2$ with those of the exact neutral singlet at dissociation. Why can no orbital rotation turn the latter into one determinant?

Solution

For a doubly occupied determinant,

$$n_g=2, \qquad n_u=0.$$

For the dissociated pure singlet,

$$n_g=n_u=1.$$

The spatial spectrum rules out a closed-shell determinant. To exclude every spin-orbital determinant, use the pure singlet’s spin-orbital occupations: $1/2,1/2,1/2,1/2$, instead of the idempotent spectrum $1,1,0,0$ required by any two-electron determinant. Those eigenvalues are unchanged by any spin-orbital basis rotation. Spatial occupations $1,1$ alone would not exclude an unrestricted determinant.

4. Correlation and dissociation conventions

At $R=1.4a_0$, use the quoted Hartree–Fock and accurate Born–Oppenheimer energies to calculate $E_c$. Then explain why subtracting the spectroscopic $D_0$ from the Born–Oppenheimer $D_e$ is not merely a measure of electronic correlation.

Solution

At the same geometry,

$$\begin{aligned} E_c &=-1.174475714220E_h -(-1.13362957E_h) \\ &\simeq-0.04084614E_h. \end{aligned}$$

This is about $-1.1115\,\mathrm{eV}$. By contrast, $D_e^{\mathrm{BO}}-D_0$ includes nuclear zero-point energy, finite nuclear masses, nonadiabatic effects, and precision relativistic and radiative corrections. It compares different physical models and energy levels, not Hartree–Fock with an exact electronic solution at one geometry.

References

• P. W. Atkins and R. S. Friedman, Molecular Quantum Mechanics, 5th ed. (Oxford University Press, 2011), Chapters 8–9.
• C. A. Coulson and I. Fischer, “Notes on the Molecular Orbital Treatment of the Hydrogen Molecule,” Philosophical Magazine 40, 386–393 (1949), doi:10.1080/14786444908521726.
• W. Heitler and F. London, “Wechselwirkung neutraler Atome und homöopolare Bindung nach der Quantenmechanik,” Zeitschrift für Physik 44, 455–472 (1927), doi:10.1007/BF01397394.
• H. M. James and A. S. Coolidge, “The Ground State of the Hydrogen Molecule,” Journal of Chemical Physics 1, 825–835 (1933), doi:10.1063/1.1749252.
• T. Kato, “On the Eigenfunctions of Many-Particle Systems in Quantum Mechanics,” Communications on Pure and Applied Mathematics 10, 151–177 (1957), doi:10.1002/cpa.3160100201.
• W. Kołos and C. C. J. Roothaan, “Accurate Electronic Wave Functions for the H₂ Molecule,” Reviews of Modern Physics 32, 219–232 (1960), doi:10.1103/RevModPhys.32.219.
• J. Liu, E. J. Salumbides, U. Hollenstein, J. C. J. Koelemeij, K. S. E. Eikema, W. Ubachs, and F. Merkt, “Determination of the Ionization and Dissociation Energies of the Hydrogen Molecule,” Journal of Chemical Physics 130, 174306 (2009), doi:10.1063/1.3120443.
• P.-O. Löwdin, “Quantum Theory of Many-Particle Systems. III. Extension of the Hartree–Fock Scheme to Include Degenerate Systems and Correlation Effects,” Physical Review 97, 1509–1520 (1955), doi:10.1103/PhysRev.97.1509.
• R. McWeeny, Methods of Molecular Quantum Mechanics, 2nd ed. (Academic Press, 1992), Chapters 4–7.
• R. S. Mulliken, “The Assignment of Quantum Numbers for Electrons in Molecules. I,” Physical Review 32, 186–222 (1928), doi:10.1103/PhysRev.32.186.
• K. Pachucki, “Born–Oppenheimer Potential for H₂,” Physical Review A 82, 032509 (2010), doi:10.1103/PhysRevA.82.032509.
• M. Puchalski, J. Komasa, and K. Pachucki, “Nonadiabatic QED Correction to the Dissociation Energy of the Hydrogen Molecule,” Physical Review Letters 122, 103003 (2019), doi:10.1103/PhysRevLett.122.103003.
• J. S. Sims and S. A. Hagstrom, “High Precision Variational Calculations for the Born–Oppenheimer Energies of the Ground State of the Hydrogen Molecule,” Journal of Chemical Physics 124, 094101 (2006), doi:10.1063/1.2173250.
• A. Szabo and N. S. Ostlund, Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory (Dover, 1996), Chapters 2–4.