Quantum Mechanics using Matrix Methods

Introduction and the simple harmonic oscillator

In this notebook we study some problems in quantum  mechanics using matrix methods. We know that we can solve quantum mechanics in any complete set of basis functions.  If we choose a particular basis, the Hamiltonian will not, in general, be diagonal, so the task is to diagonalize it to find the eigenvalues (which are the possible results of a measurement of the energy) and the eigenvectors. 

In many cases this can not be done exactly and some numerical approximation is needed. A common approach, which is the basis of a lot of quantum chemistry is to take a finite basis set and diagonalize it numerically. The ground state of this reduced basis set will not be  the exact ground state, but by increasing the size of the basis we can improve the accuracy and check if the energy converges as we increase the basis size. We will apply this approach here for an anharmonic oscillator. The reference for this material is Kinzel and Reents, p. 47-51.

We first discuss the exactly solvable case of the simple harmonic oscillator.  The Hamiltonian is given by

(Subscript[H, 0]) = ((p^2)/(2m)) + (1/2)m (ω^2)(x^2)

where p is the momentum, x the position, m the mass and ω the angular frequency of the classical oscillator.  x and p satisfiy the commutation relation [x, p] = i ℏ.  The Hamiltonian can be written in dimensionless form as

(((Subscript[H, 0]) )/(ℏω))= (1/2) ((((p/(Subscript[p, 0]))))^2)+ (1/2)((((x/(Subscript[x, 0]))))^2)

where (Subscript[p, 0])= (Sqrt[ℏm ω]) and  (Subscript[x, 0])= (Sqrt[ℏ/(m ω)]), are the basic momentum and length scales. From now on, we will give the energy in units of ℏω,  x in units of (Subscript[x, 0]), and p in units of (Subscript[p, 0]), so the reduced Hamiltonian is

(Subscript[H, 0]) = ((p^2)/2) + ((x^2)/2).

In these units the commutation relation is [x, p] = i.  In any textbook on quantum mechanics, it is shown that the energy levels are given by

(Subscript[E, n])= n + (1/2),      n = 0, 1, 2, .... .

and the wavefunctions are given by

(Subscript[ψ, n])(x) =((((2^n)n! (Sqrt[π])))^(-1/2)) (e^(-(x^2)/2)) (Subscript[H, n])(x) ,

where (Subscript[H, n])(x) is a Hermite polynomial.

The books also show that it is easier to determine the energy levels using operator methods rather than the Schrödinger equation, and that is the approach that we will take here.  In this approach, one introduces so-called  raising and lowering operators, (SuperDagger[a])and a, related to x and p by

x = (1/(Sqrt[2]))((SuperDagger[a])+ a),      p = (i/(Sqrt[2]))((SuperDagger[a]) - a),

which have the commutation relations

[ a,(SuperDagger[a])] = 1 .

and in terms of which the Hamiltonian is written

(Subscript[H, 0]) = (SuperDagger[a])a + (1/2).

Denoting eigenstates of the Hamiltonian by |n⟩, then one finds

(SuperDagger[a])|n⟩ = (Sqrt[n+1])|n+1⟩,              a|n⟩ = (Sqrt[n])|n-1⟩

and so

(SuperDagger[a])a|n⟩ = n|n⟩,

from which it follows that the energy (Subscript[E, n])is equal to (n + 1/2) as stated above.

The matrix elements of x in the basis |n⟩ are given by

(Subscript[X, nm])= ⟨n|x|m⟩ = (1/(Sqrt[2]))(Sqrt[m+1])(Subscript[δ, n,m+1])+  (1/(Sqrt[2]))(Sqrt[m])(Subscript[δ, n,m-1]) = (1/2)(Sqrt[n+m+1])(Subscript[δ, |n-m|,1]).

Hence we can conveniently define the matrix elements for x as follows:

Clear["Global`*"]

x[n_, m_] := Sqrt[(n+m+1)]/2 /; Abs[n-m] == 1

x[n_, m_] := 0 /; Abs[n-m] != 1

where we have used the symbol /; , which means "assuming that". From these we can generate a matrix for x keeping a basis of "basissize" states:

x[basissize_] := Table[x[n, m], {n, 0, basissize-1}, {m, 0, basissize-1} ]

Let's see what we get for basissize = 4

x[4] //MatrixForm

It is trivial to generate the Hamiltonian matrix of the simple harmonic oscillator, since it is diagonal, i.e.

h0[basissize_] := DiagonalMatrix [ Table[n+1/2, {n, 0, basissize-1} ] ]

h0[4] //MatrixForm

The anharmonic oscillator

Now we make the problem non-trivial  by adding an anharmonic term. We will take it to be proportional to (x^4), i.e.

H = (Subscript[H, 0 ])+ λ (x^4)

It is easy to generate the matrix for H using the matrix obtained above for x and the convenient "dot" notation in Mathematica for performing matrix products:

h[basissize_, λ_] := h0[basissize] + λ x[basissize] . x[basissize] . x[basissize] . x[basissize]

For example, with a basis size of 4 we get

h[4, λ] // MatrixForm

The eigenvalues can also be obtained numerically and then sorted. Here we give a function (with delayed assignment) for doing this:

evals[basissize_, λ_] := Sort [ Eigenvalues [ N[ h[basissize,  λ] ] ] ]

Now we get some numbers. First consider the trivial example of the simple harmonic oscillator, λ = 0, and a basis size of 10.

evals[10, 0]

Naturally we get the right answer, since the matrix is diagonal in this case. Next we consider the anharmonic oscillator, λ = 1.

evals[10, 1]

The energy levels are no longer uniformily spaced. We expect that the lowest few levels will be fairly accurate, but the higher levels will be inaccurate because they feel more severely the truncation of the basis states.

To test this, we plot the eigenvalues for a range of λ, starting with a basis size of 15.

basissize = 15;

p1 = Plot [Evaluate[evals[basissize, λ]], {λ, 0, 1}, PlotRange -> {0, 11} , PlotStyle -> {AbsoluteThickness[2]} , AxesLabel -> {"λ", "E"}]

To see how accurate these are we compare the results for the size-15  basis with those for a size-10  basis.

basissize= 10;

p2 = Plot [ Evaluate [ evals[basissize, λ] ], {λ, 0, 1}, PlotRange -> {0, 10} , PlotStyle -> {{ AbsoluteThickness[2], Hue[0], Dashing[{0.01, 0.02}]}  }, BaseStyle-> {FontSize->14,FontWeight->"Bold",FontColor -> Hue[0] }];

Show[p1, p2, PlotLabel -> "basissize=15 (solid) and 10 (dashed)"]

The lowest levels appear to have converged well, but the higher levels have not. This is not surprising because the highest of the plotted levels are some of the highest levels possible for the size-10  basis, and so the effects of the finite basis size are naturally large in that region. 

Now we focus on the ground state, which is given by the first element of evals, e.g. for basis size 20 and λ = 0.2 we have

basissize = 20;

evals[basissize, 0.2][[1]]

For a given basis size, we can plot the lowest energy level versus λ,

Plot [ evals[basissize, λ] [[1]] , {λ, 0, 1}, PlotRange -> {0.5, 0.8}, AxesLabel -> {"λ", "(Subscript[E, 0])"}]

We see that the energy level starts off, as it must, at the value 1/2 for λ = 0, and then increases with increasing λ.

Note that to plot all the eigenvalues as a function of λ it is necessary to Evaluate the evals function inside the Plot function, whereas to plot just the lowest eigenvalue it is important not to use Evaluate otherwise the wrong eigenvalue is plotted (For some reason, Plot seems to ignore Sort in this case.) The details of the functioning of the Plot command, and when one should Evaluate the function to be plotted, remain rather mysterious to me.

We really wish to check that our calculation of the ground state energy used a large enough basis size in order to be accurate. To check this, we plot the ground state energy versus the inverse of the size of the basis set. Here are some results for λ = 0.2:

ListPlot[Table[{1/basissize, evals[basissize, 0.2][[1]]}, {basissize, 9, 30}], Frame -> True, Axes -> False, PlotStyle -> {PointSize[0.02], Hue[0] }, Epilog-> Line[{{1/30, 0.602405},{1/8,0.602405}}], PlotRange -> All, FrameLabel -> {"1/(basis size)", "(Subscript[E, 0])"} , RotateLabel -> False, PlotLabel -> "λ = 0.2"]

We see that for basis size greater than about 16, the energy has converged well. As you might expect, and you can verify yourself, a larger basis set is needed if λ is larger. 

The eigenvectors can also be found, and these give the linear coefficients of the Hermite polynomials which make up the coordinate space wavefunction. Hence the wavefunction of the anharmonic oscillator can also be obtained by this method.

In another handout we calculated the lowest three energy levels of this anharmonic oscillator in a completely different way, from the Schrodinger equation. In that approach we solve a differential equation rather than diagonalize a matrix. We found the lowest three energy levels to be
         0.602405,      1.95054,      3.5363

We see that the lowest agrees very well with the ground state determined above from matrix methods. To compare all three numbers we now compute the lowest three levels from the matrix method, still with λ = 0.2:

Table[evals[30, 0.2][[n]], {n, 1, 3}]

All three energies agree with the corresponding result from Schrodinger's equation, to within the digits printed.

Double Well Potential

Next we use matrix methods to calculate the lowest energy levels in a double well potential. The Hamiltonian is given by

H = ((p^2)/2) + V(x),   where V(x) =  - ((x^2)/2) + λ ((x^4)/4).

Note that the coeficient of (x^2) is negative. The minima of the potential energy are at

(Subscript[x, 0]) = ± (1/(Sqrt[λ])),      (Subscript[V, min])  ≡  V((Subscript[x, 0])) = -(1/(4 λ)) .

We plot the potential for the case of λ = 0.2.

Plot[-x^2/2 + λ x^4 /4/. λ -> 0.2, {x, -4, 4}]

The new physics in this example is the possibility of tunneling between the two minima. The reader is referred to any quantum mechanics book for more details on tunneling.

Setup

To obtain the Hamiltonian, it is convenient to use the matrix elements of p, given in any quantum mechanics textbook, (as well those of of x which we already have). (To avoid needing p we could alternatively write H = (Subscript[H, 0]) - (x^2) + (x^4)/4 where (Subscript[H, 0])is the simple harmonic oscillator Hamiltonian which is diagonal).

p[j_, k_] := I Sqrt[(j+k+1)]/2 /; j - k == 1

p[j_, k_] := -I Sqrt[(j+k+1)]/2 /; j - k == -1

p[j_, k_] := 0 /; Abs[j-k] != 1

p[basissize_] := Table[p[n,  m],  {n, 0, basissize - 1}, {m, 0, basissize - 1} ]

h[n_, λ_] := p[n] . p[n] / 2 - x[n] . x[n] / 2 + λ x[n]. x[n] . x[n] . x[n]/4

The function evals is the same as before.

λ = 0.2

Consider first λ = 0.2.

We can find where the two minima are either analytically, which gives (Subscript[x, 0]) =± 1 / (Sqrt[λ]), and (Subscript[E, min]) = -1/4 λ, as stated above, or use Mathematica:

FindMinimum [ -x^2/2 + λ x^4 /4/. λ -> 0.2, {x, 2}]

We extract the energy levels for n (the basis size) equal to 10 and 20

evals[10, 0.2]

evals[20, 0.2]

Next we plot the lowest energy level versus 1/n.  It seems to be converging.

ListPlot[ 
Table[{1/size, evals[size, 0.2][[1]] }  , {size, 10, 30, 4}] , PlotStyle -> {Hue[2], PointSize[0.02]}, Axes -> False, Frame -> True,
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0])"}, RotateLabel -> False,
PlotLabel -> "λ = 0.2"]

Now we want to plot more than one energy level for different values of n.

We construct evalslist, a two-dimensional  list containing all the eigenvalues for different values of n between 10 and 30,

evalslist = Table[evals[n, 0.2], {n, 1, 30} ];

and from this plot the two lowest levels versus 1/n (note we plot each set of data twice, once to get the points without the lines, and once to get the lines without the points; this is done by the Joined -> {True, False, Tre, False} option.)

ListPlot[  {
Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}]},
Axes -> False, Frame -> True, Joined -> {True,False, True, False},PlotStyle -> {{Thick,Red}, {Red,PointSize[Large]}, {Thick,Blue}, {Blue,PointSize[Large]}},
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0])"}, RotateLabel -> False,
PlotLabel -> "λ = 0.2"]

As expected, both levels converge as the basis size increases.

λ = 0.1

Next we consider a smaller value, λ = 0.1, for which the minima are deeper.

Clear[λ];Plot[-x^2/2 + λ x^4 /4/. λ -> 0.1, {x, -4.7, 4.7}]

minpot = FindMinimum [ -x^2/2 + λ x^4 /4/. λ -> 0.1, {x, 3.5}]

When the potential is very deep the two the particle will only tunnel very slowly from one well to the other. The two lowest energy levels are split by this small tunneling splitting, and the eigenstates are the symmetric and anti-symmetric  combinations of the wavefunction for the particle localized in the ground state of the left well and the right well. If we represent the potential at the bottom of each well as a parabola then the lowest energy of each of these states is given by the simple harmonic ground state energy above (Subscript[V, min]), i.e.

(Subscript[ϵ, 0])= (Subscript[V, min]) + (1/2)(Subscript[ω, min])

where  (Subscript[ω, min])is classical freqency of oscillation about the minimum, i.e.  (Subscript[ω, min])=  (Sqrt[V''((Subscript[x, 0]))])  = (Sqrt[2]). Hence

(Subscript[ϵ, 0])=- (( 1)/(4 λ)) + (1/(Sqrt[2])).

We  expect then, that for deep wells, the two lowest energy levels will be given, to a good approximation by

(Subscript[E, 1])≃ (Subscript[ϵ, 0])  - ((Subscript[Δ, 0])/2),  (Subscript[E, 2])≃ (Subscript[ϵ, 0])  + ((Subscript[Δ, 0])/2)   (deep wells)

where (Subscript[Δ, 0]) is the (small) tunneling splitting.

For λ = 0.1,  (Subscript[ϵ, 0])is given by

ϵ0 = ((-1)/(4λ)) + (1/(Sqrt[2]))  /. λ -> 0.1

Now we compare this with the actual energy of the two lowest levels.

evalslist = Table[evals[n, 0.1], {n, 1, 30} ];

ListPlot[  {
Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}]},
Axes -> False, Frame -> True,Joined -> {True, False, True, False},
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0]), (Subscript[E, 1])"}, RotateLabel -> False,
PlotLabel -> "λ = 0.1"]

We see that the two lowest levels are close together as expected (the difference being the tunneling splitting), and the energies are fairly close to (Subscript[ϵ, 0])= -1.793 ... , though they are not precisely equal to that value. The difference arises because, for this value of λ,  the wells are not deep enough to accurately  represent their behavior about the minima as  parabolas.

Now we also include the next two levels. These are also fairly close together, but not as close as the lowest two because, having more energy, they can tunnel more easily from one well to the other.

ListPlot[  {
Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[3]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[3]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[4]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[4]]}, {n, 10, 30}]},
Axes -> False, Frame -> True, Joined -> {True,False,True,False,True,False,True,False},
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0]),(Subscript[E, 1]),(Subscript[E, 2]),(Subscript[E, 3])"}, RotateLabel -> False,
FrameLabel -> "λ = 0.1"]

λ = 0.05

Finally we take  λ to have the smaller value of  0.05, which  makes the wells deeper still:

Clear[λ]; Plot[-x^2/2 + λ x^4 /4/. λ -> 0.05, {x, -6.5, 6.5}]

We calculate (Subscript[ϵ, 0])for this value of λ:

minpot = FindMinimum [ -x^2/2 + λ x^4 /4/. λ -> 0.05, {x, 3.5}]

ϵ0 = ((-1)/(4λ)) + (1/(Sqrt[2]))  /. λ -> 0.05

Now we plot the lowest two levels for different basis sizes:

evalslist = Table[evals[n, 0.05], {n, 1, 30} ];

lw2 = ListPlot[  {
Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],Table[ {1/n, evalslist[ [n]][[1]]}, {n, 10, 30}],
Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}], Table[ {1/n, evalslist[ [n]][[2]]}, {n, 10, 30}]},
Axes -> False, Frame -> True, Joined -> {True,False,True,False},
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0]), (Subscript[E, 1])"}, RotateLabel -> False,
PlotLabel -> "λ = 0.05"];

pe0 = ListPlot[Table[ {1/n, ϵ0}, {n, 10, 30}],PlotStyle -> {AbsoluteThickness[2],Dashing[{0.01, 0.02}], Hue[0.65]},  Joined -> True ];

Show[lw2, pe0]

The dashed line is (Subscript[ϵ, 0 ])= -4.292893218813452`(the value obtained above assuming that potential can be represented as a parobola about each minimum and neglecting tunneling) which we see is very close to the energy of the ground state and first excited state, as expected because the wells are very deep so the approximations made to get the value -4.29289 are accurate.

Finally we show the lowest 8 levels. They are grouped into four pairs as expected, with each pair corresponding to a symmetric and anti-symmetric combination of wavefunctions localized  in the left and the right well.

lw8 = ListPlot[  Flatten[{Table[ 
Table[ {1/n, evalslist[ [n]][[i]]}, {n, 10, 30}],{i, 1, 8}],Table[ 
Table[ {1/n, evalslist[ [n]][[i]]}, {n, 10, 30}],{i, 1, 8}]}, 1],
Axes -> False, Frame -> True,Joined ->Flatten[ {Table[True,{8}], Table[False ,{8}]}],
FrameLabel -> {"1/(basis size)", "(Subscript[E, 0]),..,(Subscript[E, 7])"}, RotateLabel -> False,
PlotLabel -> "λ = 0.05", PlotRange -> {-4.8, 0.5} ];

pe1 = ListPlot[Table[ {1/n, ϵ0+ (Sqrt[2])}, {n, 10, 30}],PlotStyle -> {AbsoluteThickness[2],Dashing[{0.01, 0.02}], Hue[0.35]},  Joined -> True];

pe2 = ListPlot[Table[ {1/n, ϵ0+ 2(Sqrt[2])}, {n, 10, 30}],PlotStyle -> {AbsoluteThickness[2],Dashing[{0.01, 0.02}], Hue[0.85]}, Joined -> True];

Show[lw8, pe0, pe1,pe2]

For the case of a really deep well, the lowest levels would correspond to simple harmonic oscillator levels in each of the two wells with energy

(Subscript[ϵ, n])  = (Subscript[V, min])+ (n + (1/2)) (Subscript[ω, min])

where here we have (Subscript[V, min])=-1/(4 λ) and (Subscript[ω, min])= (Sqrt[2]).  For each  n, tunneling causes the two degenerate levels to split by a small amount (Subscript[Δ, n]), so

(Subscript[E, 2n+1])≃ (Subscript[ϵ, n]) - ((Subscript[Δ, n])/2) ,  (Subscript[E, 2n+2])≃ (Subscript[ϵ, n]) + ((Subscript[Δ, n])/2)   (n = 0, 1, 2, ...)     (very deep wells).

In the above figure, the dashed lines show (Subscript[ϵ, n]) for n=0, 1 and 2.   We see that the lowest pair of levels ((Subscript[E, 0]),(Subscript[E, 1])) is close to the expected value ((Subscript[ϵ, 0])), but higher pairs differ progressively more and more, indicating that the potential cannot be accurately represented by a parabola at these energies.