3.3. Examples of non-linear variational calculus

Linear variational calculus leads to a (generalised) matrix eigenvalue problem to be solved. Variational methods form a much wider class, including trial functions which depend non-linearly on the variational parameters. Numerically this is quite complicated to solve, but several analytic non-linear variational calculations exist which give quite good results even if only one or two variational parameters are used. A nice example is the hydrogen atom, which we try to solve using a variational wave function ('trial function') $\left| \psi_{\text{T}}\right\rangle$ of the form $$\left\langle{\bf r}| \psi_{\text{T}} \right\rangle\propto e^{-r/a}.$$ You may note that this is the form of the exact ground state of the hydrogen atom. However, to illustrate the variational method, we first relax the value of the parameter $a$ and then vary it to minimise the expectation value of the energy. We should then find the exact ground state wave function and energy.

The Schrödinger equation was already given in atomic units in the previous section: $$\left[ -\frac{1}{2} {\bf \nabla}^2 - \frac{1}{r} \right] \psi({\bf r}) = E \psi ({\bf r}).$$

It is useful to first normalise the trial wave function: $$4\pi \int r^2 e^{-2r/a} dr = \pi a^3,$$ so that we have $$\left\langle{\bf r}| \psi_{\text{T}} \right\rangle= \frac{1}{\sqrt{\pi a^3}} e^{-r/a}.$$

Its is now easy to calculate the expectation value of the kinetic energy. Using the fact that for a function $\psi({\bf r})$ in 3D which only depends on $r$, $$\nabla^2 \psi({\bf r}) = \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{d}{dr} \psi(r)\right),$$ we have $$\nabla^2 \psi_{\text{T}} ({\bf r}) = \frac{1}{\sqrt{\pi a^3}} \left( \frac{1}{a^2} - \frac{2}{ar}\right)e^{-r/a},$$ we find, after some calculation that $$- \left\langle\psi_{\text{T}} \left| \frac{1}{2} \nabla^2 \right| \psi_{\text{T}} \right\rangle= \frac{1}{2\pi a^3} \int \left( -\frac{1}{a^2}+\frac{2}{ar} \right) e^{-2r/a} 4 \pi r^2 dr = \frac{1}{2 a^2}.$$ For the potential energy, we find $$-\frac{1}{\pi a^3} \int e^{-2r/a} 4 \pi r dr = - \frac{1}{a}.$$ Therefore, the expectation value of the energy for the trial wave function is given by $$E_{\text{T}} = \frac{1}{2 a^2} - \frac{1}{a}.$$ The minimum of this expression is found at $a=1$ and yields an energy of $E_{\text{T}} = -1/2$ in units of $27.212$ eV, which is the correct ground state energy of $-13.6058$ eV.

Now we turn to a more complicated problem: the helium atom, which (when the nucleus is considered not to move because of its large mass) is described by the Hamiltonian $$H = \frac{p_1^2}{2m} + \frac{p_2^2}{2m} - \frac{2e^2}{4\pi \epsilon_0 r_1} - \frac{2e^2}{4\pi \epsilon_0 r_2} + \frac{e^2}{4\pi \epsilon_0 |{\bf r}_1 - {\bf r}_2|}.$$ In atomic units, this becomes, with $r_{12} = |{\bf r}_1 - {\bf r}_2|$: $$H = -\frac{1}{2} \frac{d^2}{dr_1^2} -\frac{1}{2} \frac{d^2}{dr_2^2} - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{r_{12}}.$$ For the trial wave function we use the form $$\left\langle{\bf r}_1, {\bf r}_2 | \psi \right\rangle= \exp\left[ -2 (r_1+r_2)/a\right].$$ This function is chosen such that it yields the ground state of two noninteracting electrons moving in the field of the helium nucleus (the nuclear charge $Z=2$ leads to a scaling of 2 in the exponent). So the trial wave function is simply the successful trial wave function (in the sense that it contains the exact solution) for the independent-electron case. In particular, when an electron approaches the nucleus, its behaviour is properly described by this wave function as the electron-nucleus interaction largely dominates the electron-electron interaction in that case.

We have taken the wave function to be symmetric in the coordinates ${\bf r}_1$ and ${\bf r}_2$. The two electrons should however form an antisymmetric wave function as they are fermions. The antisymmetry is taken care of by the spin wave function $$\frac{1}{\sqrt{2}}\left( \left| \uparrow \downarrow \right\rangle- \left| \downarrow \uparrow \right\rangle\right).$$ Later we shall go much deeper into the structure of many-body wave functions -- here we just mention that, as the Hamiltonian does not contain any spin dependence, we can forget about the spin part of the wave function and can safely assume that the ground state wave function is symmetric in ${\bf r}_1$ and ${\bf r}_2$.

We first must normalise this solution. This can be done for ${\bf r}_1$ and ${\bf r}_2$ independently and we obtain $$4 \pi \int r^2 dr e^{-4r/a} = \frac{\pi a^3}{8},$$ so that the normalised solution reads $$\label{Eq:PsiTHe} \left\langle{\bf r}_1 {\bf r}_2 | \psi \right\rangle= \frac{8}{\pi a^3} \exp\left[ -2 (r_1+r_2)/a\right].$$ In order to find the expectation value of the Hamiltonian for this wave function, it is convenient to write it in the form $$H = H_1 + H_2 + \frac{1}{|{\bf r}_1 - {\bf r}_2|},$$ where $$H_i = \frac{p_i^2}{2m} - \frac{2}{r_i}, \phantom{xxx} i=1,2.$$ This is the Hamiltonian for an electron moving in the helium potential.

We calculate the kinetic energy following the method used for the hydrogen atom. The result is (for the two electrons together) $$E_{\text{K}} = \frac{4}{a^2}.$$ A quick way to arrive at this result is by taking the kinetic energy $1/(2a^2)$ for the hydrogen atom, replacing $a\rightarrow a/2$ and multiplying by 2 because we now have two electrons. The potential energy due to the attraction between the nucleus and the electrons is found to be $$E_{\text{n-e}} = \frac{-8}{a},$$ as can be verified by taking the hydrogen result $1/a$, replacing $a\rightarrow a/2$ and multiplying by 2 because the nuclear charge is twice as large as for the hydrogen atom and again by 2 because we now have two electrons.

Now we must add to this the contribution from the electron repulsion: $$\int \frac{64 e^{-4(r_1 + r_2)/a}}{\pi^2 a^6} \frac{1}{|{\bf r}_1 - {\bf r}_2|} \; d^3r_1\; d^3r_2.$$ If we fix ${\bf r}_1$, we can evaluate the integral over ${\bf r}_2$. Choosing the $z$-axis to be the direction of ${\bf r}_1$, we have (without the prefactor): $$\int e^{-4(r_1 + r_2)/a} \frac{1}{\sqrt{r_1^2 + r_2^2 - 2 r_1 r_2 \cos\theta}} 2\pi \sin \theta d\theta r_2^2 dr_2.$$ We first evaluate the integral over $\theta$. Choosing $\cos\theta = u$, this is of the form $\int du/\sqrt{ p - qu}$ and we are left with $$\begin{gathered} 2\pi \int e^{-4(r_1 + r_2)/a} \frac{1}{r_1 r_2} \left( \sqrt{r_1^2 + r_2^2 + 2 r_1 r_2} - \sqrt{r_1^2 + r_2^2 - 2 r_1 r_2} \right) r_2^2 dr_2 = \\ 2\pi \int e^{-4(r_1 + r_2)/a} \left[ \frac{1}{r_1 r_2} \left( r_1 + r_2 - |r_1 - r_2| \right) \right] r_2^2 dr_2 \end{gathered}$$ The term in square brackets equals $2/\max(r_1, r_2)$ where the function $\max(x, y)$ returns the largest of the two numbers $x$ and $y$. Therefore, we need to split the integral over $r_2$ into a part running from $0$ to $r_1$ and a part running from $r_1$ to $\infty$: $$\begin{aligned} & 4\pi e^{-4r_1/a} \left[ \frac{1}{r_1} \int_0^{r_1} e^{ -4 r_2/a} r_2^2 \; dr_2 + \int_{r_1}^\infty e^{ -4 r_2/a} r_2 \; dr_2 \right] = \\ & \frac{\pi a^3}{8r_1}e^{-4r_1/a} \left[ 1 - e^{-4r_1/a} - 2\frac{r_1}{a} e^{-4r_1/a}\right]. \end{aligned}$$

Now it remains to multiply this by the normalisation factor $64/(\pi^2 a^6)$ and by $\exp(-4r_1/a)$, and then integrate over $4\pi r_1^2dr_1$ (there is no dependence on the two angular variables for ${\bf r}_1$). All the integrals are straightforward and the final result is $$\left\langle\psi_{\text{T}} \left| \frac{1}{|{\bf r}_1 - {\bf r}_2|} \right|\psi_{\text{T}}\right\rangle= \frac{5}{4a}$$ in units of Hartree = $27.212$ eV. So the total energy is given by $$E_T = \frac{4}{a^2} - \frac{-8}{a} + \frac{5}{4a} = \frac{4}{a^2} - \frac{27}{4a}.$$ Taking $a=1$, i.e. assuming that both electrons are in the ground state of the atom with the electron-electron interaction switched off, yields an energy of $$4 - 27/4 = - 11/4 {\text{ Hartree}} = 74.833{\text{ eV}}.$$ The experimental value is $-79$ eV, so although our result is not extremely bad, it is not impressively accurate either.

Now let us relax $a$ and find the minimum of the trial energy. This is found at $a=32/27=1.1815$ Bohr radii. The energy is then found at $-2.8477$ Hartree$=77.49$ eV.