Lecture 2 - Quantum mechanics of identical bosons

by Leonid Pryadko

Initialization

Pauli Matrices

(Subscript[σ, x])={{0,1},{1,0}};(Subscript[σ, y])={{0,-I},{I,0}};(Subscript[σ, z])={{1,0},{0,-1}}; e=IdentityMatrix[2];

Direct product

Unprotect[NonCommutativeMultiply];

NonCommutativeMultiply[A_?MatrixQ,B_?MatrixQ]:=KroneckerProduct[A,B];

Spin-1/2 spin matrices

sx=(Subscript[σ, x])/2;sy=(Subscript[σ, y])/2; sz=(Subscript[σ, z])/2;

Identical particles: bosons

Define truncated Hilbert space for bosons

NN=7;

Define boson creation and annihilation operators

bp:=Block[{i,j},Table[If[j==i-1,Sqrt[j],0],{i,1,NN},{j,1,NN}]];

NN=5;bp//MatrixForm

bp//MatrixForm

bm:=Transpose[bp];

bm.bp//MatrixForm

bp.bm//MatrixForm

Commutator

bm.bp-bp.bm//MatrixForm

This is the coordinate operator

xx:=(bp+bm)/Sqrt[2];xx//MatrixForm

Operator for non-linear oscillator

be:=IdentityMatrix[NN];

be//MatrixForm

H:=bp.bm+0.5 be +0.1 (xx.xx.xx.xx); H//MatrixForm

NN=17;H//MatrixForm

NN=7;Eigenvalues[H]//Sort

NN=5;Eigenvalues[H]//Sort

NN=35;vals=Eigenvalues[H]//Sort

Compare with the solution of the differential equation

solution=NDSolve[{-u''[x]/2+(x^2/2+0.1 x^4) u[x]==0.559(* vals[[1]]*)u[x],u'[0]==0,u[0]==1},u,{x,0,4.5}];

Plot[u[x]/.solution,{x,0,4.5`},PlotRange->{0,1}]

Now we can introduce many-body operators easily

solution=NDSolve[{-u''[x]/2+(x^2/2+0.1 x^4) u[x]==vals[[2]]u[x],u'[0]==1,u[0]==0},u,{x,0,4.5}];

Plot[u[x]/.solution,{x,0,4.5`},PlotRange->{0,1}]

solution=NDSolve[{-u''[x]/2+(x^2/2+0.1 x^4) u[x]==vals[[3]]u[x],u'[0]==0,u[0]==1},u,{x,0,4.5}];

Plot[u[x]/.solution,{x,0,4.5`},PlotRange->{-1.2`,1}]

Two decoupled oscillators

NN=5;

HH:=(bp.bm+0.5 be)**be+be**(bp.bm+0.5 be);

vals=.;vecs=.;

bp**be

HH//MatrixForm

{vals,vecs}=Eigensystem[HH];vals

vals=Reverse[vals];vecs=Reverse[vecs];vals

N1=(bp.bm)**be;N2=be**(bp.bm);

For[i=1,i≤Length[vals],i++,Print[vecs[[i]].HH.vecs[[i]], "   ",vecs[[i]].N1.vecs[[i]], "   ",vecs[[i]].N2.vecs[[i]]]]

Two linearly coupled oscillators

HH:=(bp.bm+0.5 be)**be+be**(bp.bm+0.5 be)+0.2 xx**xx;

Eigenvalues[HH]//Sort

NN=10;Eigenvalues[HH]//Sort

Try the conventional method for the ground state

Apply[Plus,Sqrt[Eigenvalues[{{1,0.2},{0.2,1}}]//Sort]]/2

{E1,E2}=Sqrt[Eigenvalues[{{1,0.2},{0.2,1}}]]

0.5 E1+1.5E2

Bosons with strong repulsion U

Introduce operator n

bn:=bp.bm

U:=20

H:=bn+U/2 bn.(bn-be)

Eigenvalues[H]//Sort

In the limit of U->∞, only two states remain.  Many bosons with hopping -- Bose Hubbard model (possibly, Infinite-U)

H=(U/2 bn.(bn-be))**be+be**(U/2 bn.(bn-be))-bp**bm-bm**bp;

Eigenvalues[N[H]]//Sort