The Variational Method with Trial Functions

Guess a wavefunction, tune its knobs, and watch the energy fall. It falls as low as your guess allows. This is the entire craft of the variational method with trial functions. It sounds almost too simple to be a cornerstone of quantum chemistry. Yet this single idea underlies how chemists estimate ionization energies of atoms nobody can solve exactly. It explains how the helium atom yields to pencil-and-paper analysis. It also shows how Hartree–Fock and density functional theory ultimately get their orbitals. The entire enterprise rests on The Variational Principle. For any normalized trial function, the expectation value of the Hamiltonian never falls below the true ground-state energy. That inequality converts an intractable differential equation into a tractable calculus problem. You minimize a function of a few parameters instead of solving the Schrödinger equation exactly. The physical motivation is concrete. Most atoms and molecules have no closed-form solution to the Schrödinger equation. This is because of electron–electron repulsion. Yet chemists routinely need quantitative energies to explain trends in ionization potentials, bond strengths, and reaction barriers. The variational method lets you build a physically motivated ansatz. It might be a screened hydrogenic orbital, a Gaussian, or a linear combination of atomic orbitals. You encode the chemistry you already understand into a handful of adjustable parameters. This includes nuclear charge, screening, and nodal structure. Then you let calculus do the rest. The output is not just a number. The value of the optimized parameter itself tells a physical story. For example, it may show an effective nuclear charge reduced by electron shielding. This article develops the method in full generality. It then covers two archetypal settings every quantum chemistry student must master. The first is a single continuously adjustable parameter, illustrated by the helium atom's effective nuclear charge. The second is a scaling parameter applied to a fixed trial shape. This is illustrated by the hydrogen atom, where the exact answer lets you check the machinery. Along the way you will see exactly how the Rayleigh quotient is built from an unnormalized trial function. You will see how differentiation under the integral sign produces the stationarity condition. You will also see how the second derivative test confirms you have found a minimum rather than a saddle point. By the end, you will be able to set up the variational integral for any parameterized trial function. You will carry the minimization through explicitly. You will interpret the optimized parameter chemically. And you will know precisely when the method is trustworthy and when it silently misleads. This groundwork feeds directly into Optimizing Variational Parameters. There, multi-parameter and multi-determinant trial functions are tuned simultaneously. It also feeds into Linear Variation and Secular Equations. There, the parameters enter linearly. The minimization becomes a matrix eigenvalue problem instead of ordinary calculus.

Key takeaways

Pick a normalizable trial function $\phi(\mathbf{r};\alpha)$ containing one or more adjustable parameters $\alpha$ .
Build the Rayleigh quotient $E(\alpha) = \langle\phi|\hat H|\phi\rangle/\langle\phi|\phi\rangle$ ; this always satisfies $E(\alpha) \ge E_0$ .
Minimize $E(\alpha)$ over $\alpha$ by solving $\partial E/\partial \alpha = 0$ and checking $\partial^2 E/\partial \alpha^2 > 0$ .
The optimal $\alpha^*$ often has direct physical meaning: an effective nuclear charge, an effective width, an effective bond order.
The method never gives the exact ground energy unless the trial function happens to contain the true eigenfunction as a special case.
Excited-state estimates require enforcing orthogonality to lower states, or using symmetry to guarantee it automatically.
Key idea: minimizing the energy of a parameterized guess over its parameters turns the variational principle into an ordinary calculus problem whose solution bounds — and often closely approximates — the true ground-state energy.

Where this lives today

The variational method with parameterized trial functions is not a historical curiosity; it is the conceptual DNA of essentially all modern quantum chemistry software. Hartree–Fock theory is itself a variational method in which the trial function is a single Slater determinant and the "parameters" are the orbital coefficients — Optimizing Variational Parameters extends the single-knob calculus shown here into the many-parameter self-consistent-field procedure used in every electronic structure package (Gaussian, ORCA, Psi4, Q-Chem). Density functional theory codes likewise build molecular orbitals from linear combinations of basis functions whose coefficients are variationally optimized, following exactly the logic of Linear Variation and Secular Equations once the parameters enter linearly rather than as exponents. In drug design, variationally optimized atomic orbital exponents (Slater-type or contracted Gaussian) underlie the basis sets — STO-3G, 6-31G*, cc-pVTZ — used to compute binding energies between candidate drug molecules and protein targets; every one of those basis functions was itself variationally fitted, historically by matching Gaussian combinations to the shape of accurately known atomic orbitals. In materials science, variational trial functions with scaling parameters describe electron localization in defects and color centers, informing the design of phosphors and semiconductor dopants. Quantum Monte Carlo methods used for benchmark-quality energies in condensed matter physics use explicitly correlated Jastrow-type trial functions — direct intellectual descendants of Hylleraas's $r_{12}$ -dependent terms — with dozens of variational parameters optimized by stochastic energy minimization on supercomputers. Even the pedagogical helium screening constant $\zeta = 27/16$ derived above reappears, in spirit, every time a chemist rationalizes ionization energy trends across the periodic table using Slater's rules for effective nuclear charge — those empirical screening constants are themselves rough variational estimates systematized into a lookup table.

Ask the tutor: How would the helium screening calculation change if you allowed the two electrons to have different effective nuclear charges instead of a shared $\zeta$ ?

The physical picture

The variational method starts from a physical intuition rather than a mathematical trick: if you already have a rough idea of what the ground-state wavefunction should look like — where it peaks, how fast it decays, its nodal structure — you can write down a functional form that captures this shape while leaving one or more parameters open. The parameters represent physical effects you have not yet quantified. In helium, the obvious first guess for each electron is a hydrogen-like 1s orbital, but each electron feels a nuclear charge $Z=2$ partially screened by the other electron. Rather than solving the coupled electron–electron problem exactly, you introduce an effective nuclear charge $\zeta$ as a free parameter and let the variational principle determine the best value from first principles — no fitting to experimental data required. Once the trial function $\phi(\mathbf{r};\alpha)$ is chosen, the machinery is completely mechanical. You compute the Rayleigh quotient,

E(\alpha) = \frac{\displaystyle\int \phi^*(\mathbf{r};\alpha)\,\hat H\,\phi(\mathbf{r};\alpha)\,d\tau}{\displaystyle\int \phi^*(\mathbf{r};\alpha)\,\phi(\mathbf{r};\alpha)\,d\tau},

as an explicit function of $\alpha$ , using whatever normalization is convenient (the quotient automatically compensates for an unnormalized $\phi$ ). Then you differentiate with respect to $\alpha$ , set the derivative to zero, and solve for the stationary point $\alpha^*$ . The variational principle guarantees $E(\alpha) \ge E_0$ for every $\alpha$ , so the minimum over $\alpha$ is the best possible upper bound obtainable from that particular functional form. Crucially, $E(\alpha^*)$ is not merely a number: it is the tightest bound your chosen shape can produce, and the gap $E(\alpha^*) - E_0$ measures how well your functional form matches physical reality, independent of the specific parameter values. Geometrically, think of $E(\alpha)$ as a landscape sitting entirely on or above the horizontal line $E = E_0$ . The trial family sweeps out a curve as $\alpha$ varies, dipping down to touch the true ground energy only if the family happens to contain the exact ground-state wavefunction. Most families never touch that floor, but a well-chosen family gets remarkably close — often within a few percent of the exact energy, using no more than one or two parameters. This is the central bargain of the variational method: physical insight, encoded in the shape of the trial function, converts an unsolvable differential equation into an elementary problem in single- or multi-variable calculus.

Figure

The variational landscape

E(\alpha)

and the exact floor

E_0

Figure

Effective nuclear charge picture for He: screening reshapes the 1s orbital

The core ideas

The single equation that defines the whole method is the Rayleigh quotient built from a trial function:

E(\alpha) = \frac{\langle \phi(\alpha)|\hat H|\phi(\alpha)\rangle}{\langle \phi(\alpha)|\phi(\alpha)\rangle} \ge E_0 .

where: $\phi(\alpha)$ is any normalizable trial function depending on parameter(s) $\alpha$ ; $\hat H$ is the exact molecular or atomic Hamiltonian; $E_0$ is the true ground-state energy; the inequality holds for every value of $\alpha$ , not only the optimum. Reading: this is The Variational Principle applied at a single point $\alpha$ in parameter space — it says nothing yet about optimization, only that any guess overestimates (or exactly equals) the true energy. The optimization step is the ordinary stationarity condition of single- or multi-variable calculus:

\left.\frac{\partial E(\alpha)}{\partial \alpha}\right|_{\alpha=\alpha^*} = 0, \qquad \left.\frac{\partial^2 E(\alpha)}{\partial \alpha^2}\right|_{\alpha=\alpha^*} > 0 .

where: $\alpha^*$ is the optimal parameter value; the first condition locates stationary points; the second condition confirms a genuine minimum (not a maximum or inflection) among the stationary points. Reading: because $E(\alpha)$ is bounded below by $E_0$ , any interior minimum found this way is the best upper bound attainable from that trial family — you are sliding down the variational landscape until it flattens out. For a trial function built by scaling a fixed reference shape, $\phi(\mathbf{r};\alpha) = \chi(\alpha \mathbf{r})$ , there is a shortcut that avoids recomputing integrals from scratch: the scaling theorem,

E(\alpha) = \alpha^2 \langle \chi|\hat T|\chi\rangle + \alpha \langle \chi|\hat V|\chi\rangle ,

where: $\hat T$ is the kinetic energy operator and $\hat V$ the potential energy operator evaluated with the unscaled reference function $\chi$ ; the kinetic contribution scales as $\alpha^2$ because differentiation twice with respect to scaled coordinates brings down two powers of $\alpha$ ; the potential (Coulombic, $\propto 1/r$ ) contributes one power of $\alpha$ from the $1/r \to \alpha/r$ substitution combined with the volume element. Reading: this quadratic-in- $\alpha$ , linear-in- $\alpha$ split is why scaling problems (like the hydrogen 1s orbital exponent, or exploring how the electron cloud "breathes" under scaling) always minimize to a clean closed form — you differentiate a simple quadratic, not a transcendental function.

Derivation AOptimizing the screening parameter ζ for the helium ground state

Goal. Using a product of hydrogen-like 1s orbitals with an adjustable effective nuclear charge $\zeta$ as trial function for the helium atom's two electrons, find $\zeta^*$ that minimizes the variational energy, and evaluate the resulting bound.

Step 1 — Write the Hamiltonian and trial function. In atomic units, the exact electronic Hamiltonian for helium ( $Z=2$ ) is

\hat H = -\tfrac{1}{2}\nabla_1^2 - \tfrac{1}{2}\nabla_2^2 - \frac{Z}{r_1} - \frac{Z}{r_2} + \frac{1}{r_{12}} .

The trial function is a product of two hydrogenic 1s orbitals sharing a common variational parameter $\zeta$ (the effective nuclear charge, not fixed at 2):

\phi(\mathbf{r}_1,\mathbf{r}_2;\zeta) = \psi_{1s}(\mathbf{r}_1;\zeta)\,\psi_{1s}(\mathbf{r}_2;\zeta), \qquad \psi_{1s}(\mathbf{r};\zeta) = \left(\frac{\zeta^3}{\pi}\right)^{1/2} e^{-\zeta r}.

This function is already normalized because each hydrogenic 1s orbital is separately normalized.

Step 2 — Split the Hamiltonian into a solvable piece and a perturbation-like piece. Rewrite $\hat H$ by adding and subtracting $\zeta/r_1$ and $\zeta/r_2$ :

\hat H = \left[-\tfrac{1}{2}\nabla_1^2 - \frac{\zeta}{r_1}\right] + \left[-\tfrac{1}{2}\nabla_2^2 - \frac{\zeta}{r_2}\right] + \frac{\zeta - Z}{r_1} + \frac{\zeta - Z}{r_2} + \frac{1}{r_{12}} .

The bracketed terms are exactly the hydrogen-like Hamiltonians with nuclear charge $\zeta$ , whose eigenvalue on $\psi_{1s}(\mathbf{r};\zeta)$ is $-\zeta^2/2$ (in hartrees, since the 1s hydrogenic energy is $-\zeta^2/2n^2$ with $n=1$ ). This algebraic trick avoids recomputing the kinetic and nuclear-attraction integrals from scratch — you recognize them as scaled hydrogen atom integrals.

Step 3 — Evaluate each expectation value. Because $\phi$ is a simple product, expectation values separate electron-by-electron except for the $1/r_{12}$ term:

\left\langle -\tfrac{1}{2}\nabla_i^2 - \frac{\zeta}{r_i}\right\rangle = -\frac{\zeta^2}{2}, \qquad \left\langle \frac{1}{r_i}\right\rangle_{\psi_{1s}(\zeta)} = \zeta .

The first result is the exact hydrogenic 1s energy at charge $\zeta$ ; the second is the standard hydrogenic expectation value $\langle 1/r\rangle = \zeta/n^2 = \zeta$ for $n=1$ . The electron repulsion integral for two 1s orbitals of exponent $\zeta$ is a classic result,

\left\langle \frac{1}{r_{12}}\right\rangle = \frac{5}{8}\zeta .

Step 4 — Assemble $E(\zeta)$ . Summing all pieces (two hydrogenic kinetic-plus-nuclear terms, two screening corrections, one repulsion term):

E(\zeta) = 2\left(-\frac{\zeta^2}{2}\right) + 2(\zeta - Z)\zeta + \frac{5}{8}\zeta = -\zeta^2 + 2\zeta^2 - 2Z\zeta + \frac{5}{8}\zeta .

Collecting terms,

E(\zeta) = \zeta^2 - 2Z\zeta + \frac{5}{8}\zeta = \zeta^2 - \left(2Z - \frac{5}{8}\right)\zeta .

Step 5 — Minimize over $\zeta$ . Differentiate and set to zero:

\frac{dE}{d\zeta} = 2\zeta - \left(2Z - \frac{5}{8}\right) = 0 \quad\Longrightarrow\quad \zeta^* = Z - \frac{5}{16}.

The second derivative is $d^2E/d\zeta^2 = 2 > 0$ , confirming a genuine minimum. For helium, $Z = 2$ , so

\zeta^* = 2 - \frac{5}{16} = \frac{27}{16} = 1.6875 .

Step 6 — Evaluate the optimized energy. Substitute $\zeta^*= Z - 5/16$ back into $E(\zeta) = \zeta^2 - (2Z - 5/8)\zeta$ :

E(\zeta^*) = -\left(Z - \frac{5}{16}\right)^2 .

With $Z=2$ : $E(\zeta^*) = -(27/16)^2 = -\dfrac{729}{256} = -2.84766$ hartree. This compares to the experimental helium ground-state energy of $-2.90372$ hartree — a variational bound that is $2\%$ too high, exactly as the variational principle requires ( $E(\zeta^*) > E_0$ ), and dramatically better than the crude $Z=2$ unscreened guess, which gives $E = -Z^2 = -4.0$ hartree for the kinetic-plus-nuclear part alone before repulsion is even included, or more precisely $E(\zeta{=}2) = 4 - (2\cdot2 - 5/8)\cdot2 = 4 - 5.75 = -2.75$ hartree, itself already an upper bound but a noticeably looser one than $\zeta^*=27/16$ gives.

Go deeper · Derivation B · 2 minScaling a trial hydrogen 1s orbital and recovering the exact answer

Derivation BScaling a trial hydrogen 1s orbital and recovering the exact answer

Goal. Take the shape of the hydrogen 1s orbital, $\chi(r) = e^{-r}$ , treat its exponent as a free scaling parameter $\alpha$ rather than fixing it at 1, and show that variational optimization exactly reproduces the true hydrogen ground state — the crucial sanity check that validates the whole machinery.

Step 1 — Define the scaled trial function. Let

\phi(r;\alpha) = \left(\frac{\alpha^3}{\pi}\right)^{1/2} e^{-\alpha r},

which is normalized for every $\alpha > 0$ (the same normalization integral that makes hydrogenic 1s orbitals normalized at any nuclear charge). The hydrogen atom Hamiltonian in atomic units is

\hat H = -\frac{1}{2}\nabla^2 - \frac{1}{r}.

Step 2 — Apply the scaling theorem. Writing $\phi(\mathbf{r};\alpha) = \chi(\alpha\mathbf{r})$ with $\chi(\mathbf{r}) = (1/\pi)^{1/2}e^{-r}$ suitably normalized, the general scaling result gives

E(\alpha) = \alpha^2 T_1 - \alpha V_1,

where $T_1 = \langle\chi|-\tfrac12\nabla^2|\chi\rangle$ and $V_1 = \langle\chi|1/r|\chi\rangle$ are the kinetic and (magnitude of) potential expectation values at $\alpha=1$ , i.e. for the unscaled reference function. For the unscaled 1s shape ( $\alpha=1$ is literally the hydrogenic solution), both integrals are elementary:

T_1 = \frac{1}{2}, \qquad V_1 = 1 .

(These follow from direct integration of $-\tfrac12\nabla^2$ and $1/r$ against $e^{-r}$ with the standard normalization; they are the same integrals that give the exact hydrogen ground energy $-1/2$ hartree when $\alpha=1$ is substituted, serving as an internal consistency check.)

Step 3 — Assemble $E(\alpha)$ . Substituting,

E(\alpha) = \frac{1}{2}\alpha^2 - \alpha .

This is an exact quadratic in $\alpha$ — no approximation beyond the choice of trial shape has been made yet.

Step 4 — Minimize. Differentiate:

\frac{dE}{d\alpha} = \alpha - 1 = 0 \quad\Longrightarrow\quad \alpha^* = 1.

The second derivative $d^2E/d\alpha^2 = 1 > 0$ confirms a minimum. Substituting back,

E(\alpha^*) = \frac{1}{2}(1)^2 - 1 = -\frac{1}{2}\ \text{hartree}.

Step 5 — Interpret. The variational optimum $\alpha^*=1$ reproduces the exact hydrogen 1s exponent, and $E(\alpha^*) = -0.5$ hartree exactly matches the known ground-state energy, $-13.6$ eV. This is not a coincidence: because the trial family itself contains the true eigenfunction (at $\alpha=1$ , $\phi$ is exactly the hydrogen 1s orbital), the variational bound becomes an equality. This example is the benchmark that proves the calculus machinery is being executed correctly — any trial family flexible enough to contain the true wavefunction as a special case will always relax to it under minimization, and no residual error remains.

Step 6 — Limiting-case check. As a stress test, consider what happens if the reference shape is not exact, e.g. a Gaussian $\chi(r) = e^{-r^2}$ instead of the exponential. Repeating Steps 2–4 with Gaussian integrals gives a different optimal exponent and an energy of $E \approx -0.424$ hartree — worse than $-0.5$ hartree, as it must be, since the Gaussian shape never coincides with the true exponential 1s orbital no matter how its width is tuned. This illustrates that the quality of the bound is set by the shape, while the scaling parameter only optimizes within that shape.

Worked examples

Example 1 — Foundational: a particle in a box with a parabolic trial function. A particle of mass $m$ is confined to an infinite square well of width $L$ , $0 \le x \le L$ . Instead of using the known exact solution $\sin(\pi x/L)$ , use the simplest polynomial trial function satisfying the boundary conditions $\phi(0)=\phi(L)=0$ :

\phi(x;\alpha) = x^\alpha (L-x)^\alpha ,

but to keep the algebra elementary, fix the functional form to the lowest nontrivial polynomial and instead vary a different parameter — a shift in normalization exponent is unnecessary here, so use the standard textbook choice

\phi(x) = x(L-x)

with no free parameter, and compare against the trial function with one adjustable parameter $c$ ,

\phi(x;c) = x(L-x) + c\,x^2(L-x)^2 .

For a clean, fully asymmetric worked calculation, set $c=0$ first (the parameter-free case) to establish the baseline, then show how introducing $c$ improves it. Baseline ( $c=0$ ). Normalize $\phi(x) = x(L-x)$ :

\int_0^L x^2(L-x)^2\,dx = \frac{L^5}{30}.

The kinetic energy integral (since $V=0$ inside the box):

\left\langle -\frac{\hbar^2}{2m}\frac{d^2}{dx^2}\right\rangle = \frac{\hbar^2}{2m}\int_0^L \left(\phi'(x)\right)^2 dx \Big/ \int_0^L \phi(x)^2\,dx.

With $\phi'(x) = L - 2x$ , $\int_0^L (L-2x)^2\,dx = L^3/3$ . Thus

E = \frac{\hbar^2}{2m}\cdot\frac{L^3/3}{L^5/30} = \frac{\hbar^2}{2m}\cdot\frac{10}{L^2} = \frac{5\hbar^2}{mL^2}.

Converting to the usual form $E = \dfrac{\pi^2\hbar^2}{2mL^2}\times(\text{factor})$ : the exact ground energy is $E_0 = \pi^2\hbar^2/(2mL^2) \approx 4.9348\,\hbar^2/(mL^2)$ , while the trial estimate gives $5\,\hbar^2/(mL^2)$ — higher by only $1.3\%$ , an excellent bound from an utterly elementary polynomial guess containing no adjustable parameter at all. Adding the parameter $c$ . This baseline already demonstrates the power of shape; introducing $c$ and minimizing $E(c)$ (a short but genuine two-term quadratic-over-quadratic optimization, entirely analogous to the derivations above) pushes the bound to $4.9348\,\hbar^2/(mL^2)$ to five significant figures — indistinguishable from the exact answer, since the quartic correction $x^2(L-x)^2$ supplies exactly the extra curvature the parabola lacks near the walls. The lesson: even one extra parameter, when it captures missing shape near the nodes, can close almost the entire remaining error. Limiting-case check. As $L \to \infty$ (box becomes free space), both $E_0$ and the trial estimate scale as $1/L^2 \to 0$ , correctly recovering the free-particle continuum limit where confinement energy vanishes — the ratio of trial to exact energy stays fixed at $\approx 1.013$ regardless of $L$ , since $L$ is an overall energy scale, not a shape parameter.

Example 2 — Applied: the variationally optimized exponent for the hydrogen molecular ion's LCAO orbital. For $\mathrm{H_2^+}$ , a natural trial function is the linear combination of atomic orbitals (LCAO) built from hydrogenic 1s functions on each nucleus, but with the orbital exponent $\zeta$ left variational rather than fixed at the free-atom value $\zeta=1$ :

\phi(\mathbf{r};\zeta) = N\left[\,\psi_{1s}(\mathbf{r}_A;\zeta) + \psi_{1s}(\mathbf{r}_B;\zeta)\,\right],

where $\mathbf{r}_A, \mathbf{r}_B$ are electron–nucleus distances to centers $A$ and $B$ separated by the (fixed, or separately optimized) bond distance $R$ . Physically, $\zeta > 1$ is expected: the single electron feels two nuclei, so the optimal atomic-like orbital should be more contracted than the free hydrogen atom's 1s orbital ( $\zeta=1$ ), because the shared electron sits, on average, in a region of enhanced positive potential from both centers. Carrying out the full two-center integral evaluation (overlap $S(\zeta,R)$ , Coulomb $J(\zeta,R)$ , resonance $K(\zeta,R)$ integrals — all textbook two-center integrals expressible via confocal elliptical coordinates) and minimizing $E(\zeta,R)$ simultaneously over both $\zeta$ and $R$ yields, at the variational optimum,

\zeta^* \approx 1.24, \qquad R^* \approx 2.00\ a_0 ,

with a bond dissociation energy of $D_e \approx 2.35$ eV, compared to the experimental $D_e = 2.79$ eV. Freezing $\zeta = 1$ (the naive free-atom exponent, no optimization) gives a noticeably worse $D_e \approx 1.76$ eV at $R^* \approx 2.5\,a_0$ . The improvement — nearly 0.6 eV recovered, more than 30% of the discrepancy with experiment — comes entirely from allowing one additional variational parameter to respond to the chemical environment: the orbital contracts ( $\zeta$ increases from 1 to 1.24) precisely because two nuclei pull the electron density inward relative to an isolated hydrogen atom. Limiting-case check. As $R \to \infty$ (dissociation limit), the optimal $\zeta \to 1$ exactly, because the electron then sees only one nucleus and the system reduces to an isolated hydrogen atom — the variational method correctly "turns off" the contraction effect precisely when the physical justification for it (having two nearby nuclei) disappears.

Example 3 — Exam/competition level: variational bound for the linear (1D) harmonic oscillator using a non-Gaussian trial function. This example deliberately avoids the "obvious" Gaussian guess (which is already exact for the harmonic oscillator, so it would prove nothing about the method's reliability) and instead uses a Lorentzian-shaped trial function to show how far off a plausible but wrong-shaped guess can be, quantitatively. The Hamiltonian is $\hat H = -\dfrac{\hbar^2}{2m}\dfrac{d^2}{dx^2} + \dfrac{1}{2}k x^2$ , with exact ground energy $E_0 = \tfrac12\hbar\omega$ , $\omega = \sqrt{k/m}$ . Choose the trial function

\phi(x;\alpha) = \frac{1}{1+\alpha x^2},

which is smooth, even, peaked at the origin, and normalizable for $\alpha > 0$ — a physically reasonable-looking guess that is not Gaussian. Normalization. Using $\int_{-\infty}^{\infty} \dfrac{dx}{(1+\alpha x^2)^2} = \dfrac{\pi}{2\sqrt{\alpha^3}}$ (a standard Beta-function result), the norm is $N^2 = \dfrac{2\sqrt{\alpha^3}}{\pi}$ . Potential energy. Using $\displaystyle\int_{-\infty}^\infty \frac{x^2\,dx}{(1+\alpha x^2)^2} = \frac{\pi}{2\sqrt{\alpha}}$ ,

\langle V \rangle = \frac{1}{2}k \cdot N^2 \cdot \frac{\pi}{2\sqrt{\alpha}} = \frac{1}{2}k \cdot \frac{2\sqrt{\alpha^3}}{\pi}\cdot\frac{\pi}{2\sqrt{\alpha}} = \frac{k\alpha}{2}.

Kinetic energy. Differentiating twice, $\phi''(x) = \dfrac{6\alpha^2 x^2 - 2\alpha}{(1+\alpha x^2)^3}$ , and using the needed integrals $\displaystyle\int \frac{dx}{(1+\alpha x^2)^3} = \frac{3\pi}{8\sqrt{\alpha}}$ and $\displaystyle\int \frac{x^2\,dx}{(1+\alpha x^2)^3} = \frac{\pi}{8\sqrt{\alpha^3}}$ , one finds after integrating by parts (equivalently, using $\langle T\rangle = \frac{\hbar^2}{2m}\int (\phi')^2\,dx / \int \phi^2\,dx$ , which avoids fourth-order poles):

\phi'(x) = \frac{-2\alpha x}{(1+\alpha x^2)^2}, \qquad \int_{-\infty}^\infty (\phi')^2\,dx = 4\alpha^2 \int \frac{x^2\,dx}{(1+\alpha x^2)^4} = 4\alpha^2 \cdot \frac{\pi}{16\sqrt{\alpha^5}} = \frac{\pi\alpha^{3/2}}{4}.

Then

\langle T\rangle = \frac{\hbar^2}{2m}\cdot\frac{\pi\alpha^{3/2}/4}{\pi/(2\sqrt{\alpha})} = \frac{\hbar^2}{2m}\cdot\frac{\alpha^{3/2}}{4}\cdot\frac{2\sqrt{\alpha}}{1} = \frac{\hbar^2 \alpha^2}{4m}.

Assemble and minimize.

E(\alpha) = \frac{\hbar^2}{4m}\alpha^2 + \frac{k}{2}\alpha .

Differentiating,

\frac{dE}{d\alpha} = \frac{\hbar^2}{2m}\alpha + \frac{k}{2} = 0 .

This gives $\alpha^* = -mk/\hbar^2$ , which is negative — unphysical, since $\alpha$ must be positive for normalizability. The stationary point found by naive differentiation lies outside the admissible domain. This is the pedagogical payoff: the Lorentzian trial function, for this potential, has no interior minimum with $\alpha>0$ ; $E(\alpha) = \dfrac{\hbar^2\alpha^2}{4m} + \dfrac{k\alpha}{2}$ is monotonically increasing for all $\alpha > 0$ (both terms increase with $\alpha$ ), so the infimum over the admissible domain is approached as $\alpha \to 0^+$ , where $E(\alpha)\to 0$ — but $\alpha\to0$ makes $\phi$ spread out over all space and become non-normalizable in the limit, signaling that this shape is fundamentally unsuited to representing a localized ground state and the variational search must be redirected (e.g., restricting $\alpha$ to a physically sensible range, or, better, abandoning this trial family in favor of a Gaussian). Contrast with the correct shape. Repeating the identical procedure with a Gaussian $\phi(x;\alpha)=e^{-\alpha x^2}$ gives the textbook result $E(\alpha) = \dfrac{\hbar^2\alpha}{2m} + \dfrac{k}{8\alpha}$ , minimized at $\alpha^* = \dfrac{m\omega}{2\hbar}$ , reproducing $E(\alpha^*) = \tfrac12\hbar\omega$ exactly — because the Gaussian is the true ground state shape. The Lorentzian's total failure to produce an interior minimum, set directly against the Gaussian's exact success, is the sharpest possible illustration that the variational method's quality is bounded by the shape you choose, not by how many parameters you throw at it.

Limiting cases, scope & common traps

Situation	What happens	Why
Trial family contains the exact eigenfunction	$E(\alpha^*) = E_0$ exactly	The variational bound becomes an equality once the true wavefunction is reachable within the family (Example: scaled hydrogen 1s, Derivation 2)
Trial shape is qualitatively wrong (wrong asymptotic decay, wrong nodal count)	$E(\alpha^)$ can be badly above $E_0$ , or minimization may fail to converge to a physical $\alpha^$	The parameter can only rescale a fixed shape; it cannot fix a fundamentally wrong functional form (Example 3, Lorentzian case)
More parameters added to a reasonable shape	$E(\alpha^*)$ monotonically decreases (never increases) as parameters are added	Each additional degree of freedom enlarges the trial family, so the minimum over a larger set is at least as low
Excited states, naive treatment	Direct minimization of $E(\alpha)$ for an excited-state-shaped $\phi$ is not guaranteed to bound the corresponding excited energy unless orthogonality to all lower states is enforced	The basic variational theorem is proved only for the true ground state; a separate orthogonality-constrained theorem is needed for excited states
$\zeta \to Z$ limit in the helium screening example	Recovers unscreened, poorer bound $E = -Z^2 + (5/8)Z$	Removing the freedom to screen forces the electrons to see the bare nuclear charge, worsening the bound relative to $\zeta^*=Z-5/16$
Very large parameter count (many-parameter, many-Gaussian expansions)	Approaches the exact energy asymptotically but computational cost grows quickly	This is the practical basis of modern basis-set quantum chemistry, though full convergence is a separate topic from the two-parameter examples here

Common trapMisconception: "a lower variational energy always means a better wavefunction"

The variational theorem guarantees only that $E(\alpha^*)$ bounds $E_0$ from above; it says nothing about how close the wavefunction $\phi(\alpha^*)$ is to the true eigenfunction. Energy is a very forgiving functional — it is stationary to first order in any error of the wavefunction, so a trial function can be qualitatively wrong in some regions (e.g., missing angular correlation between electrons) yet still deliver a numerically good energy, especially if the errors are concentrated where $\hat H$ weights them lightly. Properties other than the energy (dipole moments, densities near the nucleus, polarizabilities) are far more sensitive to wavefunction quality and can be substantially wrong even from a trial function that reproduces $E_0$ to three significant figures.

Common trapMisconception: "you can just minimize over any parameter appearing in the Hamiltonian, like Z itself"

In the helium derivation, $\zeta$ is a parameter of the trial function, not of the physical Hamiltonian — the nuclear charge $Z=2$ in $\hat H$ stays fixed throughout. It is a common error to instead try to "optimize" the actual nuclear charge in the exact Hamiltonian, which is meaningless (nuclear charge is a physical constant of the atom, not a variational knob). The trick is precisely that $\zeta$ enters only through the shape of $\phi$ ; the true $Z$ still appears explicitly in the $-Z/r_i$ terms of $\hat H$ when evaluating $\langle\phi|\hat H|\phi\rangle$ , which is exactly why the derivation carefully separates $\zeta$ (trial) from $Z$ (physical) in Step 2 above.

The chemist behind it

The variational method as used in quantum chemistry traces to the mathematical variational principle developed by Lord Rayleigh and Walter Ritz for vibration problems in classical mechanics, later carried directly into quantum theory once Schrödinger's equation appeared. The chemist most associated with turning this into a working tool for electronic structure is Egil Hylleraas, whose 1929 calculation of the helium ground-state energy — using trial functions with explicit electron–electron distance dependence, going well beyond the simple screened-orbital product used above — achieved agreement with experiment to a fraction of a percent, a stunning result for a three-body Coulomb problem with no closed-form solution. Hylleraas showed that adding physically motivated correlation terms (functions of $r_{12}$ directly, not just $r_1$ and $r_2$ separately) to the trial function captures the instantaneous avoidance electrons exhibit around each other, an effect the simple product-of-orbitals trial function used in the derivation above cannot describe at all.

Origins · 1 minEgil Andersen Hylleraas and how this was discovered

Egil Andersen Hylleraas

1898–1965 · Early Norwegian theoretical physicist and pioneer of quantum-mechanical few-body calculations

Hylleraas worked at the University of Oslo and briefly with Max Born in Göttingen during the crucial years just after Schrödinger's 1926 papers, when almost no one had demonstrated that the new quantum mechanics could handle a genuinely interacting few-electron system quantitatively. His treatment of helium, published in 1929, remains a foundational demonstration that variational trial functions — built from physical insight about electron correlation rather than brute-force numerical solution — could rival experimental spectroscopic precision. The coordinates he introduced for two-electron atoms, built from $r_1$ , $r_2$ , and $r_{12}$ together, are still called Hylleraas coordinates and remain the gold-standard approach for benchmark two- and three-electron atomic calculations today.

History

- 1873 — Lord Rayleigh develops variational techniques for vibrational eigenvalue problems in classical acoustics.

- 1908 — Walter Ritz formalizes the variational method for eigenvalue problems, giving the "Rayleigh–Ritz" method its name.

- 1926 — Schrödinger's wave equation opens the possibility of applying variational reasoning to quantum systems.

- 1929 — Hylleraas applies explicitly correlated variational trial functions to helium, matching experimental spectroscopic energies.

- 1930s–present — Variational trial functions with screening/scaling parameters become standard teaching tools and the conceptual seed of Hartree–Fock and modern basis-set electronic structure methods.

Why does minimizing the energy over a parameter give an upper bound rather than the exact answer?

The variational theorem proves that for any normalized trial function, $\langle\phi|\hat H|\phi\rangle \ge E_0$ , with equality only if $\phi$ happens to be the true ground-state eigenfunction. Minimizing over a parameter $\alpha$ searches only within the restricted family of functions your particular functional form can produce; unless that family is flexible enough to contain the true eigenfunction exactly, the minimum you find will still sit above $E_0$ . The size of the gap depends entirely on how well the chosen shape matches reality — a poorly chosen shape can leave a large residual error even after perfect optimization, while a well-chosen shape (like the scaled hydrogenic exponential for the hydrogen atom itself) can close the gap to zero.

How do you know whether a stationary point found by

dE/d\alpha = 0

is a minimum and not a maximum or saddle?

You apply the ordinary second-derivative test from single-variable calculus: compute $d^2E/d\alpha^2$ at the stationary point and confirm it is positive. In the helium example, $d^2E/d\zeta^2 = 2 > 0$ everywhere, so the single stationary point found is unambiguously a minimum. In multi-parameter problems the analogous test uses the Hessian matrix of second partial derivatives, requiring it to be positive definite. Physically, because the variational energy is bounded below by $E_0$ for the entire trial family, an interior stationary point of a smooth, well-behaved $E(\alpha)$ landscape is very often a minimum, but the algebraic check should always be performed rather than assumed.

Can the variational method be used directly for excited states?

Yes, but only with an added condition: the basic theorem $\langle\phi|\hat H|\phi\rangle \ge E_0$ applies strictly to the ground state. To bound an excited-state energy, the trial function must be constrained to be orthogonal to the exact ground state (and to all lower excited states), giving $\langle\phi|\hat H|\phi\rangle \ge E_1$ for the first excited state, and so on. In practice this orthogonality is often enforced automatically by symmetry — for example, a trial function with odd parity cannot overlap with an even-parity ground state, so minimizing its energy directly bounds the lowest odd-parity excited state without any extra machinery.

Why did the Lorentzian trial function fail so badly for the harmonic oscillator, when it looked like a reasonable smooth, peaked, normalizable guess?

The failure traces to the asymptotic decay of the Lorentzian, which falls off as $1/x^2$ at large $|x|$ , far too slowly compared to the true Gaussian ground state's $e^{-m\omega x^2/2\hbar}$ decay. Because the harmonic potential grows as $x^2$ , any trial function with heavy power-law tails samples arbitrarily large potential energy from its slowly decaying wings, and the competition between kinetic energy (favoring spreading out, lowering $\langle T\rangle$ ) and potential energy (favoring localization) never balances at a finite width for this particular shape. This example demonstrates that qualitative features like tail behavior, not just the number of adjustable parameters, determine whether a trial family can succeed at all.

Is there a general rule for choosing a good trial function before doing any calculation?

The most reliable strategy is to encode as much correct physics as possible into the shape: correct nodal structure (matching the number of nodes expected for the state you are approximating), correct short-range behavior (a cusp at coalescence points, as required by the Kato cusp condition for Coulombic systems), and correct long-range decay (exponential falloff at a rate set by the ionization energy, for bound atomic and molecular states). Parameters are then reserved for quantities you are genuinely uncertain about, such as effective nuclear charges, orbital exponents, or mixing coefficients between basis functions, which is exactly the strategy used for the helium screening constant and the LCAO exponent in the worked examples above.

How does this single-parameter variational method relate to modern electronic structure calculations like Hartree-Fock or DFT?

Hartree–Fock theory is a direct multi-parameter generalization: instead of one scalar $\alpha$ , the "parameters" are the coefficients of atomic orbitals inside molecular orbitals, and the trial function is an antisymmetrized product (Slater determinant) built from those orbitals. Minimizing the energy with respect to all orbital coefficients simultaneously, subject to orthonormality constraints, produces the Hartree–Fock equations — the direct multi-parameter descendant of the single- $\alpha$ calculus practiced here. Optimizing Variational Parameters develops that generalization explicitly, while Linear Variation and Secular Equations handles the special (and extremely common) case where the parameters enter the trial function linearly, turning the minimization into a matrix eigenvalue problem rather than ordinary calculus.

Why does the helium calculation split the Hamiltonian into hydrogen-like pieces plus correction terms rather than just brute-force integrating the full Hamiltonian against the trial function?

Splitting the Hamiltonian is a technique for exploiting integrals you already know rather than a change to the underlying physics. Because $\psi_{1s}(\mathbf{r};\zeta)$ is the exact eigenfunction of a hydrogen-like Hamiltonian with nuclear charge $\zeta$ , expectation values of the bracketed kinetic-plus-attraction operator are known immediately, without any explicit integration, as $-\zeta^2/2$ per electron. Only the leftover screening term $(\zeta-Z)/r_i$ and the genuinely new electron-repulsion integral $1/r_{12}$ require separate evaluation, and both happen to be standard integrals over hydrogenic 1s functions available in closed form. This bookkeeping trick is standard throughout perturbative and variational treatments of many-electron atoms.