The Hellmann-Feynman Theorem

Solving a quantum system means finding its energy $E$ , but an enormous amount of physics lives instead in how that energy responds when we turn a knob — the field strength, the nuclear charge, a coupling constant, or the position of an atom in a molecule. The Hellmann–Feynman theorem is the exact statement that this response is astonishingly simple: to compute $dE/d\lambda$ you do not need to know how the wavefunction shifts as $\lambda$ changes. You only need the expectation value of $\partial H/\partial\lambda$ in the unchanged eigenstate. The reason is a small miracle of first-order perturbation bookkeeping. When $\lambda$ moves, the energy changes for two reasons: the Hamiltonian itself changes, and the state $|\psi_\lambda\rangle$ rearranges. The second effect looks like it should matter — but because $|\psi_\lambda\rangle$ is a stationary eigenstate, the pieces that track its rearrangement cancel exactly against each other, leaving only the direct term $\langle\psi_\lambda|\partial H/\partial\lambda|\psi_\lambda\rangle$ . The wavefunction's motion is a first-order effect on a quantity that is already stationary with respect to it, so it drops out. This one identity is a workhorse. It turns hard expectation values like $\langle 1/r\rangle$ and $\langle 1/r^2\rangle$ for hydrogen into one-line parameter derivatives; it delivers the virial theorem for the oscillator without an integral; and in its molecular form (Feynman's 1939 "Forces in Molecules") it proves that the force on a nucleus is just the classical electrostatic pull of the surrounding electron cloud — the conceptual foundation of every geometry optimization in modern quantum chemistry. The single non-negotiable condition, stressed throughout, is that $|\psi_\lambda\rangle$ must be an exact, normalized eigenstate: applied to a crude approximation the theorem quietly fails, and knowing why it fails is as important as knowing that it holds.

Key takeaways

For an exact normalized eigenstate $|\psi_\lambda\rangle$ of $H(\lambda)$ , the Hellmann–Feynman theorem states $\dfrac{dE}{d\lambda}=\left\langle\psi_\lambda\left|\dfrac{\partial H}{\partial\lambda}\right|\psi_\lambda\right\rangle$ .
The terms tracking how $|\psi_\lambda\rangle$ rearranges cancel because the state is stationary and normalized — you never need the derivative of the wavefunction.
It converts awkward expectation values into parameter derivatives: $\langle x^2\rangle$ for the oscillator from $\partial_\omega E$ , and $\langle 1/r\rangle=Z/(n^2a_0)$ and $\langle 1/r^2\rangle=Z^2/[n^3(\ell+\tfrac12)a_0^2]$ for hydrogen from $\partial_Z E$ and $\partial_\ell E$ .
Feynman's electrostatic form ( $\lambda=$ nuclear position) shows the force on a nucleus is the classical Coulomb force of the electron density — the basis of chemical-force calculations.
It fails for approximate or unnormalized states, holding for trial functions only at a variational stationary point.
Key idea: because an eigenstate is stationary, the energy's response to any parameter is just the average of the Hamiltonian's own response — the wavefunction's rearrangement contributes nothing to first order.

The core insight

For an exact normalized eigenstate, $dE/d\lambda=\langle\partial H/\partial\lambda\rangle$ : the wavefunction's change with $\lambda$ contributes zero because the state is stationary and normalized, so the energy derivative is the pure "explicit" derivative of $H$ .

Where this lives today

The Hellmann–Feynman theorem is the silent engine of computational materials science and quantum chemistry. Every geometry optimization — finding the shape of a molecule or the lattice constant of a crystal — computes interatomic forces as $-\partial E/\partial R_I=-\langle\partial H/\partial R_I\rangle$ , the expectation value of the (analytically known) derivative of the electron–nucleus attraction, rather than by expensive finite differences of the total energy. In density-functional theory the same identity, applied to the Kohn–Sham system, gives the forces that drive ab initio molecular dynamics (the Car–Parrinello and Born–Oppenheimer schemes). The theorem also underlies the calculation of a vast range of response properties: dipole moments as derivatives with respect to an applied field, polarizabilities as second derivatives, and pressures in solids as volume derivatives. Its one caveat — the failure for incomplete basis sets, cured by adding the "Pulay force" corrections — is a daily concern in production codes, and the reason well-converged basis sets matter. Conceptually, it remains the cleanest statement of a deep fact: for a stationary state, physics responds to a perturbation only through the perturbation's own explicit form, never through the reshuffling of the state.

The physical picture

Imagine the energy $E(\lambda)$ as the height of a landscape as you slide a control $\lambda$ . Naively, $E$ changes both because the rules ( $H$ ) change and because the system's response (the state $|\psi_\lambda\rangle$ ) shifts to accommodate the new rules. The Hellmann–Feynman theorem says the second contribution vanishes for a genuine eigenstate. The mental model is a ball resting at the bottom of a valley: if you tilt the valley slightly, the ball's height changes almost entirely because the ground moved, not because the ball rolled — the ball was already at a stationary point, so its own small displacement is a second-order effect. The eigenstate is that stationary ball; the energy is stationary with respect to variations of the state (the variational principle), so state-motion is invisible to first order.

Figure

Two contributions to dE/dλ, and why one cancels

The cancellation has a precise engine. Differentiate $E=\langle\psi|H|\psi\rangle$ and three terms appear. In the two that carry $\partial_\lambda\psi$ , the Hamiltonian acts on an eigenstate and becomes a number $E$ ; what remains is $E$ times the derivative of $\langle\psi|\psi\rangle$ . But $\langle\psi|\psi\rangle=1$ is constant, so its derivative is zero. The physics is that the state's rearrangement is orthogonal, at leading order, to the state itself, and $H$ cannot see it.

Figure

The stationary-ball analogy

The core equations

The Hamiltonian depends smoothly on a real parameter $\lambda$ , and for each $\lambda$ we have a normalized eigenstate. This is the entire setup.

H(\lambda)\,|\psi_\lambda\rangle = E(\lambda)\,|\psi_\lambda\rangle, \qquad \langle\psi_\lambda|\psi_\lambda\rangle = 1 \qquad (1)

where: $H(\lambda)$ is the Hamiltonian carrying an explicit dependence on the parameter $\lambda$ (a field strength, charge, coupling, or coordinate); $|\psi_\lambda\rangle$ is the exact normalized eigenstate; and $E(\lambda)$ is the corresponding eigenvalue. Both the state and the eigenvalue inherit a $\lambda$ -dependence from $H$ . The normalization is enforced at every $\lambda$ , which is the condition that does the work below. The theorem itself is a single line: the total derivative of the eigenvalue equals the expectation value of the partial derivative of the Hamiltonian.

\frac{dE}{d\lambda} = \left\langle\psi_\lambda\left|\frac{\partial H}{\partial\lambda}\right|\psi_\lambda\right\rangle \qquad (2)

where: $dE/d\lambda$ is the total rate of change of the energy (including all implicit dependence through the state); $\partial H/\partial\lambda$ is the explicit (partial) derivative of the Hamiltonian operator, holding the coordinates fixed; and the expectation value is taken in the eigenstate at that same $\lambda$ . The remarkable content is that the left side is a total derivative while the right side needs only a partial one — the difference between them, which would involve $\partial_\lambda\psi$ , is exactly zero. For applications it helps to state the two most-used specializations. With $\lambda=\omega$ for the harmonic oscillator, $\partial H/\partial\omega=m\omega x^2$ , so equation (2) reads

\frac{dE_n}{d\omega} = m\omega\,\langle x^2\rangle \;\Longrightarrow\; \langle x^2\rangle = \frac{1}{m\omega}\frac{dE_n}{d\omega} = \frac{(n+\tfrac12)\hbar}{m\omega} \qquad (3)

where: $E_n=(n+\tfrac12)\hbar\omega$ is the oscillator spectrum; $m$ the mass; $\omega$ the angular frequency; and $\langle x^2\rangle$ the position variance in the $n$ th state. The right-hand expression is obtained purely by differentiating the closed-form $E_n$ — no wavefunction integral is performed. This is also the virial theorem for the oscillator, since $\langle\tfrac12 m\omega^2 x^2\rangle=\tfrac12(n+\tfrac12)\hbar\omega=\tfrac12 E_n$ . With $\lambda=Z$ for a hydrogenic atom in atomic units ( $H=-\tfrac12\nabla^2-Z/r$ , so $\partial H/\partial Z=-1/r$ ),

E_n=-\frac{Z^2}{2n^2}\;\Longrightarrow\; \frac{dE_n}{dZ}=-\frac{Z}{n^2}=-\left\langle\frac1r\right\rangle \;\Longrightarrow\; \left\langle\frac1r\right\rangle=\frac{Z}{n^2 a_0} \qquad (4)

where: $Z$ is the nuclear charge treated as a continuous parameter; $n$ the principal quantum number; $a_0$ the Bohr radius (equal to $1$ in atomic units, restored explicitly on the right); and $\langle 1/r\rangle$ the mean inverse radius. The last equality is the Hellmann–Feynman result: an expectation value that ordinarily requires a Laguerre-polynomial integral is read off from the energy in one step.

Derivation ADeriving the theorem from stationarity and normalization

The goal is to prove equation (2) — that the total energy derivative reduces to the expectation value of the explicit Hamiltonian derivative — with no assumption beyond that $|\psi_\lambda\rangle$ is an exact normalized eigenstate for every $\lambda$ . The route is to write $E$ as an expectation value, differentiate it by the product rule into three terms, and then use the eigenvalue equation to collapse the two "state-motion" terms into a multiple of the derivative of the norm, which vanishes. The one thing to watch is that we are allowed to let $H$ act to the left on the bra as $\langle\psi_\lambda|H=E\langle\psi_\lambda|$ , which requires $H$ to be Hermitian and $|\psi_\lambda\rangle$ to be a genuine eigenstate; this is precisely where an approximate state would break the argument.

Step 1. Write the energy as an expectation value using the normalization: $E(\lambda)=\langle\psi_\lambda|H(\lambda)|\psi_\lambda\rangle$ , valid because $\langle\psi_\lambda|\psi_\lambda\rangle=1$ .

Step 2. Differentiate with respect to $\lambda$ by the product rule, treating bra, operator, and ket as three factors: $\dfrac{dE}{d\lambda}=\langle\partial_\lambda\psi|H|\psi\rangle+\langle\psi|\,\partial_\lambda H\,|\psi\rangle+\langle\psi|H|\partial_\lambda\psi\rangle$ , where $\partial_\lambda H$ is the explicit derivative and $|\partial_\lambda\psi\rangle$ the state's response.

Step 3. Use the eigenvalue equation on the first and third terms. Acting to the right, $H|\psi\rangle=E|\psi\rangle$ ; acting to the left (Hermiticity), $\langle\psi|H=E\langle\psi|$ . So the two state-motion terms become $E\langle\partial_\lambda\psi|\psi\rangle+E\langle\psi|\partial_\lambda\psi\rangle=E\big(\langle\partial_\lambda\psi|\psi\rangle+\langle\psi|\partial_\lambda\psi\rangle\big)$ . Why this is the crux: the operator has been replaced by the number $E$ , so what is left measures only the overlap of the state with its own derivative.

Step 4. Recognize the bracket as the derivative of the norm: $\langle\partial_\lambda\psi|\psi\rangle+\langle\psi|\partial_\lambda\psi\rangle=\dfrac{d}{d\lambda}\langle\psi|\psi\rangle=\dfrac{d}{d\lambda}(1)=0$ .

Step 5. The two state-motion terms vanish, leaving $\dfrac{dE}{d\lambda}=\langle\psi_\lambda|\partial_\lambda H|\psi_\lambda\rangle$ , which is equation (2). The wavefunction derivative never had to be computed. Verified: the companion verifier confirms symbolically that $E\,\dfrac{d}{d\lambda}\langle\psi|\psi\rangle=0$ under $\langle\psi|\psi\rangle=1$ , isolating the cancellation, and then checks the full identity on the oscillator and hydrogen (below).

Go deeper · Derivation B · 1 minExtracting hydrogen expectation values by the parameter trick

Derivation BExtracting hydrogen expectation values by the parameter trick

Here the aim is to show that the theorem turns two genuinely awkward hydrogen integrals, $\langle 1/r\rangle$ and $\langle 1/r^2\rangle$ , into one-line derivatives of the energy — the first by differentiating with respect to the nuclear charge $Z$ , the second by the more subtle trick of differentiating with respect to the angular-momentum quantum number $\ell$ treated as a continuous parameter. The route for $\langle 1/r\rangle$ is direct; for $\langle 1/r^2\rangle$ we must first notice that $\ell$ enters the radial Hamiltonian only through the centrifugal term $\ell(\ell+1)/2r^2$ , and that the energy depends on $\ell$ only through the combination $n=n_r+\ell+1$ . The subtlety to watch is that $\ell$ is normally an integer, so we are extending the exact spectrum to non-integer $\ell$ — legitimate because the radial equation is analytic in $\ell$ and its bound-state solutions vary smoothly.

Step 1. For $\langle 1/r\rangle$ , take $\lambda=Z$ in the atomic-unit Hamiltonian $H=-\tfrac12\nabla^2-Z/r$ . Then $\partial H/\partial Z=-1/r$ , so equation (2) gives $dE_n/dZ=-\langle 1/r\rangle$ .

Step 2. Differentiate the closed-form energy $E_n=-Z^2/2n^2$ : $dE_n/dZ=-Z/n^2$ . Equate to $-\langle 1/r\rangle$ to get $\langle 1/r\rangle=Z/n^2$ (in units of $1/a_0$ ), i.e. $\langle 1/r\rangle=Z/(n^2 a_0)$ .

Step 3. For $\langle 1/r^2\rangle$ , write the radial Hamiltonian with the centrifugal barrier explicit: $H_\ell=-\tfrac{1}{2r^2}\partial_r(r^2\partial_r)+\dfrac{\ell(\ell+1)}{2r^2}-\dfrac{Z}{r}$ . Only the barrier carries $\ell$ , so $\partial H_\ell/\partial\ell=\dfrac{2\ell+1}{2r^2}$ .

Step 4. Apply the theorem treating $\ell$ as continuous: $\dfrac{dE}{d\ell}=\dfrac{2\ell+1}{2}\left\langle\dfrac1{r^2}\right\rangle$ . Since $E_n=-Z^2/2n^2$ with $n=n_r+\ell+1$ (so $dn/d\ell=1$ at fixed radial node number $n_r$ ), $\dfrac{dE}{d\ell}=\dfrac{dE}{dn}=\dfrac{Z^2}{n^3}$ .

Step 5. Solve: $\left\langle\dfrac1{r^2}\right\rangle=\dfrac{2}{2\ell+1}\cdot\dfrac{Z^2}{n^3}=\dfrac{Z^2}{n^3(\ell+\tfrac12)a_0^2}$ , the standard Kramers result — obtained without ever integrating a squared radial function. Verified: the verifier computes $\langle 1/r\rangle$ and $\langle 1/r^2\rangle$ directly from the normalized hydrogen radial functions for $1s,2s,2p,3d,3p$ and matches the Hellmann–Feynman formulas exactly ( $\langle 1/r\rangle_{1s}=1/a_0$ , $\langle 1/r^2\rangle_{2p}=1/[(3/2)\cdot 8]\,a_0^{-2}$ , etc.).

Limiting cases

Each limit is a consistency check or a boundary where the theorem's single hypothesis (exact normalized eigenstate) is stressed.

Case	Result	Meaning
Exact normalized eigenstate	$dE/d\lambda=\langle\partial_\lambda H\rangle$ holds exactly	The generic theorem; wavefunction response drops out.
Variational trial state at its optimum	Theorem holds for $\langle\partial_\lambda H\rangle$	At a stationary point of the variational energy, $\partial E/\partial(\text{params})=0$ removes the correction — this is why it is safe in optimized calculations.
Approximate / non-stationary state	Theorem fails	The state-motion terms no longer cancel; $dE/d\lambda\ne\langle\partial_\lambda H\rangle$ in general (see the Warning).
Degenerate level	Use the "good" basis	$\|\psi_\lambda\rangle$ must be the eigenstate that diagonalizes $\partial_\lambda H$ within the degenerate subspace, exactly as in degenerate perturbation theory.
$\lambda=$ nuclear coordinate (classical limit)	Feynman electrostatic force	The force on a nucleus is the classical Coulomb force of the electron density plus nucleus–nucleus repulsion — no kinetic-energy or exchange terms survive.

Common trapApplying the theorem to an approximate wavefunction

The Hellmann–Feynman theorem is not valid for an arbitrary approximate state. Its derivation used $H|\psi\rangle=E|\psi\rangle$ twice; if $|\psi\rangle$ is only an approximation, that step is false and the state-motion terms $E\,d\langle\psi|\psi\rangle/d\lambda$ do not fully account for the change — a residual term proportional to how badly $|\psi\rangle$ solves the Schrödinger equation survives. The one rescue is a variational trial state evaluated at its optimized parameters: there the energy is stationary against the variational parameters, so their $\lambda$ -induced drift contributes nothing and the theorem is restored (this is the "generalized Hellmann–Feynman" statement and the reason force calculations on optimized wavefunctions are trustworthy, while forces on unconverged or basis-set-incomplete wavefunctions suffer the notorious "Pulay forces"). Rule of thumb: use it freely on exact solutions and on fully optimized variational states; never on a hand-waved guess.

Worked examples

Example 1 — Mean inverse radius of hydrogen from the energy. Use the Hellmann–Feynman theorem to obtain $\langle 1/r\rangle$ for the hydrogen ground state, and evaluate it in physical units.

Step 1. Choose the parameter and identify $\partial H/\partial\lambda$ . Take $\lambda=Z$ with $H=-\tfrac12\nabla^2-Z/r$ (atomic units). The only $Z$ -dependence is in the potential, so $\partial H/\partial Z=-1/r$ .

Step 2. Differentiate the closed-form energy. With $E_n=-Z^2/2n^2$ , $dE_n/dZ=-Z/n^2$ . For the ground state $n=1$ : $dE_1/dZ=-Z$ .

Step 3. Apply equation (2): $dE_1/dZ=\langle-1/r\rangle=-\langle 1/r\rangle$ , hence $\langle 1/r\rangle=Z$ . Setting $Z=1$ for hydrogen, $\langle 1/r\rangle=1$ in atomic units, i.e. $\boxed{\langle 1/r\rangle=1/a_0}$ .

Step 4. Convert to physical units. With $a_0=0.529\ \text{Å}$ , $\langle 1/r\rangle=1/a_0=1.89\ \text{Å}^{-1}$ . Note this is not $1/\langle r\rangle$ : since $\langle r\rangle_{1s}=\tfrac32 a_0$ , $1/\langle r\rangle=0.67/a_0\ne\langle 1/r\rangle$ — the two differ because $1/r$ is convex, a small but classic trap. The direct Laguerre integral gives the same $1/a_0$ ; the takeaway is that the parameter derivative replaces the integral entirely.

Example 2 — Position spread of an oscillator state without integrating. Find $\langle x^2\rangle$ for the $n$ th harmonic-oscillator state and verify the virial split of the energy.

Step 1. Identify the parameter derivative. With $H=p^2/2m+\tfrac12 m\omega^2 x^2$ , $\partial H/\partial\omega=m\omega x^2$ .

Step 2. Differentiate the spectrum $E_n=(n+\tfrac12)\hbar\omega$ : $dE_n/d\omega=(n+\tfrac12)\hbar$ .

Step 3. Apply the theorem: $(n+\tfrac12)\hbar=m\omega\langle x^2\rangle$ , so $\langle x^2\rangle=\dfrac{(n+\tfrac12)\hbar}{m\omega}$ . For the ground state $n=0$ , $\langle x^2\rangle=\hbar/2m\omega$ — the minimum-uncertainty spread.

Step 4. Check the virial theorem. The mean potential energy is $\langle V\rangle=\tfrac12 m\omega^2\langle x^2\rangle=\tfrac12(n+\tfrac12)\hbar\omega=\tfrac12 E_n$ , so kinetic and potential energies each carry half the total — exactly the virial result for a quadratic potential, here delivered as a corollary. The Hermite-function integral reproduces $\langle x^2\rangle=(n+\tfrac12)\hbar/m\omega$ ; the point is that the theorem gets it from the spectrum alone.

Origins · 1 minrichard feynman and how this was discovered

richard feynman

1918–1988 · American theoretical physicist who, as an MIT undergraduate, proved the electrostatic form of the theorem

Richard Phillips Feynman derived the molecular-force version of this theorem in his 1939 MIT bachelor's thesis, published as "Forces in Molecules" (Phys. Rev. 56, 340). His result — that the force on a nucleus is the ordinary classical electrostatic attraction of the surrounding electron charge density, with no mysterious quantum force term — became the conceptual bedrock of computational chemistry. Feynman went on to co-invent quantum electrodynamics, for which he shared the 1965 Nobel Prize with Schwinger and Tomonaga, and to develop the path-integral formulation of quantum mechanics. The theorem's fuller pedigree is a lesson in credit: the identity had been written down independently by Güttinger, Pauli, and Hellmann before Feynman, and the modern name honors only two of at least four discoverers.

History

1932

— Paul Güttinger derives the energy-derivative identity in a paper on the Stark effect of rotating molecules

1933

— Wolfgang Pauli states the theorem in his Handbuch der Physik article on wave mechanics

1937

— Hans Hellmann presents it in his textbook "Einführung in die Quantenchemie," the first quantum-chemistry text

1939

— Richard Feynman, as an MIT undergraduate, proves the electrostatic (molecular-force) form in "Forces in Molecules" (Phys. Rev. 56, 340)

1969

— Pulay analyzes the "Pulay forces" that arise when the theorem is applied to incomplete-basis wavefunctions

Going further

References. Griffiths & Schroeter, Introduction to Quantum Mechanics, 3rd ed., §6 (perturbation theory, HF problems) · Feynman, "Forces in Molecules," Phys. Rev. 56, 340 (1939) · Hellmann, Einführung in die Quantenchemie (1937) · Sakurai & Napolitano, Modern Quantum Mechanics, 2nd ed., §5 · Cohen-Tannoudji, Diu & Laloë, Quantum Mechanics, complement on the theorem · Pulay, Mol. Phys. 17, 197 (1969).

Provenance. A verifier should reproduce: (1) the cancellation $E\,d\langle\psi|\psi\rangle/d\lambda=0$ under normalization; (2) the oscillator result $\langle x^2\rangle=(n+\tfrac12)\hbar/m\omega$ both from $\partial_\omega E_n$ and from Hermite eigenfunctions ( $n=0,1,2$ ); (3) $\langle 1/r\rangle=Z/(n^2 a_0)$ from $\partial_Z E_n$ and from the radial functions ( $1s,2s,2p$ ); (4) $\langle 1/r^2\rangle=Z^2/[n^3(\ell+\tfrac12)a_0^2]$ from the continuous- $\ell$ trick and from radial functions ( $2p,3d,3p$ ); (5) the failure of the identity for a deliberately frozen non-eigenstate. The historical dates (Güttinger 1932, Pauli 1933, Hellmann 1937, Feynman 1939) are as cited.

Builds on: Time-Independent Schrödinger Equation · Expectation Values · Quantum Harmonic Oscillator: Analytic Solution · Hydrogen Coulomb Radial Equation · The Variational Principle. Leads to: The Virial Theorem · Stationary-State Perturbation Theory · Molecular Forces and Geometry Optimization · Density-Functional Theory.

Why is the "partial" derivative of H the only thing that survives, when E clearly depends on λ through the state too?

Because the state's contribution cancels. When you differentiate $E=\langle\psi|H|\psi\rangle$ , you get three terms; the two carrying $\partial_\lambda\psi$ collapse, via $H|\psi\rangle=E|\psi\rangle$ , into $E\,d\langle\psi|\psi\rangle/d\lambda$ , which is zero because the norm is fixed at one. So although $E$ genuinely depends on $\lambda$ through the state, that dependence contributes nothing to the derivative at first order — the hallmark of a stationary quantity. The total derivative on the left and the partial (explicit) derivative on the right therefore agree.

Does the theorem require the parameter λ to be small, like in perturbation theory?

No. Unlike perturbation theory, the Hellmann–Feynman theorem is exact for any value of $\lambda$ , large or small, as long as $|\psi_\lambda\rangle$ is the exact eigenstate at that value of $\lambda$ . It is a statement about the local slope $dE/d\lambda$ of the exact energy curve, evaluated point by point. Perturbation theory approximates $E(\lambda)$ as a power series around $\lambda=0$ ; Hellmann–Feynman instead gives the exact derivative everywhere, and in fact the first-order perturbation energy $E^{(1)}=\langle\psi^{(0)}|V|\psi^{(0)}\rangle$ is just the theorem applied at $\lambda=0$ to $H=H_0+\lambda V$ .

What exactly goes wrong if I use an approximate wavefunction?

The step $H|\psi\rangle=E|\psi\rangle$ becomes false, so the two state-motion terms no longer combine into $E\,d\langle\psi|\psi\rangle/d\lambda$ . Concretely, differentiating $E_{\text{approx}}=\langle\phi|H|\phi\rangle/\langle\phi|\phi\rangle$ leaves a residual term proportional to $\langle\partial_\lambda\phi|(H-E)|\phi\rangle+\text{c.c.}$ , which vanishes only if $|\phi\rangle$ is a true eigenstate or a variational stationary point. For an unconverged calculation this residual is real and is the origin of "Pulay forces" in quantum chemistry: forces computed naively as $\langle\partial_\lambda H\rangle$ are wrong unless you add the correction or use a complete/optimized basis.

How is this related to the virial theorem?

The virial theorem falls out as a special case. Rescaling the coordinate $x\to\alpha x$ introduces a parameter $\alpha$ into $H$ ; demanding $dE/d\alpha=0$ at the physical value $\alpha=1$ (the eigenstate is stationary under this scaling too) gives a relation between $\langle T\rangle$ and $\langle r\cdot\nabla V\rangle$ — exactly the quantum virial theorem $2\langle T\rangle=\langle r\cdot\nabla V\rangle$ . For the oscillator this yields $\langle T\rangle=\langle V\rangle=\tfrac12 E$ , which we also derived directly from the $\omega$ -derivative above; for Coulomb it gives $\langle T\rangle=-\tfrac12\langle V\rangle=-E$ . The Hellmann–Feynman theorem with a scaling parameter is the cleanest proof of the virial theorem.

What is Feynman's "electrostatic theorem" and why is it surprising?

Take $\lambda$ to be a nuclear coordinate $\mathbf{R}_I$ . The only part of $H$ that depends on it is the electron–nucleus Coulomb attraction and the nucleus–nucleus repulsion, so $\partial H/\partial\mathbf{R}_I$ is a pure Coulomb-force operator. The theorem then says the force on nucleus $I$ is exactly the classical electrostatic force exerted by the electron charge density $\rho(\mathbf{r})=|\psi|^2$ and the other nuclei — with no kinetic-energy contribution, no exchange force, nothing quantum in the force expression itself. This is surprising because the energy is thoroughly quantum mechanical, yet the force is purely classical electrostatics acting on the quantum-mechanically-determined electron cloud. It is the reason molecular structure can be understood as electrostatic equilibrium of nuclei in an electron glue.

The theorem needs a normalized eigenstate — what happens in a degenerate level?

You must use the "good" states: within a degenerate eigenspace, the eigenstate to which the theorem applies is the particular linear combination that diagonalizes $\partial H/\partial\lambda$ inside that subspace. If you pick an arbitrary member of the degenerate space, $dE/d\lambda$ is ambiguous because the level splits and each branch has its own slope. This is exactly the same requirement as in degenerate perturbation theory, where one first diagonalizes the perturbation in the degenerate subspace before applying first-order formulas. Choose the branch, use its state, and the theorem gives that branch's slope.

Can I use it to get ⟨1/r²⟩ even though ℓ is an integer?

Yes, by analytic continuation in $\ell$ . The radial Schrödinger equation is analytic in the centrifugal parameter, and its bound-state energies vary smoothly if you pretend $\ell$ is a continuous real variable, with the physical spectrum recovered at integer $\ell$ . Since $\ell$ enters only through $\ell(\ell+1)/2r^2$ , we have $\partial H/\partial\ell=(2\ell+1)/2r^2$ , and $dE/d\ell$ (through $n=n_r+\ell+1$ ) gives $\langle 1/r^2\rangle=Z^2/[n^3(\ell+\tfrac12)a_0^2]$ . The trick is legitimate because we only evaluate the final formula at integer $\ell$ , where it matches the direct radial integral exactly — the verifier confirms this for several states.