Matrix mechanics and Perturbation Theory

Quantum Mechanics
 K. Konishi, G. Paffuti

Use a finite matrix formulation based on harmonic oscillator states to implement perturbation theory. This approach can give exact results for polynomial interaction potentials to arbitrary order of perturbation theory. Formulate the same approach using operators instead of numerical matrices and write a program which can compute in principle in analytic form the perturbation theory up to an arbitrary order for a quartic oscillator.

Utilities

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

Introduction

General strategy

In notebook [NB-3.3] we have introduced a simple matrix formulation of Quantum Mechanics, useful for one dimensional problems. The idea was that a basis like the one of an harmonic oscillator is a basis the Hilbert space, albeit in general the Hamiltonian is not diagonal in this basis. A truncation of the basis turn out  to be a simple approximation to several Quantum Mechanics computations.

This idea is developed in this notebook from this point of view: let us suppose we have a basis which diagonalize some Hamiltonian (Subscript[H, 0])and let us number the basis vectors as |1⟩, |2⟩ ...  . Consider a perturbation V to (Subscript[H, 0]) and suppose that V has matrix elements which vanishes between basis vectors more distant than k,  ie. ⟨i|V|j⟩ = 0 if  |j-i| > k. Then at a finite order n of perturbation theory on a state |i⟩, only states up to |i+ k n⟩ can give a contribution, i.e. we can work with finite matrices and compute analytically an arbitrary order in perturbation theory. A typical example is a (x^4)perturbation on an harmonic oscillator. As the operator change at most by 4 the principal quantum number of the states, matrices (4 n+s)×(4 n+s) suffice to compute n-th order in perturbation theory.

All this work can be accomplished just by writing sufficiently large matrices, and using exact arithmetic. The drawback is that we have an answer for each computation, i.e. a function of couplings, but we have to do again the computation if we change starting state, i.e. we have not an answer as a function of n, the starting level.

We can use Mathematica ability to perform symbolic computations to translate matrix operations in abstract relations between operators, as is usually done in Quantum Mechanics, and this will provide analytical results both in couplings and in quantum numbers.

Perturbation theory, for non degenerate levels, can be summarized in two recursive formulas (see text). Let |(Subscript[ψ, 0])⟩ the unperturbed eigenstate relative to the energy level (Subscript[E, 0]), let (Subscript[H, 0])the unperturbed Hamiltonian and V the perturbation. Define 𝒫 and 𝒬 the projectors on the state |(Subscript[ψ, 0])⟩ and its orthogonal, i.e. 𝒬 projects out on level (Subscript[E, 0]). For the n-th order of perturbation theory we have

ψ = (Subscript[ψ, 0])+ (Subscript[ψ, 1])+ (Subscript[ψ, 2]) ... ;      E = (Subscript[E, 0])+(Subscript[ϵ, 1])+(Subscript[ϵ, 2])+ ...

Where in general (Subscript[ψ, n]) and (Subscript[ϵ, n])are of order (g^n) if g is the strength of the perturbation V. The corrections are given by

(Subscript[ψ, n]) = G V (Subscript[ψ, n-1]) - (∑_(k=1)^(n-1)) (Subscript[ϵ, n-k]) G  (Subscript[ψ, k])

(Subscript[ϵ, n]) =⟨(Subscript[ψ, 0])|V|(Subscript[ψ, n-1])⟩=⟨(Subscript[ψ, 0])|V G V|(Subscript[ψ, n-2])⟩- (∑_(k=1)^(n-2)) (Subscript[ϵ, n-k-1])   ⟨(Subscript[ψ, 0])|V G |(Subscript[ψ, k]) ⟩=⟨(Subscript[ψ, 1])|V|(Subscript[ψ, n-2])⟩-(∑_(k=1)^(n-2)) (Subscript[ϵ, n-k-1])  ⟨(Subscript[ψ, 1])|(Subscript[ψ, k]) ⟩

G is the operator

G = 𝒬 (1/((Subscript[E, 0])-(Subscript[H, 0])))𝒬

Then all we have to do is to write explicitly these operators along with operators necessary to write H, V etc.

In particular the two first orders are given by

(Subscript[ϵ, 1]) = ⟨(Subscript[ψ, 0])|V|(Subscript[ψ, 0])⟩;    (Subscript[ϵ, 2]) = ⟨(Subscript[ψ, 0])|V G V|(Subscript[ψ, 0])⟩.

The model system

We study as a model the anharmonic oscillator, but the technique can be extended to other cases. In this case the building blocks to construct matrices are those given by matrix elements of p and q, or a and (SuperDagger[a]).  These elements also can be computed in Mathematica but we take them for granted.

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 the notebook  [NB-3.3] we gave a function XPmatrices which compute matrices above up to a maximum value of n, we refer the reader to that notebook for more informations. We report the calling conventions from that notebook:

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.

Eg.

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

give matrix elements of fixed parity (positive) up to state dim/2.

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

Give matrices elements for both parity states up to dim. Here we copy the code for the matrix from notebook [NB-3.3]

XPmatrices

XPmatrices[dim_,parity_Integer:-1]:=
Module[{x2array,p2array,x4array,x2,p2,x4,x1,x1array,p1,p1array}, 
Block[{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}]]];

G matrix, numbering of states and units

We construct a G matrix for the state k, simply writing 1/ ((Subscript[E, k])- H + ε) for all states and then put to zero the element (k,k) and ε. As for XPmatrices we allow an optional argument to specify the parity.

Gmatrix[k_,dim_,parity_Integer:-1]:= Module[{factor,ε}, 
Clear[nn];
If[parity == 0 || parity == 1,nn =2k-2; factor = 2;,nn=k-1;factor =1;];
mat= (1/(ℏ ω ))DiagonalMatrix[1/(nn-factor(Range[dim]-1)+ε)];
mat[[k,k]]=0; mat = mat/.ε-> 0]

N.B. We follow usual Mathematica convention for vectors, i.e. the first element is v[[1]]. The fundamental state of the oscillator is k=1, the first excited state k=2 etc.

If we use all states the connection is simply

nOscillator = kComponent - 1

If instead we use parity fixed states, it is easy to see that

kComponent = (nOscillator + 2 - parity)/2

Eg. for positive parity (parity = 0) oscillator numbers 0,2,4,6... correspond to components 1,2,3,4... etc. To simplify the life to the reader we write a function to do the conversion

fromOscillatorToVector[n_,parity_]:= (n+2-parity)/2;

We can choose also if we want or not explicit ℏ, ω, m etc. factors

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

If we call fixnaturalUnits[False] we have all right factors and usual symbols are considered positive numbers in Simplify.

Perturbation theory for anharmonic oscillator (matrices)

Direct computation

The Hamiltonian is

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

The term proportional to g is our perturbation. The known first two perturbative correction to energy of state n are (see text or do it by yourself)

((ΔE)^((1)))= g (3/2)((((ℏ/(2m ω))))^2)(1 + 2n +2 (n^2));
Δ (E^((2)))= - (1/(4ℏ ω))(g^2) ((((ℏ/(2m ω))))^4)(68 (n^3)+102 (n^2)+118 n + 42);

checkPert[n_]:= {g (3/2)((((ℏ/(2m ω))))^2)(1 + 2n +2 (n^2)), - (1/(4ℏ ω))(g^2) ((((ℏ/(2m ω))))^4)(68 (n^3)+102 (n^2)+118 n + 42)}

Just to warm up let us check these results for ground state and first excited states using simple formulas (). We use fixed parity matrices to show hot to work with them. We take dim = 20, which is much more than we need

fixnaturalUnits[False];

dim = 20;

Ground state

stateOscillator = 0;
parity = Mod[stateOscillator,2];
n = fromOscillatorToVector[stateOscillator,parity];

{p2, q2, q4} = XPmatrices[dim/2, parity];
G = Gmatrix[n,dim/2,parity];

This is the perturbation in our basis

V = (g/2) q4;

The first two orders are given by, see ()

ris1={V[[n,n]], (V.G.V)[[n,n]]}

Let us compare with known results

ris1 == checkPert[stateOscillator]

First excited state

stateOscillator = 1;
parity = Mod[stateOscillator,2];
n = fromOscillatorToVector[stateOscillator,parity];

{p2, q2, q4} = XPmatrices[dim/2, parity];
G = Gmatrix[n,dim/2,parity];

V = (g/2) q4;

ris2={V[[n,n]], (V.G.V)[[n,n]]}

Let us compare with known results

ris2 == checkPert[stateOscillator]

Both states

We can compute matrices without fixing parity and obtain both results at once

{p,p2,q,q2,q4,aOp,adagOp}=XPmatrices[dim];
V = (g/2) q4;

n=1;
G = Gmatrix[n,dim];
ris1={V[[n,n]], (V.G.V)[[n,n]]}

n=2;
G = Gmatrix[n,dim];
ris2={V[[n,n]], (V.G.V)[[n,n]]}

Recursive computation

It is much easier, for higher orders, to define a recursive procedure based on () and (). Let us suppose for example that we want the first 10 orders for E for the ground state. Dimesion 50 for matrices is large enough (it would be sufficient 21 for fixed parity, we use 50 to allow user to check other states and higher orders).

In recursive computations a standard construct is

f[x_] := f[x] = expression

in such a way that Mathematica "remembers" computed values. It is necessary, if you do another computations to remove old values with the following instructions.

Clear[i,n,p2,q2,q4];
Off[Remove::"rmnsm"];
Remove["ψ*","ε*"];

fixnaturalUnits[True];

stateOscillator = 0;
parity = Mod[stateOscillator,2];
nState = fromOscillatorToVector[stateOscillator,parity];

dim = 50;
{p2,q2,q4} = XPmatrices[dim,parity]; 
V = (g/2)q4;
Λ = Gmatrix[nState,dim,parity];

The ground state, with g = 0, is

ψ0 = Table[0,{dim}]; ψ0[[1]] = 1;

First order corrections are

ψ[1] = Λ.V.ψ0;
ε[1] = ψ0.V.ψ0;

Higher orders follow from recursion :

ψ[n_]:=ψ[n]=Λ.V.ψ[n-1] - Sum[ε[n-k]Λ.ψ[k],{k,1,n-1}];
ε[n_]:=ε[n] = ψ0.V.ψ[n-1] ;

Here we write the first 10 coefficients in the expansion of E (  (g^k)(Subscript[c, k])=ε[k] )

coeff=Table[ε[k]/g^k,{k,1,10}]

In another notebook we compute with a similar, but more efficient method, the first several orders more for the energy of the ground state, the first 20 being

ene20={(3/8),-(21/32),(333/128),-(30885/2048),(916731/8192),-(65518401/65536),(2723294673/262144),-(1030495099053/8388608),(54626982511455/33554432),-(6417007431590595/268435456),(413837985580636167/1073741824),-(116344863173284543665/17179869184),(8855406003085477228503/68719476736),-(1451836748576538293163705/549755813888),(127561682802713500067360049/2199023255552),-(191385927852560927887828084605/140737488355328),(19080610783320698048964226601511/562949953421312),-(4031194983593309788607032686292335/4503599627370496),(449820604540765836160529697491458635/18014398509481984),-(211491057584560795425148309663914344715/288230376151711744)};

Let us compare the results

coeff == ene20[[1;;10]]

Operators

General formalism

All the computations given above can be done for an arbitrary state n in a purely algebraic way by defining directly operators. Let us denote by ket[n] the n-th state vector in the basis of our oscillator. The operators q and p are defined on individual ket and extended to linear combinations via linearity, i.e. using Distribute (see Mathematica help). We write operators with script fonts, to distinguish from others objects.

Script letters as ℊ or 𝒬 are written with the Escape key:   Escape scg Escape  or Escape scQ Escape  for the two cases above.

Clear[i,n,k,j]

fixnaturalUnits[False];

𝓆[x_?NumberQ] = 0;
𝓆[x_.ket[n_]]:=x  (Sqrt[ℏ/(2 m ω)]) (Sqrt[n])ket[n-1] +x  (Sqrt[ℏ/(2 m ω)]) (Sqrt[n+1])ket[n+1];
𝒬[x_]:= Distribute[𝓆[x]];

𝓅[x_?NumberQ] = 0;
𝓅[x_.ket[n_]]:=-I x  (Sqrt[ℏ m ω/2]) (Sqrt[n])ket[n-1] +I x  (Sqrt[ℏ m ω/2]) (Sqrt[n+1])ket[n+1];
𝒫[x_]:= Distribute[𝓅[x]];

The first instruction in each operator is needed because operators act on ket, not on numbers. The name ket is just used as a tag for an expression.

In manipulating operators is more natural to use infix notation @, so 𝒬@𝒬@ket[1] means 𝒬[𝒬[ket[1]]], but is much more readable.

We need also the functional form for matrix G, i.e. Q1/(E - (Subscript[H, 0]))Q for a given state. This is simply to write as the kets are eigenstates of (Subscript[H, 0]):

ℊ[k_][x_?NumberQ] = 0;
ℊ[k_][x_.ket[j_]]:=If[j===k,0,(x/(ℏ ω(k-j)))   ket[j]]
𝒢[k_][x_]:= Expand[Distribute[ℊ[k][(x)]]];

If after a calculation we want to project on a given state ket[n] this is simply done taking the coefficient of ket[n] in the expression:

proj[n_,expr_]:= Coefficient[expr,ket[n]];

Some checks

To see how the notation works let us see a couple of trivial examples. First of all the commutations relations, on an arbitrary ket of the basis

Simplify[𝒬@𝒫@ ket[n] - 𝒫@𝒬@ ket[n]]

A projection works in this way

proj[2,3 ket[1] + 7 ket[2]]

Anharmonic oscillator

First orders

Let us proceed as in the matrix case  first looking ad lowest orders in perturbation theory.

Potential terms are easily defined

V4[x_]:=Expand[ (1/2)g 𝒬@𝒬@𝒬@𝒬[x]];
V3[x_]:=Expand[ (1/2)g3 𝒬@𝒬@𝒬[x]];

First perturbative corrections to an arbitrary state ket[n] are

Simplify[proj[n,V4@ket[n]]]

Simplify[proj[n,V4@𝒢[n]@V4@ket[n]]]

A comparison with checkPert confirms the correctness of the result.

Simplify[checkPert[n]]

All orders

Now let us construct our iterative procedure. The starting point is an arbitrary ket[n]. The first corrections are given by

Clear[i,n,p2,q2,q4];
Off[Remove::"rmnsm"];
Remove["ψ*","ε*"];

ψ[1] = 𝒢[n] @ V4@ket[n];
ε[1] = Simplify[proj[n,V4@ket[n]]];

and recursive relations

ψ[k_]:=ψ[k]= Expand[𝒢[n] @ V4@ψ[k-1] - Sum[ε[k-s]𝒢[n]@ψ[s],{s,1,k-1}]]

ε[k_]:= ε[k]=Expand[proj[n,V4@(ψ[k-1])]]

Then for example the first eight orders for an arbitrary state ar given by

(ris8 = Block[{$RecursionLimit = 1000},Simplify[Table[ε[k],{k,1,8}]]])//Timing

For the ground state in natural units

ris8/.{n-> 0, m-> 1,ℏ-> 1,g-> 1,ω-> 1}

which is identical to previous result :

ene20[[1;;8]]

Scalar products and a different recursion relation

It can be useful to have a scalar product function in this framework. It can be most simply defined by doing the usual product and substitute 1 for product of equal ket, zero in other cases:

scalarproduct[a_,b_]:= Expand[a*b]/.{ket[i_]*ket[j_]-> If[i=== j,1,0],((ket[i_])^2)-> 1}

It would be necessary to require that a and b in the function are linear combinations of ket, but we do this for granted. We assume also that coefficients are real, the trivial change to adapt the function to the general case are left to the reader. We used the operator = = = as the two expression must be equal independently of the numerical value of the index.

Let us make an example with previous computed vectors

ψ[1]

scalarproduct[ψ[1],ket[2+n]]

With a scalar product function we can also implement the last form of recursive relations (), which turn out to be a bit faster

Clear[i,n,p2,q2,q4];
Off[Remove::"rmnsm"];
Remove["ψ*","ε*"];

ψ[0] = ket[n];
ψ[1] = 𝒢[n] @ V4@ket[n];
ε[1] = Simplify[proj[n,V4@ket[n]]];

and recursive relations

ψ[k_]:=ψ[k]= Expand[𝒢[n] @ V4@ψ[k-1] - Sum[ε[k-s]𝒢[n]@ψ[s],{s,1,k-1}]]

ε[k_]:= ε[k]=
Expand[scalarproduct[ψ[1],V4@ψ[k-2]]] - Expand[Sum[ε[k-s-1]scalarproduct[ψ[1],ψ[s]],{s,1,k-2}]]

(ris8bis = Block[{$RecursionLimit = 1000},Simplify[Table[ε[k],{k,1,8}]]])//Timing

ris8 == ris8bis

Exercises