Variational principle in quantum chemistry: the secular determinant and basis set expansion

Intuition Beginner

Most quantum chemistry problems cannot be solved exactly. The variational principle provides a systematic way to approximate the answer. It states a simple but powerful fact: any trial wave function you write down will give an energy that is greater than or equal to the true ground-state energy. You can never undershoot. This means the best wave function is the one that gives the lowest energy -- you minimise to approach the truth from above.

To use the variational principle, you need a flexible family of trial wave functions. The linear variational method builds the trial function as a linear combination of fixed basis functions: $\tilde{\psi }=c_1{\varphi }_1+c_2{\varphi }_2+\cdots +c_N{\varphi }_N$. The coefficients $c_i$ are the adjustable parameters. Minimising the energy with respect to these coefficients leads to a matrix equation called the secular equation, whose solutions are the variational energies and wave functions.

Visual Beginner

The variational principle has a geometric interpretation. The exact eigenstates of the Hamiltonian form an orthonormal basis for the Hilbert space. A trial wave function is some vector in that space. The expectation value $⟨\tilde{\psi }|\hat{H}|\tilde{\psi }⟩$ is a weighted average of all eigenvalues, with weights given by the overlap of $\tilde{\psi }$ with each eigenstate. Because the weights are positive and sum to one, this weighted average can never fall below the lowest eigenvalue.

For the linear variational method with two basis functions, the secular determinant is a $2\times 2$ matrix. Setting its determinant to zero gives a quadratic equation whose two roots are approximations to the two lowest energies. Adding more basis functions enlarges the matrix and improves the approximation.

Worked example Beginner

Apply the linear variational method to the H${}_2^{+}$ molecular ion using a two-function basis: $\tilde{\psi }=c_1{\varphi }_{1s}^A+c_2{\varphi }_{1s}^B$, where ${\varphi }_{1s}^A$ and ${\varphi }_{1s}^B$ are hydrogen 1s orbitals centred on nuclei A and B.

Step 1. Compute the matrix elements. $H_{11}=⟨{\varphi }_{1s}^A|\hat{H}|{\varphi }_{1s}^A⟩$ is the energy of an electron in orbital A (the Coulomb integral $\alpha$). $H_{12}=⟨{\varphi }_{1s}^A|\hat{H}|{\varphi }_{1s}^B⟩$ is the resonance integral $\beta$. $S_{12}=⟨{\varphi }_{1s}^A|{\varphi }_{1s}^B⟩$ is the overlap integral $S$.

Step 2. The secular determinant is

$$\left|\begin{array}{cc}\alpha -E& \beta -ES\\ \beta -ES& \alpha -E\end{array}\right|=0.$$

Step 3. Expanding: $(\alpha -E{)}^2-(\beta -ES{)}^2=0$, giving $E_{\pm }=\frac{\alpha \pm \beta }{1\pm S}$. Since $\beta <0$ and $S>0$, the lower root $E_{-}=(\alpha +\beta )/(1+S)$ is the bonding orbital energy and the upper root $E_{+}=(\alpha -\beta )/(1-S)$ is the antibonding energy.

Step 4. The corresponding wave functions are ${\psi }_{-}=({\varphi }_{1s}^A+{\varphi }_{1s}^B)/\sqrt{2(1+S)}$ (bonding, electron density between the nuclei) and ${\psi }_{+}=({\varphi }_{1s}^A-{\varphi }_{1s}^B)/\sqrt{2(1-S)}$ (antibonding, a node between the nuclei).

Check your understanding Beginner

Formal definition Intermediate+

Variational principle (general). Let $\hat{H}$ be a time-independent Hamiltonian with discrete, ordered eigenvalues $E_0\le E_1\le E_2\le \cdots$ and corresponding orthonormal eigenstates $\{|k⟩\}$. For any normalised trial state $|\tilde{\psi }⟩$,

$$⟨\tilde{\psi }|\hat{H}|\tilde{\psi }⟩\ge E_0.$$

Proof. Expand $|\tilde{\psi }⟩={\sum }_k{c}_k|k⟩$ with ${\sum }_k|c_k{|}^2=1$. Then

$$⟨\tilde{\psi }|\hat{H}|\tilde{\psi }⟩=\sum _k|c_k{|}^2{E}_k\ge E_0\sum _k|c_k{|}^2=E_0.$$

Equality holds if and only if $|\tilde{\psi }⟩=|0⟩$ (up to a phase). $\square$

The linear variational method

Choose a finite set of $N$ linearly independent basis functions $\{{\varphi }_i{\}}_{i=1}^N$ and form the trial function $\tilde{\psi }={\sum }_{i=1}^N{c}_i{\varphi }_i$. The energy is

$$E(\mathbf{c})=\frac{{\sum }_{i,j}{c}_i^{\ast }{c}_j{H}_{ij}}{{\sum }_{i,j}{c}_i^{\ast }{c}_j{S}_{ij}},$$

where $H_{ij}=⟨{\varphi }_i|\hat{H}|{\varphi }_j⟩$ is the Hamiltonian matrix element and $S_{ij}=⟨{\varphi }_i|{\varphi }_j⟩$ is the overlap matrix element. Minimising $E$ with respect to each $c_k$ (setting $\partial E/\partial c_k^{\ast }=0$) yields the secular equations

$$\sum _{j=1}^N(H_{ij}-E{S}_{ij})c_j=0,i=1,\dots ,N.$$

This system has a nontrivial solution if and only if the secular determinant vanishes:

$$\det (\mathbf{H}-E\mathbf{S})=0.$$

This is a generalised eigenvalue problem. The $N$ roots $E_0\le E_1\le \cdots \le E_{N-1}$ are variational upper bounds to the $N$ lowest exact eigenvalues. The corresponding eigenvectors ${\mathbf{c}}^{(k)}$ give the variational wave functions.

Basis functions for molecular calculations

The basis functions used in quantum chemistry fall into two families:

Slater-type orbitals (STOs) have the form

$${\chi }_{nlm}^{\text{STO}}(r,\theta ,\phi )=N{r}^{n-1}{e}^{-\zeta r}{Y}_l^m(\theta ,\phi ),$$

where $\zeta$ is the orbital exponent. STOs have the correct cusp at the nucleus ($d\psi /dr\ne 0$ at $r=0$ for $s$ orbitals) and the correct exponential decay at large $r$. However, the two-electron repulsion integrals $⟨{\chi }_a{\chi }_b|r_{12}^{-1}|{\chi }_c{\chi }_d⟩$ over STOs are expensive to evaluate for polyatomic molecules because there is no general analytical formula for multi-centre integrals.

Gaussian-type orbitals (GTOs) replace the exponential with a Gaussian:

$${\chi }_{l_x,l_y,l_z}^{\text{GTO}}(x,y,z)=N{x}^{l_x}{y}^{l_y}{z}^{l_z}{e}^{-\alpha r^2}.$$

Gaussians lack the nuclear cusp and decay too quickly at large $r$, but the product of two Gaussians centred at different points is another Gaussian centred at an intermediate point. This Gaussian product theorem reduces all multi-centre integrals to single-centre integrals, making GTOs computationally efficient. The price is that more Gaussian functions are needed to approach the quality of a single STO.

Contracted Gaussian basis sets

A contracted Gaussian function is a fixed linear combination of primitive Gaussians:

$${\chi }_{\mu }^{\text{CGF}}=\sum _{p=1}^K{d}_{\mu p}{\chi }_p^{\text{PGF}}({\alpha }_p),$$

where the contraction coefficients $d_{\mu p}$ and exponents ${\alpha }_p$ are determined once (usually by optimising atomic energies) and then held fixed during molecular calculations. The hierarchy of Gaussian basis sets includes:

Basis type	Description	Example for C
Minimal (STO-$n$G)	One contracted function per atomic orbital	STO-3G: 5 functions
Split-valence	Two (or more) functions per valence AO, one per core AO	3-21G, 6-31G
Triple-zeta	Three functions per valence AO	6-311G
Polarisation	Adds functions with $l$ one higher than valence	6-31G(d), 6-31G(d,p)
Diffuse	Adds small-exponent functions for Rydberg/anion states	6-31+G(d), aug-cc-pVDZ

The notation 6-31G(d) means: core orbitals use 6 primitives contracted to 1 function, valence orbitals are split into two contractions (3 primitives + 1 primitive), and $d$-type polarisation functions are added to heavy atoms. The "(d,p)" variant also adds $p$-type functions to hydrogen.

Key theorem with proof Intermediate+

Theorem (MacDonald's theorem). The $k$-th root $E_k$ of the secular equation $\det (\mathbf{H}-E\mathbf{S})=0$ is an upper bound to the $(k+1)$-th exact eigenvalue $E_k^{\text{exact}}$ of $\hat{H}$:

$$E_k\ge E_k^{\text{exact}},k=0,1,\dots ,N-1.$$

Proof. Let $\hat{H}$ have exact eigenstates $|j⟩$ with eigenvalues $E_j^{\text{exact}}$ in non-decreasing order. The variational subspace $\mathcal{V}=\mathrm{span}\{{\varphi }_1,\dots ,{\varphi }_N\}$ has dimension $N$. Consider the subspace $\mathcal{W}$ spanned by the exact eigenstates $|0⟩,|1⟩,\dots ,|k⟩$. Since $\dim \mathcal{W}=k+1$, the intersection $\mathcal{V}\cap {\mathcal{W}}^{\perp }$ has dimension at least $N-k$ (where ${\mathcal{W}}^{\perp }$ is the orthogonal complement of $\mathcal{W}$ in $\mathcal{V}$).

For any $|\psi ⟩\in \mathcal{V}\cap {\mathcal{W}}^{\perp }$, the expansion $|\psi ⟩={\sum }_{j\ge k+1}{a}_j|j⟩$ contains only terms with $j\ge k+1$, so

$$⟨\psi |\hat{H}|\psi ⟩=\sum _{j\ge k+1}|a_j{|}^2{E}_j^{\text{exact}}\ge E_k^{\text{exact}}\sum _{j\ge k+1}|a_j{|}^2=E_k^{\text{exact}}.$$

The minimum of $⟨\psi |\hat{H}|\psi ⟩$ over all normalised $|\psi ⟩$ in $\mathcal{V}$, subject to being orthogonal to the first $k$ variational eigenvectors, equals $E_k$ by the Rayleigh-Ritz variational characterisation. Since $\mathcal{V}\cap {\mathcal{W}}^{\perp }$ is a subset of the admissible states, $E_k\ge E_k^{\text{exact}}$. $\square$

Corollary (Bracketing theorem). If the variational eigenvalues $E_0,E_1,\dots ,E_{N-1}$ satisfy $E_k<E_{k+1}$ and there exists a gap between exact eigenvalues such that $E_k^{\text{exact}}<E_{k+1}^{\text{exact}}$, then $E_k^{\text{exact}}\le E_k<E_{k+1}^{\text{exact}}$ is possible but $E_k^{\text{exact}}\le E_k\le E_{k+1}^{\text{exact}}$ always holds.

Worked example at intermediate level

Construct and solve the secular determinant for a particle in a one-dimensional box ($V=0$ for $0<x<L$, $V=\infty$ otherwise) using the two-function basis ${\varphi }_1=x(L-x)$ and ${\varphi }_2=x^2(L-x)$.

Both functions vanish at $x=0$ and $x=L$, satisfying the boundary conditions. The Hamiltonian is $\hat{H}=-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}$.

The overlap matrix element is $S_{12}={\int }_0^{L}x(L-x)\cdot x^2(L-x)dx=L^6/60$.

The Hamiltonian matrix elements involve the second derivatives: $\hat{H}{\varphi }_1=-\frac{{\hslash }^2}{2m}(-2)=\frac{{\hslash }^2}{m}$ and $\hat{H}{\varphi }_2=-\frac{{\hslash }^2}{2m}(2L-6x)=\frac{{\hslash }^2}{2m}(6x-2L)$.

Computing $H_{11}=⟨{\varphi }_1|\hat{H}|{\varphi }_1⟩=\frac{{\hslash }^2}{m}{\int }_0^{L}x(L-x)dx=\frac{{\hslash }^2{L}^3}{6m}$. Similarly, $H_{12}=\frac{{\hslash }^2{L}^4}{10m}$ and $S_{11}=\frac{L^5}{30}$.

Setting up the $2\times 2$ secular equation and solving gives the lowest variational energy $E_0=0.1250{\hslash }^2{\pi }^2/(2m{L}^2)$, compared with the exact $E_0^{\text{exact}}={\pi }^2{\hslash }^2/(2m{L}^2)$. The variational result exceeds the exact by $E_0/E_0^{\text{exact}}=1.267$, a 26.7% overestimate. This is expected for a small, unoptimised basis.

Exercises Intermediate+

The Roothaan equations and computational cost Master

For closed-shell molecules, the Hartree-Fock equations $\hat{f}{\chi }_a={\epsilon }_a{\chi }_a$ are a set of coupled integro-differential equations with no closed-form solution. Roothaan (1951) and Hall (1951) independently showed that expanding each molecular orbital in a finite set of $K$ basis functions $\{{\varphi }_{\mu }\}$ converts the Hartree-Fock equations into a matrix algebraic equation. Writing ${\psi }_a={\sum }_{\mu =1}^K{C}_{\mu a}{\varphi }_{\mu }$, the Fock equation becomes the Roothaan equation

$$\mathbf{F}\mathbf{C}=\mathbf{S}\mathbf{C}\mathbf{\epsilon },$$

where $F_{\mu \nu }=⟨{\varphi }_{\mu }|\hat{f}|{\varphi }_{\nu }⟩$ is the Fock matrix, $S_{\mu \nu }=⟨{\varphi }_{\mu }|{\varphi }_{\nu }⟩$ is the overlap matrix, $\mathbf{C}$ is the matrix of orbital coefficients, and $\mathbf{\epsilon }$ is the diagonal matrix of orbital energies. This is a generalised eigenvalue problem, solved in each SCF iteration.

The Fock matrix elements depend on two-electron integrals over the basis functions:

$$(\mu \nu |\lambda \sigma )=\iint {\varphi }_{\mu }({\mathbf{r}}_1){\varphi }_{\nu }({\mathbf{r}}_1)\frac{1}{r_{12}}{\varphi }_{\lambda }({\mathbf{r}}_2){\varphi }_{\sigma }({\mathbf{r}}_2)d{\mathbf{r}}_{1}d{\mathbf{r}}_2.$$

The number of unique two-electron integrals (using all permutational symmetry) scales as $O(K^4/8)$ for $K$ basis functions. For a moderate-sized molecule with $K=500$ basis functions (roughly a 30-atom system at the 6-31G(d) level), this gives approximately $5.2\times 10^9$ integrals. At 8 bytes per double-precision number, this requires about 40 GB of storage. This $O(K^4)$ scaling was the primary bottleneck of quantum chemistry for decades.

Modern implementations avoid storing all integrals. The resolution-of-the-identity (RI) approximation replaces the four-centre integral with a product of three-centre integrals expanded in an auxiliary basis, reducing the cost to $O(K^3)$. Density fitting is the most common form of RI. For large systems, linear-scaling methods exploit the locality of Gaussian functions: integrals involving distant basis functions are negligibly small and can be neglected, reducing the effective scaling to $O(K)$ for extended systems.

The SCF cycle itself converges the density matrix $\mathbf{P}$ through iterative updates. The convergence is not guaranteed to be monotonic and can oscillate. Damping (mixing old and new density matrices), level shifting (raising virtual orbital energies to reduce mixing), and DIIS (direct inversion in the iterative subspace, Pulay 1980) are standard acceleration techniques. DIIS extrapolates from previous iterations to predict the converged solution, typically achieving convergence in 6--15 iterations for well-behaved systems.

Formal scaling of quantum chemical methods

Method	Formal scaling	Typical system size limit
Hartree-Fock (HF)	$O(K^4)$	~500 atoms (with RI)
MP2	$O(K^5)$	~100 atoms
CCSD	$O(K^6)$	~30 atoms
CCSD(T)	$O(K^7)$	~20 atoms
Full CI	$O(e^K)$	~10 electrons, small basis

These scalings are the reason DFT ($O(K^3)$ with RI, or $O(K)$ with local methods) dominates routine computational chemistry despite its approximations.

Basis set convergence and the complete basis set limit Master

The variational principle guarantees that the Hartree-Fock energy decreases monotonically as the basis set grows. The complete basis set (CBS) limit is the energy obtained in the limit $K\to \infty$. In practice, basis set convergence is systematic and can be extrapolated.

Dunning's correlation-consistent basis sets (cc-pV$n$Z, where $n=$ D, T, Q, 5, 6) are designed for systematic convergence of the correlation energy. The cardinal number $n$ controls the flexibility: cc-pVDZ (double-zeta), cc-pVTZ (triple-zeta), cc-pVQZ (quadruple-zeta), cc-pV5Z, cc-pV6Z. Each step adds a full shell of higher-angular-momentum functions, and the correlation energy converges approximately as $E_{\text{corr}}(n)\approx E_{\text{corr}}^{\text{CBS}}+A/n^3$ (Helgaker et al., 1997).

For the Hartree-Fock energy, the convergence is exponential in $n$ (much faster). The slow convergence of the correlation energy dominates computational cost, because each additional angular momentum shell roughly triples the basis set size while increasing the integral count by a factor of $3^4=81$.

Basis set superposition error (BSSE) arises when two molecules interact: each monomer "borrows" basis functions from its partner, artificially lowering its energy. The counterpoise correction (Boys and Bernardi, 1970) computes each monomer in the full dimer basis to remove this artifact. BSSE decreases with basis set quality and vanishes at the CBS limit.

The practical hierarchy for basis set selection follows the problem. For geometry optimisations and vibrational frequencies, 6-31G(d) or def2-SVP suffices. For accurate energies, cc-pVTZ or def2-TZVP is the minimum. For benchmark-quality thermochemistry, cc-pVQZ or cc-pV5Z with CBS extrapolation is standard. The choice of basis set is often the single most important decision in a quantum chemical calculation.

The Roothaan equations for H${}_2^{+}$ Master

The simplest molecular application of the Roothaan formalism is H${}_2^{+}$, which has only one electron. The Fock operator reduces to the core Hamiltonian (no Coulomb or exchange terms), so the Roothaan equation becomes a simple $K\times K$ generalised eigenvalue problem. With a minimal basis of two 1s STOs:

$$\left(\begin{array}{cc}{H}_{AA}& H_{AB}\\ H_{BA}& H_{BB}\end{array}\right)\left(\begin{array}{c}{c}_A\\ c_B\end{array}\right)=E\left(\begin{array}{cc}1& S\\ S& 1\end{array}\right)\left(\begin{array}{c}{c}_A\\ c_B\end{array}\right),$$

where $H_{AA}=H_{BB}=\alpha$ (by symmetry of the two identical nuclei). The solutions are the bonding (${\sigma }_g$) and antibonding (${\sigma }_u$) combinations derived in the beginner worked example.

Enlarging the basis to include 2$s$ and 2$p_{\sigma }$ functions on each centre gives a $6\times 6$ problem. The additional flexibility allows the bonding orbital to contract and polarise toward the internuclear region, recovering more of the binding energy. With this extended basis, the calculated equilibrium bond length and dissociation energy of H${}_2^{+}$ approach the exact values ($R_e=1.997a_0$, $D_e=0.1026$ hartree) to within chemical accuracy.

The convergence pattern is instructive. The minimal (two 1s) basis recovers about 60% of the binding energy. Adding 2$p_{\sigma }$ functions (four-function basis) recovers about 90%. A full triple-zeta plus polarisation basis recovers over 99%. The remaining 1% requires high-angular-momentum functions that describe the cusp and polarisation with increasing fidelity.

Connections Master

• Many-electron atoms: Hartree-Fock 14.01.03 develops the atomic Hartree-Fock equations, the Fock operator, and the self-consistent field cycle. The present unit extends those ideas to molecules by introducing basis set expansions and the Roothaan equations. The Fock operator structure (core Hamiltonian + Coulomb + exchange) is the same; only the basis representation changes.

• Hydrogen atom quantum chemistry 14.04.01 provides the 1s, 2s, 2p atomic orbitals that serve as the basis functions in minimal-basis molecular calculations. The shapes and energies of hydrogenic orbitals determine the overlap and resonance integrals that drive chemical bonding.

• Molecular orbital theory 14.05.02 presents the qualitative MO picture (bonding, antibonding, $\sigma$, $\pi$) that emerges naturally from the linear variational method. The present unit provides the quantitative machinery that underlies those qualitative arguments.

• Perturbation theory 12.07.01 (physics sequence) provides the mathematical framework for Moller-Plesset perturbation theory (MP2, MP3, MP4), which treats electron correlation as a perturbation to the Hartree-Fock reference. The basis set machinery developed here is shared between HF and all post-HF methods.

Historical & philosophical context Master

The variational principle in quantum mechanics was first applied by Schrodinger himself in 1926, who used it to derive the wave equation from an energy minimisation principle. The linear variational method (Rayleigh-Ritz) has older roots in classical mechanics and the theory of vibrations. Rayleigh (1877) used it to find approximate vibration frequencies of continuous systems, and Ritz (1909) formalised the procedure for arbitrary linear expansions.

The secular determinant takes its name from the "secular equation" of celestial mechanics, where it arose in the study of long-period (secular) perturbations of planetary orbits. The mathematical structure -- a determinant set to zero whose roots are eigenvalues -- is identical, whether the context is planetary dynamics or molecular orbitals.

Roothaan's 1951 paper ("New Developments in Molecular Orbital Theory", Rev. Mod. Phys. 23, 69--89) transformed quantum chemistry from a field of clever hand calculations to one amenable to systematic computation. By recasting the Hartree-Fock equations as a matrix eigenvalue problem, Roothaan made it possible to program the entire SCF procedure on the electronic computers that were just becoming available. Hall's independent derivation appeared the same year. The Roothaan-Hall equations remain the foundation of every Hartree-Fock and Kohn-Sham DFT program in use today.

The development of Gaussian basis sets followed a parallel trajectory. Boys (1950) first proposed using Gaussian functions, recognising their computational advantages. The price -- poor representation of the nuclear cusp and exponential tail -- was judged acceptable given the speed of integral evaluation. The contracted Gaussian basis sets of Pople and coworkers (STO-$n$G, 3-21G, 6-31G, 6-311G with various polarisation and diffuse extensions) became the standard for routine calculations from the 1970s onward. Dunning's correlation-consistent basis sets (cc-pV$n$Z), introduced in 1989, were designed for systematic convergence of electron correlation and are now the standard for high-accuracy work.

The philosophical significance of the variational principle lies in its relationship to exact solutions. The variational energy is always an upper bound, but there is no general lower bound. A calculation yielding $E=-76.02$ hartree for water tells you the true energy is below this value, but not by how much. The convergence of the variational energy with increasing basis set size provides empirical evidence for approach to the exact result, but convergence itself is not a proof of correctness. This asymmetry -- firm upper bounds but no firm lower bounds -- is a structural feature of the variational method that distinguishes it from methods with two-sided error bounds.

Bibliography Master

• Roothaan, C. C. J., "New Developments in Molecular Orbital Theory", Rev. Mod. Phys. 23 (1951), 69--89.
• Hall, G. G., "The Molecular Orbital Theory of Chemical Valency VIII. A Method of Calculating Ionization Potentials", Proc. R. Soc. Lond. A 205 (1951), 541--552.
• Boys, S. F., "Electronic Wave Functions I. A General Method of Calculation for the Stationary States of Any Molecular System", Proc. R. Soc. Lond. A 200 (1950), 542--554.
• Hehre, W. J., Stewart, R. F., & Pople, J. A., "Self-Consistent Molecular-Orbital Methods I. Use of Gaussian Expansions for Molecular Orbitals", J. Chem. Phys. 51 (1969), 2657--2664.
• Dunning, T. H., "Gaussian Basis Sets for Use in Correlated Molecular Calculations I. The Atoms Boron Through Neon and Hydrogen", J. Chem. Phys. 90 (1989), 1007--1023.
• Helgaker, T., Klopper, W., Koch, H., & Noga, J., "Basis-Set Convergence of Correlated Calculations on Water", J. Chem. Phys. 106 (1997), 9639--9646.
• Pulay, P., "Convergence Acceleration of Iterative Sequences. The Case of SCF Iteration", Chem. Phys. Lett. 73 (1980), 393--398.
• Szabo, A. & Ostlund, N. S., Modern Quantum Chemistry (Dover, 1996), Ch. 3.
• Levine, I. N., Quantum Chemistry, 7e (Pearson, 2014), Ch. 8.
• Helgaker, T., Jorgensen, P., & Olsen, J., Molecular Electronic-Structure Theory (Wiley, 2000), Ch. 3.
• MacDonald, J. K. L., "Successive Approximations by the Rayleigh-Ritz Variation Method", Phys. Rev. 43 (1933), 830--833.
• Hylleraas, E. A. & Undheim, B., "Numerische Berechnung der 2S-Terme von Ortho- und Par-Helium", Z. Physik 65 (1930), 759--772.