3. The variational method for the Schrödinger equation

Quantum systems are governed by the Schrödinger equation. In particular, the solutions to the stationary form of this equation determine many physical properties of the system at hand. The stationary Schrödinger equation can be solved analytically in a very restricted number of cases -- examples include the free particle, the harmonic oscillator and the hydrogen atom. In most cases we must resort to computers to determine the solutions to this equation. It is of course possible to integrate the Schrödinger equation using discretisation methods but in realistic electronic structure calculations for example we would need a huge number of grid points, leading to important computer time and memory requirements. On the other hand, the variational method enables us to solve the Schrödinger equation much more efficiently in many cases. In this chapter we introduce the variational method for solving the Schrödinger equation.

3.1. Variational calculus

In the variational method, the possible solutions are restricted to a subspace of the Hilbert space, and in this subspace we seek the best possible solution (below we shall define what is to be understood by the 'best' solution). To see how this works, we first show that the stationary Schrödinger equation can be derived by a stationarity condition of the functional $$E[\psi] = \frac{\int dX \psi^*(X) H \psi(X)}{\int dX \psi^*(X) \psi(X)} = \frac{\left\langle\psi |H| \psi \right\rangle}{\left\langle\psi |\psi\right\rangle},$$ which is recognised as the expectation value of the energy for a stationary state $\psi$ (in order to keep the analysis general we are not specific about the form of the generalised coordinate $X$ -- it may include the space and spin coordinates of a collection of particles). The stationary states of this energy functional are defined by postulating that if such a state is changed by an arbitrary $\delta \psi,$ the corresponding change in $E$ vanishes to first order. Formally, this means that $$\left. \frac{d}{d\lambda}\left( E\left[ \psi + \lambda \delta \psi \right] - E[\psi] \right) \right|_{\lambda =0} \equiv 0$$ for all normalised vectors $\delta \psi$. Defining $$\begin{split} P & = \left\langle\psi | H | \psi \right\rangle\phantom{xxx} \text{and} \\ Q & = \left\langle\psi | \psi \right\rangle, \end{split}$$ we can write the change $\delta E$ in the energy to first order in $\delta \psi$ as $$\begin{aligned} \delta E &= \frac{\left\langle\psi+\delta \psi| H |\psi+\delta\psi\right\rangle}{ \left\langle\psi+\delta \psi|\psi+\delta\psi\right\rangle} - \frac{\left\langle\psi |H| \psi \right\rangle}{\left\langle\psi |\psi\right\rangle} \nonumber \\ & \approx \frac{\left\langle\delta \psi | H | \psi \right\rangle-\frac{P}{Q}\left\langle\delta \psi | \psi \right\rangle}{Q} + \frac{\left\langle\psi | H | \delta \psi \right\rangle-\frac{P}{Q}\left\langle\psi | \delta \psi \right\rangle}{Q}. \end{aligned}$$ As this should vanish for an arbitrary but small change in $\psi$, we find, using $E=P/Q$: $$H\psi = E\psi , \tag{3.5}\label{Schrod}$$ together with the Hermitian conjugate of this equation, which is equivalent.

In variational calculus, stationary states of the energy functional are found within a subspace of the Hilbert space. An important example is linear variational calculus, in which the subspace is spanned by a finite set of basis vectors $\left|\chi_p\right\rangle$, $p=1,\ldots,N$, that we take to be orthonormal at first, that is, $$\left\langle\chi_p | \chi_q \right\rangle= \delta_{pq},$$ where $\delta_{pq}$ is the Kronecker delta-function which is 0 unless $p=q$, in which case it is 1.

For a state $$\left|\psi\right\rangle= \sum_{p=1}^N C_p \left|\chi_p\right\rangle,$$ the energy functional is given by $$E = \frac{\sum_{p,q=1}^N C_p^* C_q H_{pq}}{ \sum_{p,q=1}^N C_p^* C_q \delta_{pq}} \tag{3.8}\label{lin_var},$$ with $$H_{pq} = \left\langle\chi_p | H | \chi_q \right\rangle.$$ The stationary states follow from the condition that the derivative of this functional with respect to the $C_p$ vanishes, which leads to $$\sum_{q=1}^N \left( H_{pq} - E \delta_{pq} \right) C_q = 0 \phantom{xx} {\text{for}} \phantom{xx} p=1,\ldots,N. \tag{3.10}\label{lin_vareq}$$ Equation \eqref{lin_vareq} is an eigenvalue problem which can be written in matrix notation: $${\bf H C } = E {\bf C}. \tag{3.11}\label{eigeq}$$ This is the Schrödinger equation formulated for a finite, orthonormal basis.

Linear parametrisations are often used because the resulting method is simple, allowing for numerical matrix diagonalisation techniques to be used. The lowest eigenvalue of \eqref{eigeq} is always higher than or equal to the ground state energy of Eq. \eqref{Schrod}, as the ground state is the minimal value assumed by the energy functional over the full Hilbert space. If we restrict ourselves to a part of this space, then the minimum value of the energy functional must always be higher than or equal to the ground state of the full Hilbert space. Adding more basis functions to our set, the subspace becomes larger, and consequently the minimum of the energy functional will decrease (or stay the same). For the specific case of linear variational calculus, this result can be generalised to stationary states at higher energies: the higher eigenvalues are always higher than the equivalent solution to the full problem, but approximate the latter better with increasing basis set size. The formal statement of this is the Hylleraas-Undin-MacDonald theorem (see for example Springer Handbook of Atomic, Molecular, and Optical Physics, Volume 1, Gordon Drake (ed.), 2006). The behaviour of the spectrum found by solving \eqref{eigeq} with increasing basis size is depicted in Figure below.

The behaviour of the spectrum of Eq. \eqref{eigeq} with increasing basis set size in linear variational calculus. The upper index is the number of states in the basis set, and the lower index labels the spectral levels.

We now describe how to proceed when the basis consists of non-orthonormal basis functions, as is often the case in practical calculations. In that case, we must reformulate \eqref{eigeq}, taking care of the fact that the overlap matrix $\bf{S}$, whose elements $S_{pq}$ are given by $$S_{pq} = \left\langle\chi_p | \chi_q\right\rangle\label{defS}$$ is not the unit matrix. This means that in Eq. \eqref{lin_var} the matrix elements $\delta_{pq}$ of the unit matrix, occurring in the denominator, have to be replaced by $S_{pq}$, and we obtain $${\bf H C} = E {\bf S C}. \tag{3.13}\label{geneig}$$ This looks like an ordinary eigenvalue equation, the only difference being the matrix $\bf S$ in the right hand side. It is called a generalised eigenvalue equation and there exist computer programs for solving such a problem.