Variational Methods

The variational technique represents a completely different way of getting

approximate energies and wave functions for quantum mechanical systems. It
is most frequently used to compute the ground state, but can be extended to
compute the low lying excited states. The technique involves guessing a reason-
able, parametric form for a trial ground state wave function. An example of such
a parametric form for a symmetric well ground state centered about the origin
might be a Gaussian distribution (simple harmonic oscillator ground state) of
the form:

 1/2
a 2
ψ̃(x) = e−ax /2 (1)
π

The adjustable parameter for this wave function is a which is related to the
inverse of the width of the wave function. We will argue later, that choosing
a trial wave function such as the harmonic oscillator ground state which is the
exact solution for another potential is frequently a wise choice since it eliminates
considerable drudge work. The harmonic oscillator ground state is often a good
choice for one dimensional square wells, and the ψ100 (r) hydrogen ground state
is often a good choice for radially symmetric, 3-d problems.

The variational procedure involves adjusting all free parameters (in this case
a) to minimize Ẽ where:

Ẽ =< ψ̃|H|ψ̃ > (2)

As you can see Ẽ is sort of an expectation value of the actual Hamiltonian using
the trial wave function. The minimum Ẽ is generally an excellant estimate of
the true ground state energy, and parameter choice which minimizes Ẽ produces
a trial wave function which sometimes reasonably approximates the true ground
state wave function.

Why does it work?

1
 Because of the completeness of bound state wave functions we can expand
any function, including our trial wave function, as a linear combination of the
true bound state wave functions which we’ll write as ψn (x) , n ∈ 1, 2, 3.... We
thus know that our trial wave function can be written as:


ψ̃ = b n ψn (3)
n

where bn are the relative amplitudes for the trial wave function. Since the ψn are
bound state wave functions, they must be eigenstates of the Hamiltonian with
eigenvalues given the the bound state energies or H ψn = En ψn . Let us evaluate
Ẽ in terms of the bn expansion coefficients.

   
Ẽ = ψ̃ |Ȟ| ψ̃ =  bm ψm |Ȟ| b n ψn 
m n

  
= b∗m bn En ψm |ψn  = |bn |2 En (4)
m n n

where the last equality follows from the orthonormality of the bound state wave
†
functions: < ψm |ψn >= 0 unless m = n in which case < ψm |ψm >= 1. Let
us assume that there is a choice of parameter a which allows the the trial wave
function to closely approach the actual wave function. As ψ̃ → ψ1 it must be
true that bn → 0 for all n = 1. In this limit Eq.(4) tells us that Ẽ → E1 .

But how do we know that the trial wave function ψ̃ which provides the best
estimate of E1 is the one with the lowest value for Ẽ? If we can show that Ẽ ≥ E1
for all trial wave wave functions, it becomes clear that the lowest possible Ẽ must
represent the best guess for E1 .

† Eq. (4) is also obvious from the measurement axiom since |bn |2 gives the probability of
measuring En .

2
 The proof of this inequality is quite simple. We note that for all n, En ≥ E1
since E1 is the ground state or lowest energy state.

 
  
Ẽ = |bn |2 En ≥ |bn |2 E1 ≥ |bn |2 E1 ≥ E1 (5)
n n n

where the last inequality follows from condition that the trial wave function is
normalized:
 

|bn |2 =1
n

Since the true energy is always lower than the Ẽ, the minimum Ẽ for all trial
wave functions is clearly the best estimate.

To the extent that the trial wave function is close to the wave function
Ẽ → E1 . All we can actually prove is that Ẽ ≥ E1 but experience shows that
reasonable trial wave functions typically result in extremely close estimates to
the ground state energy. Often even extremely crude trial wave functions give
good ground state energy estimates.

Variational Technique Example

As an example of the variational technique, consider solving for the ground

state of the double ramp well illustrated below using a the harmonic ground

3

a 1/2  2

state: ψ0 = π exp − ax2 where a is the adjustable parameter which is
inversely related to the width of the wave function.

V(x)
ψ~ = exp(-a 2 x 2 /2)
F|x| F | x|

1/a

Our job is to evaluate H = T  +V . To evaluate V  we will need integals

 +∞ 2
of the form Jn = 0 dy y n e−ay . These differ from our previous integrals that
appear in the Bound States chapter since they involve odd powers of x and are
integrated from 0 → ∞ rather than from −∞ → ∞. We can get any ”odd”
integrals we need by starting with

1 ∂Jn
J1 = a−1 and the recursion relation Jn+2 = −
2 ∂a

We begin with

  ∞
 
a  a F
V  = F |x| = 2 dxF x exp −ax = 2
2
F J1 = √ a−1/2
π π π
0

We next compute the T . This can in principle be done by brute force

  2  
h̄2 ∂ 
T  = ψ0  2  ψ0
2m ∂x

It is apparent, that T  only involves derivatives and integrals of the trial wave
function and thus can only depend on the shape parameter a and the mass of

4
the particle m and will thus be independent of the potential well. We can thus
exploit the fact that ψ0 is the ground state of a harmonic oscillator which allows
us to compute the kinetic energy very easily by the virial theorem for a harmonic
oscillator wave function: T  = Eo/2 = h̄ω/4. But what ω corresponds to our
trial wave function a parameter? Fortunately this is easy since a = mω/h̄ for the
harmonic oscillator wave function. Hence



h̄ω h̄ h̄a h̄2 a
T  = = =
4 4 m 4m

We can now write the Ẽ for the trial wave function

h̄2 a F
Ẽ = H = + √ a−1/2
4m π

The formation of a minimum Ẽ is very easy to understand. As a is increased, the

ψ̃0 becomes narrower and more peaked. Hence its double derivative or kinetic
energy increases. At the same time, the narrowing of ψ̃0 means that ψ̃0 probes
closer to the bottom of the well where the potential is lower. There is clearly a
happy balance between lowering the kinetic energy at the expense of raising the
potential.

Minimizing with respect to a we have:

 
2/3
2
∂ Ẽ ∂ h̄ a F h̄2 F −3/2 2mF
= + √ a−1/2 =0→ = √ a∗ → a∗ = √
∂a ∂a 4m π 4m 2 π h̄2 π

Inserting a∗ into our expression for Ẽ and performing some algebra we have:

1/3
2/3
27 h̄F
Ẽ = √
16π m

It would certainly be great to check our variational estimate for ground state
energy against the true ground state energy. Unfortunately I don’t know the exact

5
answer. We can get some confidence in the answer by comparing the inflection
point of our ψ̃0 solution to the classical turning point for the two ramp system.
If we had the true solution, the inflection point would equal the turning point.
Recall at the classical turning point E = V and ∂ 2 ψ/∂x2 = 0. The inflection
point of our trial wave function occurs at:

∂ 2 −ax2 /2 2 2

2
e = (−a + a2 x2 ) e−a x /2 = 0
∂x

√ 1/3  1/3
1 πh̄2 h̄2
→ xinflect =√ = = 0.961 (6)
a 2mF mF

Lets next compute the x coordinate of the classical turning point:

1/3
2/3
27 h̄F
Ẽ = √ = V = F xturn
16π m

1/3  1/3  1/3

27 h̄2 h̄2
or xturn = = 0.813
16π mF mF

Both the classical turning point and inflection point are in reasonable agreement.
It would be very surprising if they were in perfect agreement since a Gaussian
wave function is the correct ground state for an V (x) = K x2 /2 potential which
goes to infinity much faster than the linear potential V (x) = F |x| being consid-
ered here.

Just a quick review of the variational strategy. The best bet is to pick a
trial wave function that that is the normalized, ground state of a physical system
that solves a potential similar to that of the novel problem at hand. The virial
theorem will then allow one to quickly write down the kinetic energy part of the
Hamiltonian and avoid all of the work of taking the double derivative of ψ̃ and
then integrating. One is left with doing a single integral over the new potential
in order to set up Ẽ and the minimization can frequently be done analytically.

6
One could get better estimates of the true ground state by allowing for more
adjustable parameters. For example in this problem, one could add in another
ground state oscillator contribution, with another adjustable shape parameter
and a third parameter that gives the relative contribution of each ground state
oscillator wave function.

Perturbation versus Variation

Recall that the first order perturbation estimate of the energy when a “small”
perturbation ∆V is added to a Hamiltonian with a known solution Ȟo is En ≈
En0 + n|∆V |n where |n is the unperturbed state that satisfies Ȟo |n = En0 |n.
We can cast this solution in a form that resembles the variational method:

Ẽn = n|Ȟ|n = n|Ȟo + ∆V |n (7)

Now since Ȟo |n = En0 |n we have n|Ȟo|n = En0 we can write Eq. (7) in the
form:

Ẽn = ψ̃|Ȟ|ψ̃ = n|Ȟ|n = En0 + n|∆V |n where ψ̃ ≡ |n (8)

The interpretation of Eq. (8) is that first order perturbation theory is equivalent
to a variational estimate where the trial wave function is the solution to the un-
perturbed Hamiltonian. The usual case is we are trying to find the ground state
using the variational technique and as discussed above, we always overestimate
the true ground state energy. Eq. (8) then tells us that since first order pertur-
bation theory is equivalent to a variational estimate, the first order perturbation
theory calculation is always an overestimate of the true ground state energy. If
we allow a shape parameter in ψ̃ instead of fixing it to |n we should nearly
always get closer to the true answer than the first order perturbative estimate.

Using ψ100 (r) as the ψ̃(r) in 3-d problems

Often trial wave function based on the hydrogen ground state wave function
is a convenient choice in 3 dimensional problems with spherically sysmmetric

7
potentials.
1 r
ψ̃100 (r) = √ exp(− ) (9)
π a3 a

In Eq. (9) we think of a as our adjustable parameter rather than setting a = a0

for hydrogen or a = a0 /Z for single electron ion. The main convenience of this
form besides its simplicity is the ability to calculate the kinetic energy term Ť .
If we were to set a = a0 , to get the actual hydrogen solutions, we would have the
following relationships between a0 and the ground state energy E1 :

h̄ 1 e2 1
a= , E1 = − mc2 α2 where α = = (10)
αmc 2 4πo h̄c 137

We know that Ť  can really only depend on a, m and h̄ since:

h̄2   2  
T  = ψ̃ ∇ ψ̃
2m

and a is the only variable in ψ̃100 . Since E = −T  as we show shortly from
the virial theorem, it follows that it must be possible to eliminate α and write Ẽ
as a function of a, m and h̄:

2
h̄ 1 1 h̄ h̄2
α= , Ẽ = − mc2 α2 = − mc2 =− (11)
amc 2 2 amc 2ma2

Our trial wave function is a solution to the Coulomb potential and hence we can
use the virial theorem to relate T  to the total energy Ẽ. According to the radial
version of the virial theorem < T >=< r dV
dr > /2. For the case of a 1/r potential

dr > /2 = −V /2. Hence < T > −2 < T >= Ẽ or T  = −Ẽ.

< T >=< r dV
Hence:
h̄2 h̄2  a0 2  a 2
0
T  = −Ẽ = = = −E 1 (12)
2ma2 2ma20 a a

The last equality, which is particularly useful in atomic problems, follows from
the observation that if a = a0 then Ẽ = E1 = −13.6 eV . Another useful result

8
that we will need shortly also follows from the virial theorem.

 
1 1
= (13)
r a

This is based on the observation that:

e2
Ẽ = − (14)
4πo (2a)

which you can obtain by manipulating the results in Eq. (10). You then use the
virial theorem result V  = −2T  and the potential expression for a hydrogen
atom [ V = e2 /(4πo r)] to obtain:

   
e2 V  V  1 e2 1 1
− = T  + V  = − + V  = = → =
4πo (2a) 2 2 2 4πo r r a

Again this will be true of any solution of the hydrogen ground state and will thus
be true for our trial wave function. We are now very well set up to tackle helium.
The Helium Ground State

The problem of estimating the helium ground state energy is becoming a

classic illustrative problem in the variational technique in quantum mechanics
courses since one gets a value close to the experimental answer, and the problem
can be done analytically - albeit using some clever tricks.

For the chemically challenged, helium has two electrons orbiting about a Z =
2 nucleus. The ground state energy is the total energy of these two electrons which
is less than zero. We can measure this energy by measuring the minimum energy
required to remove both electrons to leave a helium nucleus and two electrons
at rest at infinity. Experimentally it takes 24.6 eV to remove one electron from
a neutral He atom leaving singly ionized He+ . Since He+ is a single electron
atom with n = 1 we can reliably compute its energy from the Bohr formula. The

9
ionization energy for the He+ ion is then Z 2 13.6 = 4 × 13.6 = 54.4 eV . Hence it
takes 54.4+24.6 = 79 eV to remove both electrons. This means the experimental
value for the He ground state is −79 eV .

The helium Hamiltonian is given by:


h̄2 2  e2 2 2 1
Ȟ = − ∇1 + ∇22 − + − (15)
2m 4πo r1 r2 |r1 − r2 |

The factors of 2 come from the nuclear charge. You can easily spot the kinetic
energy terms. The fact that there are 3 particles in the final state means that
there is no obvious reduced mass so we work in the limit of an infinite mass
nucleus. The last term represents the Coulomb repulsion energy between the two
electrons.

In the absence of the 1/|r1 − r2 | term, the two electron Hamiltonian is just
Ȟ = Ȟ1 + Ȟ2 , the ground wave function would be just the direct product:

↑1 ↓2 − ↓1 ↑2
ψ(r1 , r2 ) = ψ100 (r1 ) ψ100 (r2 ) √
2

and the energy would just be E = E(1) + E(2) where E(1),(2) are the energies
of a single electron helium ion in the ground state or −54.4 eV for a total of
†
−108.8 eV . Since the true ground state energy of helium is -79 eV the electron-
electron repulsion must be important.

Our approach will be to construct a trial, two electron wave function of the
form: Ψ̃(r1 , r2) = ψ̃100 (r1 ) ψ̃100 (r2 ) where we use the floating a form given by
Eq. (9). If we set a = a0 /Z = a0 /2 and use this to evaluate Ψ̃|Ȟ(r1 , r2 )|Ψ̃ we
will be effectively doing a first order perturbation calculation as described in the
previous section.

† Since the two electrons are identical Fermions, we put the electrons in an antisymmetric
spin state with S = 0.

10
 If we vary a , to get the minimal Ẽ, we are pretty much guaranteed to get a
more accurate answer since the Ẽ will be lower with an adjustable compared to
a fixed a parameter and thus closer to the true value.

If we use a = a0 /2, we would estimate the ground state energy as:

e2 1
Eg = 8E1 + < ψ100 (r1 )ψ100 (r2 )| |ψ100 (r1 )ψ100 (r2 ) >= 8E1 + < ∆V >
4πo |r1 − r2 |
(16)

Although it is surprising how well this approach will work, we clearly seem to
be in a regime where perturbation theory is unreliable given we know the correct
result is -79 eV and the zero order answer is -108.8 eV.

Essentially the bulk of the drudge work is in the calculation of Ψ̃|∆V |Ψ̃.
We will pursue this calculation and then show how we can use the Ψ̃|∆V |Ψ̃ in
a true variational calculation where we put some adjustable parameters in Ψ.
Computing the ∆V integral

I would have guessed that this 6 dimensional integral is intractable like most
integrals but it can actually be done using an elegant technique. The ∆V 
integral is just ∆V  = [e2 /(4π0 )] |r1 − r2 |−1  We can be fairly confident that
|r1 − r2 |−1  = ξ/a where ξ is a dimensionless number since a is the only length
scale involved with our trial wave function. We will work out the integral in full
detail, just to prove it can be done, but the integral is considerably beyond the
skills of most undergraduates, graduates, or Professor’s who aren’t working in
this branch of computational atomic physics. Our conclusion will be ξ = 5/8 but
getting there is half the fun.

To compute < |r1 − r2 |−1 > we write |r1 − r2 |2 in the spherical variables we
will use for the integration.

|r1 − r2 |2 = (r1 − r2 ) · (r1 − r2 ) = r12 − 2r1 · r2 + r22 = r12 − 2r1 r2 cos θ12 + r22

11
where θ12 is the angle between r1 and r2 . We next insert the ψ̃100 wave functions

1
ψ̃100 = √ e−r/a
π a3



e2 1 e−2r1 /a e−2r2 /a e2 1
< ∆V >= d r1 d r2 
3 3
= I
4πo a π2
6
r12 − 2r1 r2 cos θ12 + r22 4πo a π2
6

We note that the ψ100 expression appears four times in the integral (twice in the
bra and twice in the cket) hence some of the surprising numerical factors. We
begin by partitioning the integral as an integral over r1 and r2 .
 
−2r1 /a e−2r2 /a
I= 3
d r1 e I2 (r1 ) , I2 (r1 ) = d 3 r2  (17)
r12 − 2r1 r2 cos θ12 + r22

We note that the inner I2 integral will have a dependence on r1 and no
dependence on r2 since that is being integrated over. Part of this dependence
√
will come from the r1 (1st electron radial coordinates) under the , the other
piece looks like it will come from cos θ12 piece. Now we introduce a famous trick.
We have a double integral. The first integral involves varying r1 over all of its
space. However in performing the inner integral, one is essentially varying r2 over
all space but keeping r1 fixed. Now the trick is since r1 is fixed, why not use it to
define ẑ direction during the evaluation of the I2 integral? If we do this, we can
equate cos θ12 with the cos θ that desribes r2 since if r1 is in the same direction
as ẑ , θ12 will be the usual angle of r2 with respect to the z axis. We then have:

∞ 1 2π
d cos θ
I2 (r1 ) = r 2 dr e−2r/a  dφ (18)
0 −1 r12 − 2r1 r cos θ + r 2 0

The φ integral is just 2π. The θ integral is very easy. Just substitute v =
r12 + r 2 − 2rr1 cos θ. Then d cos θ = dv/(−2rr1) and the limits become v1 =
r12 + r 2 + 2rr1 = (r1 + r)2 and v2 = (r1 − r)2

12
 1 v2
d cos θ −1
Iθ =  → v −1/2 dv
r12 − 2r1 r cos θ + r 2 2rr 1
−1 v1


−1 √ √ (r1 + r) − |r1 − r| 2/r1 if r < r1
= ( v1 − v2 ) = =
rr1 rr1 2/r if r > r1

The last equality is sort of interesting and follows from |r1 − r| is either r1 − r
or r − r1 depending on the magnitude of r versus r1 . The conditional is easy to
insert by breaking up Eqn. (18) into 2 pieces:

r 
1 −2r/a ∞ −2r/a 3  
e e  = πao 1 − 1 + r1 e−2r1 /a
I2 (r1 ) = 4π  r 2 dr + r 2 dr
r1 r r1 a
0 r1

The exponential integrals can be found in any integral table. We can now insert
I2 (r1 ) into Eqn. (17).

 ∞
−2r1 /a −2r1 /a 1   r1 −2r1 /a 
I= 3
d r1 e I2 (r1 ) = 4π a
2 3
r12 dr1 e 1− 1+ e
r1 a
0

Since I integrand has no angular dependence, the angular pieces integrate to 4π.
∞
The remaining r1 integral is very easy using 0 dy y n exp(−by) = n!/bn+1 to
get:



5a2 e2 1 5a2 5 e2
I = 4π a 2 3
= or < ∆V >= × 4π a
2 3
=
32 4πo a6 π 2 32 8a 4π0


 
5 e2 5 e2 a0
< ∆V >= − − =− −
4 4π0 (2a) 4 4π0 (2a0 ) a

5  a0
Thus < ∆V >= − E1 (19)
4 a
If we insert a = a0 /2 we get < ∆V >= −5E1 /2 = +34 eV . We estimate from the
forgoing that the helium ground state energy is Eg < −108.8 + 34 = −74.8 eV .

13
We are thus much closer to the true answer of -79 eV which is remarkable given
either the size of the perturbation or the fact that we used a very naive wave
function which ignores the electron-electron repulsion.

Improving our Bound

We can do better by allowing a to vary. I find it more instructive to write

a = a0 /Z and allow Z = Zeffective to be the adjustable parameter. This is a
nice approach since we can think of the Coulomb repulsion of the other electron,
as creating a “shielded” nuclear charge for a given electron. As such we might
expect an effective nuclear charge of Z < 2 and which will place the electrons
further away that we were assuming before with our Z = 2 Bohr wave functions.

Our goal is to compute Ẽ from Eqn(15) which gives the true helium Hamil-
tonian:

       
e2 1 1   1 
Ẽ = Ψ|Ť1 + Ť2 |Ψ − 2 Ψ  +  Ψ − Ψ  Ψ (20)
4πo r1 r2 |r1 − r2 | 

h̄2 2
where Ť1,2 = − ∇
2m 1,2
You will recognize the last term as < ∆V > which we computed in the previous
section in Eq. (19):

5  a0 5Z
< ∆V >= − E1 = − E1 (21)
4 a 4
We next compute one of the kinetic energy terms:

Ψ(r1 , r2 )|Ť1 |Ψ(r1 , r2 ) = ψ100 (r1 )|Ť1 |ψ100 (r1 )ψ100 (r2 )|ψ100 (r2 ) = ψ100 (r)|Ť |ψ100 (r)

h̄2 2 h̄
2
T1  = = Z = −Z 2 E1
2ma2 2ma20

T1 + T2  = 2T1  = −2 Z 2 E1 (22)

where we have used Eq. (11) to equate E1 with −h̄2 /(2ma20 ) and the fact that
the two electrons are indistinquishable to conclude T1  = T2 .

14
 The same reasoning tells us us 1/r1 = 1/r2 . Using Eq. (13) we have
1/r1 = 1/r2  = 1/a = Z/a0 . Hence the second term of Eq. (20) representing
the attraction of the two electrons to the nucleus is:

   
e2 1 1  e2
− 
2 Ψ + Ψ = −4Z = 8ZE1 (23)
4πo r1 r2 4πo a0

Assembling Eq. (21) - Eq. (23), we have:


5Z 27
Ẽ = −2Z E1 + 8ZE1 −
2
E1 = −2Z + Z E1
2
4 4

We thus obtain a classic concave up parabola (remember E1 is negative) with

a minimum at dẼ/dZ = (−4Z + 27/4)E1 = 0 or Z = 27/16 = 1.69 or Ẽmin =
(36 /27 )E1 = −77.5 eV which is getting pretty close to the experimental value of
-79 eV. Our first order calculation where we set Z = 2 placed us 4.2 eV above
the true value. By allowing Z to be an adjustable parameter in a variational
calculation, we are now only 1.5 eV above the experimental value. We have thus
shrunk the gap between calculation and observation by a factor of 2.8 with very
little additional work.

The Z value is quite sensible as well. Without the additional electron, a given
electron would see Z = 2. If the additional electron provided perfect shielding
we would have Z = 1. A value of Z = 1.69 seems quite reasonable.

Although this problem is solveable analytically , it would not be that difficult

to perform even the 6 dimensional < ∆V > integral on a computer by summation.
Essentially we know that |r1 − r1 |−1  which is the essence of the < ∆V >
calculation must be of the form |r1 − r1 |−1  = ξ/a where ξ is a numerical
constant since a is the only thing in our Ψ̃ that carries dimensions of length. All
of the clever analytic integration did was to establish that ξ = 5/8. It is likely to
be simpler and faster to write such a program to do the calculation numerically
from the get go. But there is a sense of accomplishment in doing the seemingly

15
impossible analytic calculation. I do hope you find you are able to follow the
basic QM that we used in this example.
The H2+ Molecule

In this section we summarize a variational calculation on the binding energy

H2+ where an electron “orbits” about two protons to form the simplest molecule.
This was one of the earliest attempts to understand the nature of a covalent
bond performed within a few years after the development of the Schrod̈inger
Equation. I will leave out many of the mathematical details which are discussed
in Griffiths and other quantum mechanics books such as Quantum Mechanics
by Goswami and Quantum Physics by Gasiorowicz. We are describing a
“covalent” bond. There are also ionic bonds in chemistry.

This calulation is concerned with whether or not the H2+ molecule will even
exist. How does the “orbiting” electron’s attraction to either proton, overcome
the Coulomb repulsion between the two protons to form a bound state i.e. a
state with negative energy? Actually the existence of a negative energy state is
really not enough by itself since the molecular configuration must be energetically
favored over a final state consisting of a hydrogen atom plus a free proton. This
means the system must have a total energy of less than E1 or -13.6 eV. It is
really not all that obvious a priori that the H2+ ion should exist. For example
in the upper part of the below figure we placed the electron half way between
the two protons which are separated by 2a0 . Adding the potential energies we
find a net potential energy of 3E1 = −40.8 eV which is considerably smaller
than E1 = −13.6 eV suggesting that this configuration lies considerably lower in
energy than a hydrogen atom with an additional proton infinitely far away and
at rest. Of course we need to worry about the fact that the elecron has a kinetic
energy which will raise its energy, perhaps by T  = −E1 . It is still plausible
that if the electron were often between the two protons as shown, the hydrogen
molecule might bind even when kinetic energy is considered.

Of course, the electron cannot be precisely located between the two protons

16
but is described by a probability function that extends over all space.

For example the electron could easily be in the position shown in the lower
half of the above figure. In this position the sum of the potential energies is just
E1 which is the same energy as single hydrogen atom with a distant proton at
rest. When kinetic energy is taken into account the electron would bind into
a hydrogen atom rather than forming a molecule with the two protons in this
case. To really know if the H2+ molecule exists we need a legitimate quantum
mechanical calculation which we will base on a variational approach.

Here is the essence of our variational approach. We will place the two protons
a distance R apart and assume they are “nailed” down at these positions. This
will give a positive proton-proton repulsion contribution of

e2 −2ao −e2 2ao

Vpp = = =− E1
4πo R R 4πo 2ao R

†
to the total binding energy. The electron will then have a Hamiltonian given

† It also seems by nailing down the protons we are ignoring their own quantum physics such
as the fact that need to move about to conform to the uncertainty principle. My guess (as
a non-expert) is this sort of kinetic energy < p2 /2M > may be may be ignorable owing to
the large proton mass.

17
by


h̄2 2 e2 1 1
Ȟ = − ∇ − + (24)
2m 4πo r1 r2

where r1 and r2 are shorthands for r1 = |r − R  1 | and r2 = |r − R

 2 |, r is the
 1 , and R
coordinate of the electron where R  2 are the positions of the two protons
2 − R
and R = |R  1 |. Our approach will be to use a trial wave function (Ψ),
compute Ẽ =< Ψ|Ȟ|Ψ > +Vpp and vary R looking for the minimum of Ẽ. If
we agree with the fixed R assumption, the true energy will lie below Ẽ hence if
Ẽ < −13.6 eV the molecule will form.

We will use a linear combination of the atomic orbits (known LCAO) as our
trial wave function based on the hydrogen ground state. Specifically we will use:


Ψ(r) = N ψ100 (|r − R
 1 |) + ψ100 (|r − R
 2 |) = N (ψ(r1 ) + ψ(r2 ))

This Ψ is again motivated for mathematical convenience and physical intu-

ition. We are saying the electron has an equal probability of being in a ground
state hydrogen centered about either proton in the molecule. These wave func-
tions satisfy the Schrödinger equation of a very similar Hamiltonian which cuts
down on the mathematical drudgery. The factor N reflects the fact that we have
a non-trivial normalization. Now it is apparent from symmetry that we want an
equal probability for being attached to either proton but we could also achieve
this with any arbitrary unimodular phase: Ψ = N (ψ(r1) + exp(iδ)ψ(r2 )).

Actually the only two possible choices are exp(iδ) = ±1 because the Hamilto-
nian has r1 ↔ r2 symmetry which we will define in terms of a symmetry operator:
χ̌ψ(r1 , r2) = ψ(r2 , r1 ). This suggests that the wave function can be chosen to be
a simultaneous eigenfunction of χ̌ and x̌. Writing the eigenvalue of χ̌ as λ we
have: χ̌2 ψ(r1, r2 ) = ψ(r1, r2 ) = λ2 ψ(r1 , r2 ) which implies λ2 = 1 or λ = ±1. This

18
means

χ̌ (ψ(r1 ) + exp(iδ)ψ(r2)) = (ψ(r2 ) + exp(iδ)ψ(r1)) = ± (ψ(r1 ) + exp(iδ)ψ(r2))

Solving this we see the only two phase choices are exp(iδ) = ±1. It turns out
that our initial phase choice (δ = 0) gives the best bound. You will explore the
other choice ( δ = π ) in the exercises.

Lets begin with the normalization

N 2 ψ(r1) + ψ(r2) | ψ(r1 ) + ψ(r2 ) = 1 = N 2 (2 + 2 < ψ(r1 )|ψ(r2) >)

We use the fact that

< ψ(r1 )|ψ(r1) >=< ψ(r2)|ψ(r2) >= 1 , < ψ(r1 )|ψ(r2) >=< ψ(r2 )|ψ(r1) >

The first relation follows for the atomic orbital normalization. The second follows
from the fact that the wave functions are real. We often call < ψ(r1 )|ψ(r2) >≡ I
the “overlap” integral. One form for this, where we put one proton at the origin
and the other on the ẑ axis , would be:




1 r |r − Rẑ|
I= d r exp −
3
× exp −
πa3o ao ao

Clearly the overlap integral I can only depend on R and ao and in fact only
depends on the ratio ξ = R/ao. One can show (although its not trivial) that:


ξ2 1
I = exp(−ξ) 1 + ξ + and N =
(25)
3 2(1 + I)

We can now proceed to construct Ẽ using the Hamiltonian given by Eqn.

(24). However, our wave function is constructed from hydrogen wave functions

19
which are eigenfunctions of pieces of the Hamiltonian. Hence for example:

 
h̄2 2 e2 1 e2 1
Ȟψ(r1) = − ∇ ψ(r1) − ψ(r1 ) − ψ(r1 )
2m 4πo r1 4πo r2

The piece of the Hamiltonian within the {} is the Hamiltonian for a hydrogen
 1 acting on a similarly displaced ψ100 (r − R
atom displaced by R  1 ). It will thus
return the usual energy E1 ψ(r1 ). Hence:

e2 1 e2 1
Ȟψ(r1) = E1 ψ(r1 )− ψ(r1) and similarly Ȟψ(r2) = E1 ψ(r2 )− ψ(r2 )
4πo r2 4πo r1

Thus ȞΨ = N Ȟ [ψ(r1 ) + ψ(r2 )] or:


e2 1 1
ȞΨ = E1 Ψ − N ψ(r1 ) + ψ(r2 )
4πo r2 r1

Hence:

Ẽ =< Ψ|Ȟ|Ψ >= E1 < Ψ|Ψ >

       
e2 1 1
−N 2    
ψ(r1) + ψ(r2 )   ψ(r1 ) + ψ(r1 ) + ψ(r2 )   ψ(r2 )
4πo r2 r1

Taking advantage of some r1 ↔ r2 symmetries:

       
e2 1 1
Ẽ = E1 − 2N 2
ψ(r1 )   ψ(r1) + ψ(r1 )   ψ(r2 )
4πo r2 r1

       
e2  ao   ao 
Ẽ = E1 − 4N 2    
ψ(r1)   ψ(r1) + ψ(r1 )   ψ(r2 )
4πo 2ao r2 r1

       
 ao   ao 
Ẽ = E1 +4E1 N (D+X) where D =
2
ψ(r1)   ψ(r1) and X = ψ(r1 )   ψ(r2 )
r2 r1

where D and X are dimensionless positive integrals which can be worked out

20
analytically but are non-trivial:


1 1 −2ξ
D = − 1+ e , X = (1 + ξ)e−ξ where ξ = R/ao (26)
ξ ξ

We note that the two integrals X , I that are based on < ψ(r1 )| anything |ψ(r2) >
fall off as exp(−R/ao) as ξ → ∞ as one might expect since the Bohr ground
state wave functions fall off exponentially as R/ao . We can take D(ξ) to the
 → ∞ by realizing r2 − r1 = R
R = |R|  →R
 thus r2 = r1 + R  and
     a   a
 ao   o 1
ψ(r1)   ψ(r1) → ψ(r1)   ψ(r1) =
o
D= =
r2 R R ξ

which is borne out by Eq. (26).

Another check is to go to the opposite limit where ξ → 0 or R → 0. In this

limit r2 → r1 . Hence
a 
   
 o 1
I(ξ) → ψ(r)|ψ(r) = 1 ; D(ξ) → I(ξ) → ψ(r)   ψ(r) = a0 =1
r r

where we use the result from Eq. (13) that r−1 = a−1
0 . Indeed the forms for all
three functions Eq. (25) and Eq. (26) approach 1 as ξ → 0, although one must
take some care in taking D(ξ) to this limit.

Including the Vpp term, and inserting our expression for N 2 in terms of third
non-trivial integral (I) we have our total energy E

2 D+X
E = Ẽ + Vpp = − E1 + E1 + 2E1
ξ 1+I

Collecting the expressions for the D , X , I integrals we have:

  
E 2 1 − (2/3)ξ 2 e−ξ + (1 + ξ) e−2ξ
f (ξ) = = −1 +
−E1 ξ 1 + (1 + ξ + (1/3)ξ 2) e−ξ

As ξ → ∞ all three integrals D , X , I vanish as does Vpp and E → E1 . This

makes perfect sense since as R → ∞ we are essentially left with a superposition

21
of a hydrogen atoms and a free proton. As R → 0 the Vpp term dominates and
E → +∞. But between these two extremes we have E < E1 with a minimum near
E = −13.6 − 1.76 eV . The 1.76 eV would be our bound on the molecular binding
energy for the H + ion. Because this is a variational calculation we know this ion
is more tightly bound and in fact has a binding energy of 2.8 eV implying an
actual energy lower than E. Our estimate was not a great estimate of the binding
energy, indicating our trial wave function wasn’t that realistic. Our minimizing
R is at 1.27 Ȧ whereas experimentally the number is closer to 1.06 Ȧ. We could
do much better by putting some adjustability in out trial wave functions. But
at least its plausible that molecules will form as stable entities. Given we are
constructed from molecules its good we proved the simplest one exists! Here is
a crude sketch of ours and a much more sophisticated molecular calculation that
agrees with the experimental values.

energy

.
1.3 A

-13.6 eV R

1.76 eV

2.8 eV

.
1.06 A

In the exercises I ask you to consider binding of an alternative trial wave

function Ψ = N (ψ100 (r1 ) − ψ100 (r2 )) where the normalization is different. I
think you will find that this wave function will not bind. I think the reason is
that with this relative phase, the trial wave function vanishes on the plane where
r1 = r2 and therefore cannot take advantage of the large binding contribution
near this plane. These are called “anti-bonding” orbitals and play an important
role in understanding molecular dynamics in chemical processes.

22
 Although we know our variational calculation is not terribly accurate it does
resemble more accurate calculations of the binding energy versus the atomic
separation R. At some risk, I have crudely represented the energy near the
minimum by the dashed parabola shown above to find an effective spring constant
of k = 500 eV /nm2 where the parabola is given by Ẽ = kx2 /2 − Emin . We
can roughly think of the two protons as interacting through the approximate 1
dimensional simple harmonic oscillator Hamiltonian of the form

−h̄2 ∂ 2 1 2
Ȟ = + kx
2µ ∂x2 2

where x = x1 − x2 is a relative coordinate and ω = k/µ where µ = M/2 is

the reduced mass for two protons. We thus expect a set of vibrational molecular
levels separated by energies of ∆E = h̄ω or

 
2k 2 × 500 eV /nm2
∆E = h̄c = 197eV nm = 0.203 eV
M c2 940 × 106 eV

where we used 940 × 106 eV as M c2 for a proton. Indeed one does typically
get a series of roughly evenly spaced infrared absorption and emission lines at
energies on the order of a few hundredths to a few tenths of an electron volt due

23
to molecular vibrations. Typically , as in this case, the level spacing will get
smaller for the more highly excited states since the well is only approximately
parabolic.

Another feature of molecular aborption or emission spectra are rotational

excitations. We can think of the H2 molecule as a rigid rotator with a moment
of inertia given by µR2 that can be in various L2 states. The kinetic energy will
then be of the form:

( + 1)h̄2 (197 eV nm)2

E = = (1+) = (+1) 2.4×10−3 eV
2µR2 2 12 (940 × 106 eV ) (0.13 nm)2

where we have used the minimum of our admittedly crude potential well calcula-
tion as a bond length of 0.13 nm. In fact measuring the energy level spacing due
to rotational excitations of molecules is a frequently used method for measuring
bond lengths – which of course play a crucial role in constructing chemical models
of molecules.

Hence superimposed on the infrared absorption spectra due to vibrational

modes will be a series of even closer spaced rotational spectrs.

Important Points

1. We discussed the use of a variational technique to compute the approximate

energy of a quantum mechanical ground state. The technique requires
you to guess a “reasonable” trial wave function, ψ̃(x), with at least one
adjustable (shape) parameter. The wisest choice is to pick an actual of
ground state of a related system such as a harmonic oscillator ground state.
We proved the inequality:

Ẽ = < ψ̃|H|ψ̃ > > Eground

In computing Ẽ one must be careful to include both the kinetic and po-
tential energy terms. If one chooses a ground state of a related system the

24
 kinetic energy can be written down very easily using the virial theorem.
The minimum value of Ẽ (as its parameters are varied) generally provides
a good estimate for the true ground state energy of the system.

2. We applied the variational technique to the helium ground state. The

complication to the quantum mechanics of helium is the electron-electron
repulsion term. For this problem we used the product of two shielded
hydrogen atom ground state orbit to act as a two particle trial wave function
and got within 2% of the experimental result and a reasonable result for
the shielded nuclear charge.

3. We computed the binding energy of the H2+ molecule where an electron

orbits about two protons. We did the calculation by placing the two protons
a distance R apart and then tried a trial wave function of the form Ψ ∝
ψ100 (r1 )+ψ100 (r2 ) where r1 is the distance of the electron to one proton and
r2 is the distance of the electron to the other. This is sometimes called a
LCAO approach or linear combination of atomic orbitals. We found that it
was indeed possible to get a variational estimate that had the total energy
of the system (including the Coulomb repulsion between the two protons)
less than -13.6 eV which means the molecule has a smaller energy than
a hydrogen atom and free proton. We varied the R separation between
the two protons but did not vary the LCAO wave function. The molecular
binding energy was rather small and not terribly accurate but we concluded
that molecules exist.

4. Finally we showed how the Ẽ versus R curve can be used to crudely es-
timate bond lengths and effective spring constants that can be used to
estimate the vibrational and rotational energy levels of diatomic molecules.
Typically the vibrational levels absorb infrared photons and the rotational
levels absorb far infrared photons.

25