# 21: Linear Variational Theory

- Page ID
- 40108

Last lecture continued the discussion of variational method approach to approximate the solutions of systems that we cannot analytically solve the Schrodinger equation. The system that motivates us is multi-electron atoms, and we focused on the He atom primarily. This approach requires postulating a trial wavefunction and then calculating the energy of that function as a function of the parameters that describe the trial wavefunction. Then we can minimize the energy as a function of these parameters and the closer the wavefunction "looks" like the true wavefunction (that we do not know), the closer the trial energy matches the true energy (however, the trial energy is ALWAYS higher in energy than the true energy). We went over several example trial wavefunctions for the He atom showing the more complex wavefunctions give better results than the simple ones (including the "Ignorance is Bliss" approximation with an average effective charge). We then transition into the Heisenberg's matrix representation of Quantum mechanics which was the segway to the linear variational method, which addresses trial functions that are a linear combination of a basis functions. We will continue that discussed next lecture.

## Overview (again) of Variational Method Approximation

We can always construct a variational energy for a trial wavefunction given a specific Hamilitonian

\[E_{trial} = \dfrac{\langle \psi_{trial}| \hat{H} | \psi_{trial} \rangle }{\langle \psi_{trial}| \psi_{trial} \rangle} \ge E_{true}\label{7.3.1b}\]

The variational energy is an upper bound to the true ground state energy of a given molecule. The general approach of this method consists in choosing a "trial wavefunction" depending on one or more parameters, and finding the values of these parameters for which the expectation value of the energy is the lowest possible.

Using the variational method approximation, find the ground state energy of a particle in a box using this trial function:

\[| \phi \rangle = N\cos\left(\dfrac{\pi x}{L}\right) \nonumber \]

How does is it compare to the true ground state energy?

**Solution**-
The problem asks that we apply variational methods approximation to our trial wavefunction.

\[ E_{\phi} = \dfrac {\langle \phi | \hat{H} | \phi \rangle} { \langle \phi | \phi \rangle} \ge E_o \nonumber \]

Let's look at the

**denominator:**\[\langle \phi | \phi \rangle = 1 = \int_{0}^{L}N^2\cos^2\Big(\dfrac{\pi x}{L}\Big)\nonumber \]

Performing this integral and solving for N yields

\[N = \sqrt{\dfrac{2}{L}}\nonumber \]

Now let's look at the

**numerator:**The Hamiltonian for a particle in a one dimensional box is \(\hat{H} = \dfrac{-\hbar^2}{2m}\dfrac{d^2}{dx^2}\)

\[ \begin{align*} \langle \phi | \hat{H} | \phi \rangle &= \langle N\cos\Big(\dfrac{\pi x}{L}\Big) \Big| \dfrac{-\hbar^2}{2m}\dfrac{d^2}{dx^2} \Big| N\cos\Big(\dfrac{\pi x}{L}\Big) \rangle\nonumber \\ &= \int_{0}^{L}N\cos\Big(\dfrac{\pi x}{L}\Big)\dfrac{-\hbar^2}{2m}\dfrac{d^2}{dx^2}N\cos\Big(\dfrac{\pi x}{L}\Big)dx\nonumber \\ &= \dfrac{N^2\pi^2 \hbar^2}{2mL^2}\int_{0}^{L}\cos^2\Big(\dfrac{\pi x}{L}\Big) dx\end{align*} \]

where \(N = \sqrt{\dfrac{2}{L}}\). The above equation after the integral becomes

\[\dfrac{\pi^2 \hbar^2}{mL^3}\Big(\dfrac{L}{2}\Big)\nonumber \]

Now the variational energy for this trial wavefunction is

\[E_\phi = \dfrac{\pi^2 \hbar^2}{2mL^2}\nonumber \]

This is equal to the ground state energy of the particle in a box that we calculated from the Schrodinger equation using

\[\psi = \sqrt{\dfrac{2}{L}}\sin(\dfrac{n\pi x}{L})\nonumber \]

This is not a good trial wavefucntion since we are not varying any parameter, and moreover it is the eigenstate of energy. So the guess (trial wavefunction) is accurate... this almost never happens.

## Let's look at a Different Type of a Trial Wavefunction

Let's expand \(\phi\) into a linear combination of basis functions:

\[ | \psi_{trial} \rangle = \sum_j^N a_j | \phi_j \rangle \label{Ex1}\]

and similarly

\[ \langle \psi_{trial} | = \sum_j^N a_j^* \langle \phi_j | \label{Ex2}\]

In these cases, one says that a 'linear variational' calculation is being performed.

The set of functions {\(\phi_j\)} are called the 'linear variational' basis functions and are usually selected:

- to obey all of the boundary conditions that the exact state \(| \psi _{trial} \rangle\) obeys,
- to be functions of the the same coordinates as \(| \psi _{trial} \rangle\), and
- to be of the same symmetry as \(| \psi _{trial} \rangle\).

Beyond these conditions, the {\(\phi_j\)} are nothing more than members of a set of functions that are convenient to deal with (e.g., convenient to evaluate Hamiltonian terms elements \(\langle \phi_i|H|\phi_j \rangle\) that can, in principle, be made complete if more and more such functions are included in the expansion in Equations \(\ref{Ex1}\) and \(\ref{Ex2}\) (i.e., increase \(N\)).

For example, For the Hydrogen atom wavefunctions, \(\phi\) could be expanded into a linear combination of Scalable Gaussian functions

\[|\phi_{trial} \rangle = \sum_{j=1}^{N} c_{j} e^{-\alpha_{j} r^2}\label{6B}\]

or for \(H_2\), (this is call the Linear Combination of Atomic Orbitals approximation discussed in detail in later sections)

\[| \phi_{trial} \rangle = c_1 \psi_{1s_A} + c_2 \psi_{1s_B}\]

The goal is to solve for the set of all \(c\) values that minimize the energy \(E_{trial}\).

\[ E_{trial} = \dfrac{ \langle \psi _{trial}| \hat {H} | \psi _{trial} \rangle}{\langle \psi _{trial} | \psi _{trial} \rangle} \label{7.1.8}\]

Substituting Equation \ref{Ex2} into Equation \ref{7.1.8} involves addressing the numerator and denominator individually. For the numerator, the integral can be expanded thusly:

\[\begin{align} \langle\psi_{trial} |H| \psi_{trial} \rangle &= \sum_{i}^{N} \sum_{j} ^{N}a_i^{*} a_j \langle \phi_i|H|\phi_j \rangle. \\[5pt] &= \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j \langle \phi_i|H|\phi_j \rangle. \label{MatrixElement}\end{align}\]

We often rewrite the following integral in Equation \ref{MatrixElement} (as a function of the basis elements, not the trial wavefunction) as

\[ H_{ij} = \langle \phi_i|H|\phi_j \rangle\]

So the numerator of the right side of Equation \ref{7.1.8} becomes

\[\langle\psi_{trial} |H| \psi_{trial} \rangle = \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j H_{ij} \label{numerator}\]

Similarly, the denominator of the right side of Equation \ref{7.1.8} can be expanded

\[\langle \psi_{trial}|\psi_{trial} \rangle = \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j \langle \phi_i | \phi_j \rangle \label{overlap}\]

We often rewrite the following integral in Equation \ref{overlap} (as a function of the basis elements, not the trial wavefunction) as

\[ S_{ij} = \langle \phi_i|\phi_j \rangle \]

where \(S_{ij}\) are overlap integrals between the different {\(\phi_j\)} functions. So Equation \ref{overlap} is

\[\langle \psi_{trial}|\psi_{trial} \rangle = \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j S_{ij} \label{denominator}\]

There is no explicit rule that the {\(\phi_j\)} functions have to be **orthogonal** and **normal** functions, although they often are selected that was for convenience. Therefore, *a priori*, \(S_{ij}\) does not have to be \(\delta_{ij}\).

Substituting Equations \ref{numerator} and \ref{denominator} into the variational energy formula (Equation \ref{7.1.8}) results in

\[ E_{trial} = \dfrac{ \displaystyle \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j H_{ij} }{ \displaystyle \sum_{i,\,j} ^{N,\,N}a_i^{*} a_j S_{ij} } \label{Var}\]

## Minimizing the Variational Energy

The optimum coefficients are found by searching for minima in the variational energy landscape spanned by varying the \(\{a_i\}\) coefficients (Figure \(\PageIndex{1}\)).

We want to minimize the variation energy with respect to the linear coefficients \(\{a_i\}\), which requires that

\[\dfrac{\partial E_{trial}}{\partial a_i}= 0\]

for all \(i\).

For the "normal" variational method discussed previously, the variation energy can be a nonlinear function of one or more parameters. Minimizing the variation energy therefore may require linear or non-linear regression. However, the linearity of the trial function and the nature of the variational energy allows for more simplistic linear regression.

The expression for variational energy (Equation \ref{Var}) can be rearranged

\[E_{trial} \sum_{i,\,j} ^{N,\,N} a_i^*a_j S_{ij} = \sum_{i,\,j} ^{N,\,N} a_i^* a_j H_{ij} \label{7.2.9}\]

Differentiating both sides of Equation \(\ref{7.2.9}\) for the \(k^{th}\) coefficient gives,

\[ \underbrace{ \dfrac{\partial E_{trial}}{\partial a_k} \sum_{i,\,j} ^{N,\,N} a_i^*a_j S_{ij}+ E_{trial} \sum_i \sum_j \left[ \dfrac{ \partial a_i^*}{\partial a_k} a_j + \dfrac {\partial a_j}{\partial c_k} a_i^* \right ]S_{ij} }_{\text{product rule}}= \sum_{i,\,j} ^{N,\,N} \left [ \dfrac{\partial a_i^*}{\partial a_k} a_j + \dfrac{ \partial a_j}{\partial a_k}a_i^* \right] H_{ij} \label{7.2.10}\]

Since the coefficients are linearly independent (Equation \ref{Ex1})

\[\dfrac{\partial a_i^*}{ \partial a_k} = \delta_{ik}\]

and

\[S_{ij}^* = S_{ji}\]

and also since the Hamiltonian is a Hermetian Operator

\[H_{ij}^* =H_{ji}\]

then Equation \(\ref{7.2.10}\) simplifies to (after collapsing the double sum to a single sum and bringing together the two other sums into one)

\[ \dfrac{\partial E_{trial}}{\partial a_k} \sum_i \sum_j a_i^*a_j S_{ij}+ 2E_{trial} \sum_i a_i S_{ik} = 2 \sum_i a_i H_{ik} \label{7.2.11}\]

Hermitian operators are operators that satisfy the general formula

\[ \langle \phi_i | \hat{A} | \phi_j \rangle = \langle \phi_j | \hat{A} | \phi_i \rangle \label{Herm1}\]

If that condition is met, then \(\hat{A}\) is a **Hermitian operator.** For any operator that generates a real eigenvalue (e.g., observables), then that operator is Hermitian. The Hamiltonian \(\hat{H}\) meets the condition and a Hermitian operator. Equation \ref{Herm1} can be rewritten as

\[A_{ij} =A_{ji}\]

where

\[A_{ij} = \langle \phi_i | \hat{A} | \phi_j \rangle\]

and

\[A_{ji} = \langle \phi_j | \hat{A} | \phi_i \rangle\]

Therefore, when applied to the Hamiltonian operator

\[\boxed{H_{ij}^* =H_{ji}.}\]

At the minimum variational energy, when

\[\dfrac{\partial E_{trial}}{\partial a_k} = 0\]

then Equation \(\ref{7.2.11}\) gives

\[ {\sum _i^N a_i (H_{ik}–E_{trial} S_{ik}) = 0} \label{7.2.12}\]

for all \(k\). The equations in \(\ref{7.2.12}\) are call the** Secular Equations**.

## Secular Determinant

From the secular equations with an orthonormal functions (Equation \ref{7.2.12}), we have \(k\) simultaneous secular equations in \(k\) unknowns. These equations can also be written in matrix notation, and for a non-trivial solution (i.e. \(c_i \neq 0\) for all \(i\)), the determinant of the secular matrix must be equal to zero.

\[ { | H–E_{trial}S| = 0} \label{7.2.13}\]

For a 2x2 system (i.e., expanded to three terms) Equation \ref{7.2.13} looks like:

\[\begin{vmatrix} H_{11}-E_{trial}S_{11}&H_{12}-E_{trial}S_{12} \\ H_{12}-E_{trial}S_{12}&H_{22}-E_{trial}S_{22}\end{vmatrix}=0\]

when the determinate is expanded, it will give a polynomial that will have N roots (solutions of E_{trials}).

- The determinant is a real number, it is not a matrix.
- The determinant can be a negative number.
- It is not associated with absolute value at all except that they both use vertical lines.
- The determinant only exists for square matrices (2×2, 3×3, ... n×n). The determinant of a 1×1 matrix is that single value in the determinant.
- The inverse of a matrix will exist only if the determinant is not zero.

### Determinants can be expanded using Minors and Cofactors

The method is called expansion using minors and cofactors. Before we can use them, we need to define them. It is the product of the elements on the main diagonal minus the product of the elements off the main diagonal.

In the case of a 2 × 2 matrix, the specific formula for the determinant is

\[{\displaystyle {\begin{aligned}|A|={\begin{vmatrix}a&b\\c&d\end{vmatrix}}=ad-bc.\end{aligned}}}\]

Similarly, suppose we have a 3 × 3 matrix *A*, and we want the specific formula for its determinant \(|A|\):

\[{\displaystyle {\begin{aligned}|A| &={\begin{vmatrix}a&b&c\\d&e&f\\g&h&i\end{vmatrix}} \\[4pt] &=a\,{\begin{vmatrix}e&f\\h&i\end{vmatrix}}-b\,{\begin{vmatrix}d&f\\g&i\end{vmatrix}}+c\,{\begin{vmatrix}d&e\\g&h\end{vmatrix}}\\[4pt] &=aei+bfg+cdh-ceg-bdi-afh.\end{aligned}}}\]

### Solving the Secular Determinant

To solve this determinate in Equation \ref{7.2.13}, it should be expanded to generate a polynomial (a **characteristic equation**) that can be directly solved with linear methods (i.e., find the roots - different \(E_{trial}\) values that satisfy the secular equations).

If \(|\psi_{trial} \rangle\) is a linear combination of two functions. In math terms,

\[|\psi_{trial} \rangle= \sum_{n=1}^{N=2} a_n |f_n\rangle = a_1 |\phi_1 \rangle + a_2 | \phi_2 \rangle\]

then the secular determinant (Equation \(\ref{7.2.13}\)), in matrix formulation would look like this

\[\begin{vmatrix} H_{11}-E_{trial}S_{11}&H_{12}-E_{trial}S_{12} \\ H_{12}-E_{trial}S_{12}&H_{22}-E_{trial}S_{22}\end{vmatrix}=0\]

###### Solution

Solving the secular equations is done by finding \(E_{trial}\) and putting the value into the expansion of the secular determinant

\[a_1^2 H_{11} + 2a_1 a_2 H_{12}+ a_2^2 H_{22}=0\]

and

\[a_1(H_{12} - E_{trial}S_{12}) + a_2(H_{22} - E_{trial}S_{22}) = 0\]

Equation \(\ref{7.2.13}\) can be solved to obtain the energies \(E\). When arranged in order of increasing energy, these provide approximations to the energies of the first \(k\) states (each having an energy higher than the true energy of the state by virtue of the variation theorem). To find the energies of a larger number of states we simply use a** greater number** of basis functions \(\{\phi_i\}\) in the trial wavefunction (Example \ref{Ex1}). To obtain the approximate wavefunction for a particular state, we substitute the appropriate energy into the secular equations and solve for the coefficients \(a_i\).

Using this method it is possible to find all the coefficients \(a_1 \ldots a_k\) in terms of one coefficient; normalizing the wavefunction provides the absolute values for the coefficients.

The potential well with infinite barriers is can be written thusly:

\[V(x) = \infty \; \text{for} | x| > |L|\]

and

\[V(x) = 0 \; \text{for} \; |x| \leq |L| \]

and it forces the wave function to vanish at the boundaries of the well at \(x=\pm a\). The exact solution for this problem is known and treated previously. Here we discuss a linear variational approach to be compared with the exact solution. We take \(a=1\) and use natural units such that \(\hbar^2/2m=1\).

As the all variational methods problems with a basis set, the trial wavefunction is expanded

\[| \varphi \rangle = \sum_n f_n(x) \]

As basis functions we take simple polynomials that vanish on the boundaries of the well:

\[ \psi_n(x)=x^n(x-1)(x+1)\]

with \(n=0,1,2,3...\)

The reason for choosing this particular form of basis functions is that the relevant matrix elements can easily be calculated analytically.

###### Solution

We start we the overlap matrix:

\[ S_{mn}=\langle \psi_n\vert\psi_m \rangle = \int_{-1}^1 \psi_n(x) \psi_m(x) dx.\]

Working out the integrals, one obtains

\[S_{mn}=\dfrac{2}{n+m+5} - \dfrac{4}{n+m+3} + \dfrac{2}{n+m+1}\]

for when \(n+m\) even, and zero otherwise.

We need to calculate the Hamiltonian matrix elements:

\[\begin{align*} H_{mn} &=\langle \psi_n \vert p^2 \vert \psi_m \rangle = \int_{-1}^1 \psi_n(x) \left(-\frac{d^2}{dx^2} \right) \psi_m(x) dx \\[4pt] &= -8 \left[ \dfrac{1-m-n-2mn} {(m+n+3)(m+n+1)(m+n-1)} \right] \end{align*}\]

for when \(n+m\) even, and zero otherwise.

## Linear Combination of Trial Wavefunctions (e.g., Gaussian functions)

In electron calculations, the trial function approximation can be constructed in terms of the ground-state wavefunctions (remember that even the radial component of higher orbitals decay exponentially when far from the nucleus:

\[\phi = \sum_{j=1}^{N} c_j e^{-\alpha_j r}\label{10A}\]

An alternative to this trial wavefunction is the combination of Gaussian functions:

\[\phi = \sum_{j=1}^{N} c_j e^{-\alpha_j r^2}\label{10}\]

We know the energy of the hydrogen atom, but using a set of \(N\) Gaussian functions gives...

\(N\) | \(E_{min}\) in Hartree units |
---|---|

1 | -0.42441 |

2 | -0.48581 |

3 | -0.49697 |

: | : |

5 | -0.49976 |

: | : |

16 | -0.49998 |

A quick note on units - \(1 Hartree = 2 R_{\infty}\)

These sets of Gaussian functions are used in computer programs for molecular quantum mechanic calculations, because the integrals are much simpler than using the one-electron Hydrogen atomic functions. Besides, for multi-electron atoms, the H-atom functions are not as accurate.

- The hydrogen atom is the only atom with an exact solution.
- Hydrogen wavefunctions are used as the approximation for atomic wave functions in multielectron atoms for atomic wavefunctions in multielectron atoms (aka orbitals).
- The variational principle states that any wavefunction we choose that satisfies the Schrödinger equation will give an energy
**greater**than the true energy of the system. - The variation method provides a general prescription for improving on any wave function with a parameter by minimizing that function with respect to the parameter minimizing that function with respect to the parameter

## Contributors and Attributions

Michael Fowler (Beams Professor, Department of Physics, University of Virginia)