Matrix mechanics

Quantum Mechanics
 K. Konishi, G. Paffuti

Solve numerially the eigenvalue problem H ψ = E ψ  by using the matrix mechanics.
Write a finite-dimensional approximation to the variables $p, q$ by using harmonic oscillator basis functions. Study such an approximation for polynomial potentials (harmonic and non harmonic oscillator) and for the Morse potential

V[x]= A(Exp[-2α x]- 2 Exp[-α x]).

Utilities

ClearAll["Global`*"];
Off[General::"spell",General::"spell1",InterpolatingFunction::"dmval"];

Introduction

Algebraic relations between operators hold in general, independently of their particular representation, i.e. in any concrete realization of Hilbert space. This suggests the possibility of realizing an approximate formulation of quantum dynamics by a suitable approximation of Hilbert space, in particular in finite dimensions spaces. The space can be generated by an arbitrary finite set of independent base vectors. A particular and very useful example is a finite dimensional space generated by the lowest N eigenstates of an harmonic oscillator. In this space matrix operators are realized as matrices (matrix elements of operators) and full dynamics will be recovered in the limit N→∞.

Let  |n> the n-th energy eigenstate of an harmonic oscillator of frequency ω. It is known (see text) that matrix elements for position operator and momentum are given by

<n| x | n-1> = <n-1| x |n> = (Sqrt[(ℏ/(mω))])(Sqrt[(n/2)]);
<n-1|p|n> = -i (Sqrt[m ω ℏ])(Sqrt[(n/2)])  = (SuperStar[(<n |p|n-1>)])

all other elements being zero. From these relations it is easy, by matrix  multiplication, write down matrix elements for powers of x and p. In particular

<n| (x^2) | n> =  (ℏ/(mω))((2n+1)/2);
<n| (x^2) | n+2> = <n+2| (x^2) |n> = (ℏ/(mω))(Sqrt[(((n+1)(n+2))/4)]);
<n| (x^4) | n> =  ((((ℏ/(mω))))^2)(3/4)(2(n^2)+2n+1);
<n| (x^4) | n+2> = <n+2| (x^4) |n> = ((((ℏ/(mω))))^2)((2n+3)/2)(Sqrt[(n+1)(n+2)]);
<n| (x^4) | n+4> = <n+4| (x^4) |n> = ((((ℏ/(mω))))^2)(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)]);
<n| (p^2) | n> =  ℏ mω((2n+1)/2);
<n| (p^2) | n+2> = <n+2| (p^2) |n> = ℏ mω(Sqrt[(((n+1)(n+2))/4)]);

all other elements being zero. In studying harmonic oscillators it is often useful to use annihilation and creation operators (see text).

a = (Sqrt[((m ω)/(2ℏ))])( x + i (p/(m ω)));    (SuperStar[a]) = (Sqrt[((m ω)/(2ℏ))])( x - i (p/(m ω)))

Their non zero matrix elements are

<n-1| a | n> = (Sqrt[n]) ;    <n+1| (SuperStar[a]) |n> = (Sqrt[n+1]) .

In this notebook we write a function XPmatrices which compute matrices above up to a maximum value of n, we will verify approximately algebraic relations between operators and we will use these finite dimensional matrices to approximate some simple systems.

Finite dimensional matrices

The function XPmatrices needs a parameter, dim, which fixes the dimension of the matrices. The usual principal number of harmonic oscillator starts from zero, so dim - 1 = nMax. The matrices are constructed using the Mathematica function SparseArray.

The function XPmatrices admits an optional parameter parity: if one study only even operators only states with the same parity are to be considered, this halves the dimension of the matrices.  parity = 0 means even states, parity = 1 odd states. Any other value or the absence of this parameter compute matrices for all states. The output of the function is a list of matrices, respectively:
1) p, (p^2),q, (q^2),(q^4),  a, (SuperStar[a]) for all states computations
2) (p^2),(q^2), (q^4)  for fixed parity computations.

Function SparseArray 
Most elements of our matrices are 0. SparseArray gives an efficient way to store such objects, keeping only elements different from 0. Mathematica help gives detailed informations on the subject, see also Advanced Documentation. 
Let us only remember that a standard way to fill the matrix is through patterns:  eg. {i_,i_} ->  -2 means that diagonals elements are -2,                                   the pattern {i_,j_}/;j-i== 1 means: for each i, j under the condition j-i = 1 the element is ... i.e. we define upper main diagonal etc.
For the matrix matV we have defined the whole set of diagonal indices and then apply the rule indDiag-> V[x]. This is faster then {i_,i_} -> V[x[[i]]], as the reader can check.
The function Normal transforms sparsearray structures in usual matrices.
The function SparseArray is available only with Mathematica 5. Users with previous version of the software can define matrices using Table, and / or other Mathematica functions. Eg. (for a discretized laplacian operator)
<<LinearAlgebra`MatrixManipulation`;
D2 = ZeroMatrix[nP];
Table[D2[[i,i]] = -2,{i,1,nP}];
Table[D2[[i,i+1]]= 1,{i,1,nP-1}];
Table[D2[[i+1,i]]= 1,{i,1,nP-1}];

Units and assumptions

In numerical computations it is often practical to use units such that  ℏ = m = ω = 1. If we keep these constants in our calculations it is convenient to state once for all that are positive, this make easier some algebraic simplifications. Below n will denote always a non negative integer.

fixnaturalUnits[naturalUnits_]:= If[naturalUnits== True,m=1;ℏ=1;ω=1;,
Clear[ω,m,ℏ];
old$Assumptions = $Assumptions;$Assumptions = {n∈Integers && n≥ 0,m >0, ω > 0,ℏ >0};];

The function XPmatrices

The function is simply a transcription in Mathematica of relations ()-()

XPmatrices[dim_,parity_Integer:-1]:=
Module[{x2array,p2array,x4array,x2,p2,x4,x1,x1array,p1,p1array}, 
Clear[i,n];
singlepower= False;
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;singlepower = True];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[singlepower == True,
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
p1= SparseArray[{{i_,k_}/;k-i==1->-I(Sqrt[(n+1)/2])},{dim,dim}];
p1 = (p1-Transpose[p1]);
aOp = SparseArray[{{i_,k_}/;k-i==1->(Sqrt[n+1])},{dim,dim}];
adagOp = Transpose[aOp];];
If[parity == 0 || parity == 1,
Normal[#1]&/@{m ω ℏ p2,(ℏ/(m ω))x2,((((ℏ/(m ω))))^2)x4},
Normal[#1]&/@{(Sqrt[m ω ℏ])p1,m ω ℏ p2,(Sqrt[(ℏ/(m ω))])x1,(ℏ/(m ω))x2,((((ℏ/(m ω))))^2)x4,aOp,adagOp}]];

Some algebraic checks

Let us start checking in which sense operator algebra is realized. Each power of x or p at most connect states differing by 1 in n, this means that in a computation involving k powers all rows and columns up to dim - k must be correct, the errors are due to lost states, with n > dim-1.

dim = 6;
fixnaturalUnits[False];

{P,P2,Q,Q2,Q4,aOp,adagOp} = XPmatrices[dim];

Commutation relations

Q.P - P.Q//Simplify//MatrixForm

We see that only last elements is wrong.

Algebraic consistency

Q.Q - Q2//Simplify//MatrixForm

{Q.Q.Q.Q- Q4//Simplify//MatrixForm, Q2.Q2 - Q4//Simplify//MatrixForm}

P.P - P2//Simplify//MatrixForm

Parity

Let us give a first example of fixed parity use of XPmatrices.

{p2Even,q2Even,q4Even}= XPmatrices[dim/2,0];

MatrixForm/@(Simplify/@{p2Even,P2})

The reader can check that matrix elements of p2Even correspond to matrix elements between even states of P2, i.e. n = 0,2,4 ... (indices 1,3,5 ... of the matrix).

The Hamiltonian

The Hamiltonian for an harmonic oscillator of frequency ω is:

H = (1/(2 m))P2 + (1/2)m (ω^2)Q2 //Simplify//MatrixForm

as expected.

States

In Schrödinger representation eigenstates of an harmonic oscillator are given by (in natural units):

ψH[n_,x_]:= (1/(π^(1/4)))(1./(Sqrt[(2^n)(n)!]))HermiteH[n,x]Exp[-x^2/2.];

We can construct approximate wave functions using this form of ψH[n,x]. To construct a generic state we have to perform a sum of the form

∑(Subscript[c, n]) ψBasis[n,x]

then it is useful to construct lists of length dim (the dimension we choose for matrices). If we consider all possible states a use of Table is sufficient, if we work, as often happens, with states of definite parity, is better to sum only states of given parity. For even parity n runs over 0,2,4,... (2 dim-2) while for odd parity runs over 1,3,5 ... (2 dim-1), If we want a general purpose function we can make a call to ψ in the following way

N.B. The choice of phases in ψH is related to the off diagonal matrix elements of Q, P etc. Let us check that the phases have been chosen consistently by comparing matrix elements to integral of x:

N@Q/.{ℏ-> 1,m-> 1,ω-> 1}//MatrixForm

Table[Chop[NIntegrate[ψH[n,x]ψH[j,x] x,{x,-∞,∞},AccuracyGoal-> 10]],{n,0,5},{j,0,5}]//MatrixForm

There is agreement, the same is true for P

Chop[N@P/.{ℏ-> 1,m-> 1,ω-> 1}]//MatrixForm

Table[Chop[NIntegrate[ψH[n,x](1/I)D[ψH[j,x],x] ,{x,-∞,∞},AccuracyGoal-> 10]],{n,0,5},{j,0,5}]//MatrixForm

A simple numerical application

It is convenient to use natural units and construct bigger matrices.

dim = 50;
fixnaturalUnits[True];

Harmonic oscillator

It must be clear to the reader that the choice of the base is completely arbitrary. In particular with a basis for an oscillator of frequency ω (ω = 1) we can compute approximately eigenvalues, matrix elements etc. of any system. Let us consider in particular an harmonic oscillator of frequency β · ω and compute its energy eigenvalues.

β = 0.7;

Consider only even states for simplicity:

{P2,Q2,Q4}= XPmatrices[dim,0];

Hβ = (1/2)P2 + ((β^2)/2)Q2;

approximateE = Sort[Eigenvalues[Hβ]]

We found dim eigenvalues by diagonalizing a dim×dim matrix. We expect only lower states to be reasonably approximated (higher states are more mixed with states n > dim ). For the first 10 (even) levels the approximate and expected eigenvalues are

approximateE[[Range[1,10]]]

expectedE = (β/2)+ β*Range[0,18,2]

Numerical applications

Anharmonic oscillator

Let us consider the system, in natural units:

H = ((p^2)/2)+ ((x^2)/2)+ (g/2)(x^4)

H describes an anharmonic oscillator, the coupling g express the anharmonicity. It is easy to found eigenvalues with the previous method. Let us take the extreme case g = 10 (very strong anharmonicity):

g = 10.;

H = P2/2 + Q2/2 + g/2 Q4;
ene = Sort[Eigenvalues[H]][[Range[5]]]

A quite precise determination of the energy of the ground state will be given in a subsequent notebook

ene0precise = 1.224587036059193459134396953093865213110138999308673549845104365885027406655096801688754290649704430;

The discrepancy is

(ene[[1]]-ene0precise)/ene0precise

We see that in this case we have an exceptional good agreement with a little effort.

Double well

An interesting example is given by the so called double well

g=0.1`;
V[x_]:=-((x^2)/2)+((g (x^4))/2);
Plot[V[x],{x,-1.5` (Sqrt[(1/g)]),1.5` (Sqrt[(1/g)])},AxesLabel->{"x","V[x]"}]

We can uses the matrices already computed and found the energy eigenvalues. It is interesting to find also eigenvectors, this can be done using the Mathematica function Eigensystem. The function Ordering[list] gives the order corresponding to Sort[list] and must be used to connect the order of eigenvalues to the sequence of eigenvectors.

H = P2/2 - Q2/2 + g/2 Q4;

{ene,eigV} = Eigensystem[H];
ord = Ordering[ene];
{ene,eigV} = {ene[[ord]],eigV[[ord]]};

Using the function ψH for the basis of harmonic oscillator it is easy to write the eigenfunction for the ground state of the system

ψGroundState[x_]:= ψH[Range[0,2dim-2,2],x].eigV[[1]]

We plot together the wave function and the potential, to emphasize the symmetry of the state. The sign is fixed positive by computing the function at x=0.2.

sign=Sign[ψGroundState[0.2]];

Show[Plot[V[x],{x,-1.5` (Sqrt[(1/g)]),1.5` (Sqrt[(1/g)])},AxesLabel->{"x","V[x]"},PlotStyle->RGBColor[1,0,0],PlotRange->{{-4,4},{-1.5`,1.2`}}],Plot[sign*ψGroundState[x],{x,-((1.5`)/(Sqrt[g])),((1.5`)/(Sqrt[g]))}]]

We leave to the reader the analogous computation for the first excited state (remember to call XPmatrices with odd parity).

Other examples

Mathematica computes easily power and exponential of a matrix, in this way other more exotic examples can be studied.

Morse potential

The Morse potential is defined by

V[x_]:= A (Exp[-2α x] - 2 Exp[-α x ])

Plot[V[x]/.{A->2,α->1},{x,-1,5},AxesLabel->{"x","V[x]"}]

The known eigenvalues are

(Subscript[E, n])=-A ( 1 - ((α ℏ)/(Sqrt[2 m A]))(((n + (1/2)))^2));   n+(1/2)< ((Sqrt[2 m A])/(α ℏ))

eneMorse = Table[-A ((( 1 - ((α ℏ)/(Sqrt[2 m A]))(n + (1/2))))^2)/.{A-> 2,α-> 1},{n,0,2}]

We will choose A=2;α = 1, we have only two bound states, from the limitation on n:

A=2; α = 1;
eneMorse =N@Table[-A ((( 1 - ((α ℏ)/(Sqrt[2 m A]))(n + (1/2))))^2),{n,0,((Sqrt[2  A])/(α ))-(1/2),1}]

Using the function MatrixExp we can construct an approximate Hamiltonian.

We need all the H.O. states (the potential has no definite parity). We use N[] to convert everything to machine numbers

dim = 50;
{P,P2,Q,Q2,Q4,aOp,adagOp} = N[XPmatrices[dim]];

H = (1/2)P2 + 2*(MatrixExp[-2 Q] - 2 MatrixExp[- Q]);

ene=Sort[Eigenvalues[H]][[Range[5]]]

We see that we obtain the correct number of bound states (negative energy) and a quite good approximation for the eigenvalues.

(eneMorse-ene[[Range[2]]])/eneMorse

Cosh-squared potential

As a second example let us consider

V[x_]:= -A/((Cosh[α x])^2);

Plot[V[x]/.{A->2,α->1},{x,-4,4},AxesLabel->{"x","V[x]"}]

Bound states have energies

(Subscript[E, n])= - (((ℏ^2)(α^2))/(8m))((( - (1+ 2n) + (Sqrt[1+ ((8 m A)/((α^2)(ℏ^2)))])))^2);  for  n < (1/2)(-1 + (Sqrt[1+ ((8 m A)/((α^2)(ℏ^2)))]))

We choose A=2, α = 1, then

nMax=N[(1/2)(-1 + (Sqrt[1+ ((8 m A)/((α^2)(ℏ^2)))]))/.{A-> 2,α-> 1}]

We expect two bound states with energy:

eneCh = Table[- (((ℏ^2)(α^2))/(8m))((( - (1+ 2n) + (Sqrt[1+ ((8 m A)/((α^2)(ℏ^2)))])))^2)/.{A-> 2,α-> 1},{n,0,nMax}]//N

Using the same matrices as before(we could do better using  fixed parity matrices but for full matrices we obtain both bound states)

H = (1/2)P2 - 2*4*MatrixPower[Inverse[MatrixExp[ Q] + MatrixExp[ -Q]],2] ;

ene=Sort[Eigenvalues[H]][[Range[5]]]

Again the agreement is very good.

Exercises

The reader can explore different bases in Hilbert space. A very useful and efficient one is given by Fourier series: the system can be put in a large box and free particle energy eigenstates can be taken as a basis.

A potential of the following form

Plot[-Tanh[α (x+L)]+Tanh[α (x-L)]/.{L->5,α->10},{x,-10,10},PlotRange->All]

can be taken as an analytical version of a finite wall. Study with the methods of this Notebook the bound states and explore the variations in the spectrum due to a small perturbation.