Ladder operators for the anharmonic oscillator

An example of working with non-commutative variables and symbolic evaluation.

Using CenterDot[] for non-commutative multiplication and normal ordering

ClearAll@CenterDot;
SetAttributes[CenterDot,Flat];
CenterDot[x___,u_+v_,y___]:=CenterDot[x,u,y]+CenterDot[x,v,y];
CenterDot[x___,c_ u_,y___]/;FreeQ[c,a]&&FreeQ[c,ad]:=c CenterDot[x,u,y];
CenterDot[x___,a,ad,y___]:=CenterDot[x,ad,a,y]+CenterDot[x,y];

Checking our implementation of CenterDot[]

x=(1/(Sqrt[2]))(ad+a);
p=(I/(Sqrt[2]))(ad-a);

x·p-p·x//Expand

H0=((p·p)/2)+((x·x)/2)//Expand

H=((p·p)/2)-((x·x)/2)+((x·x·x·x)/4)//Expand

The explicit Fock basis

ClearAll/@{State,Costate};
State[n_]:=(1/(Sqrt[n!]))CenterDot@@Table[ad,{n}]·Ket[0];
Costate[n_]:=(1/(Sqrt[n!]))Bra[0]·CenterDot@@Table[a,{n}];

CenterDot[x___,a,Ket[0]]:=0;
CenterDot[Bra[0],ad,x___]:=0;
CenterDot[Bra[0],Ket[0]]:=1;

Checking orthonormalization and the harmonic Hamiltonian

Table[Costate[i]·State[j],{i,0,7},{j,0,7}]//MatrixForm

Table[Costate[n]·H0·State[m],{n,0,6},{m,0,6}]//MatrixForm

Table[{n,Costate[n]·H0·State[n]},{n,0,6}]

FindSequenceFunction[%,n]//FullSimplify

The anharmonic oscillator’s energy levels

Computing Hamiltonian matrix elements and energy levels

Dynamic@{n,m}

Table[Costate[n]·H·State[m]//Simplify,{n,0,6},{m,0,6}]//MatrixForm

Sort@Eigenvalues@N@%

Guessing the general form

Table[{n,Costate[n]·H·State[n]}//Simplify,{n,0,5}]

FindSequenceFunction[%,n]//FullSimplify

Table[{n,Costate[n]·H·State[n+2]//Simplify},{n,0,5}]

(Sqrt[{First@#,(((Last@#))^2)}&/@%376//FindSequenceFunction[#,n]&])//FullSimplify

Table[{n,Costate[n]·H·State[n+4]//Simplify},{n,0,5}]

(Sqrt[{First@#,(((Last@#))^2)}&/@%385//FindSequenceFunction[#,n]&])//FullSimplify

Using the general form to define (Subscript[H, nm]) for larger matrix sizes

Hnm[k_]:=DiagonalMatrix@Table[(3/16) (1+2 n (1+n)),{n,0,k}]+DiagonalMatrix[Table[(1/8) (2 n-1)(Sqrt[(1+n) (2+n)]),{n,0,k-2}],2]+DiagonalMatrix[Table[(1/8) (2 n-1)(Sqrt[(1+n) (2+n)]),{n,0,k-2}],-2]+DiagonalMatrix[Table[(1/16) (Sqrt[(1+n) (2+n) (3+n) (4+n)]),{n,0,k-4}],4]+DiagonalMatrix[Table[(1/16) (Sqrt[(1+n) (2+n) (3+n) (4+n)]),{n,0,k-4}],-4]

Hnm[6]==%259

Computing energy levels with a <1% error

En[k_]:=
Block[{a,b},
a=Take[Sort@Eigenvalues@N@Hnm@k,k];
b=Take[Sort@Eigenvalues@N@Hnm@(k+2),k];
First/@Select[Transpose@{a,((Abs[a-b])/a)},Last@#<.01&]
]

En[100]

%-En[1000][[;;36]]

A more general approach to compute (Subscript[H, nm]) for arbitrary (n,m)

CenterDot[x___,a,Ket[n_]]:=(Sqrt[n])CenterDot[x,Ket[n-1]];
CenterDot[x___,ad,Ket[n_]]:=(Sqrt[n+1])CenterDot[x,Ket[n+1]];
CenterDot[Bra[n_],Ket[m_]]:=KroneckerDelta[n,m];

Check the harmonic Hamiltonian

Bra[m]·H0·Ket[n]//FullSimplify

On to the anharmonic case

Bra[m]·H·Ket[n]//FullSimplify