Mathematica on Simple harmonics
Masatsugu Sei Suzuki and Itsuko S. Suzuki
Department of Physics, SUNY at Binghamton
Binghamton, New York 13902-6000, U.S.A.
(Date: January 10, 2023)
We use Matrhematica Vesion 13.1.
_______________________________________________
Mathematica program-1
_______________________________________________

Creation operator and annihilation operator
(OverHat[(SuperPlus[a])])  (CR; creation operator)
 (OverHat[a]) (AN: annihilation operator)

Clear["Global`*"];

AN:=(1/(Sqrt[2]))(ξ # +D[#,ξ])&;CR:=(1/(Sqrt[2]))(ξ # -D[#,ξ])&;

Ground state wave function

f21=AN[ψ1[ξ]]==0//Simplify

f22=DSolve[f21,ψ1[ξ], ξ]

ψ1[0,ξ_]:=ψ1[ξ]/.f22[[1]]//.{(C[1])->(π^(-1/4))}

ψ1[0,ξ]

Wave function ψ2[n,ξ] in the real ξ  space (in general case)

ψ1[n_,ξ_]:=(1/(Sqrt[n!]))Nest[CR, ψ1[0,ξ],n]//FullSimplify

Table[{n,ψ1[n,ξ]},{n,0,10}]//TableForm[#,TableHeadings->{None,{"n","ψ[n,ξ]"}}]&

wave function ψ2[n,ξ] with the use of Hermite polynomial

ψ2[n_,ξ_]:=(1/(Sqrt[(Sqrt[π])  n! (2^n)]))HermiteH[n, ξ] Exp[((-(ξ^2))/2)];

Table[{n,ψ2[n,ξ]}//Simplify,{n,0,10}]//TableForm

Equivalence between ψ1[n,ξ] and ψ2[n,ξ]    (proof)

Table[{n,(ψ1[n,ξ]-ψ2[n,ξ])}//Simplify,{n,0,10}]//TableForm

___________________________________________________
Mathematica program - 2
___________________________________________________

Simple Harmonics wave function

Plot of probability ; P[n,ξ]=ψ[n,ξ]^2

Clear["Global`*"];ψ[n_,ξ_]:=(2^(-n/2))  (π^(-1/4)) (((n!))^(-1/2)) Exp[-((ξ^2)/2)]  HermiteH[n,ξ];P[n_,ξ_]:=((ψ[n,ξ])^2);

pt1[n_]:=Module[{f1, f2},f1=Plot[Evaluate[P[n,ξ]],{ξ,-6,6},PlotPoints->100,PlotRange->All,PlotStyle->{Red, Thick},DisplayFunction->Identity,Frame->True]; f2=Graphics[{Text[Style["n="<>ToString[n], Black, 15], {4, 0.4}], Arrowheads[0.05],Blue, Thick, Arrow[{{0,0},{0, 0.6}}],Arrow[{{-6,0},{6, 0}}], Text[Style["ξ", Black, Italic, 15],{5,0.05}] , Text[Style["P[n,ξ]", Black, Italic, 12],{-3,0.5}]}];Show[f1, f2]]; pt2=Evaluate[Table[pt1[n],{n,0,12}]];

Show[GraphicsGrid[Partition[pt2,2]]]

___________________________________________________
Mathematica Program - 3
___________________________________________________

Formula
(((-1))^n) Exp[(ξ^2)] D[Exp[-(ξ^2)],{ξ,n}]==HermiteH[n,ξ]
Proof

Clear["Global`*"];H[n_,ξ_]:=(((-1))^n) Exp[(ξ^2)] D[Exp[-(ξ^2)],{ξ,n}]

f1=Table[{n,H[n,ξ]},{n,0,5}]//Simplify//TableForm[#,TableHeadings->{None,{"n","H[n,ξ]"}}]&

f2=Table[{n,HermiteH[n,ξ]},{n,0,5}]//TableForm[#,TableHeadings->{None,{"n","HermiteH[n,ξ]"}}]&

Table[{n,(H[n,ξ]-HermiteH[n,ξ])}//Simplify,{n,0,10}]//TableForm

___________________________________________________
Mathematica Program-4

Clear["Global`*"];

AN[ψ[ξ]]=(1/(Sqrt[2]))(ξψ[ξ]+D[ψ[ξ],ξ])
CR[ψ[ξ]]=(1/(Sqrt[2]))(ξψ[ξ]-D[ψ[ξ],ξ])
CR1[ψ[ξ]]=ξψ[ξ]-D[ψ[ξ],ξ]

AN:=(1/(Sqrt[2]))(ξ # +D[#,ξ])&;CR:=(1/(Sqrt[2]))(ξ # -D[#,ξ])&;CR1:=(ξ # -D[#,ξ])&;

Proof
J[n,ξ]=[n,(((-1))^n) Exp[(ξ^2)/2] D[ Exp[-(ξ^2)/2] ψ[ξ] ,{ξ,n}]
 
is equal to 		

HE1[n,ξ]=(((CR1))^n)ψ[ξ] = Nest[CR1,ψ[ξ],n]

J[n_,ξ_]:=(((-1))^n) Exp[(ξ^2)/2] D[ Exp[-(ξ^2)/2] ψ[ξ] ,{ξ,n}]//FullSimplify;

HE1[n_,ξ_]:=Nest[CR1,ψ[ξ],n]//FullSimplify;

Table[{n,J[n,ξ]}//Simplify,{n,0,4}]//TableForm[#, TableHeadings->{None,{"n" ,"J[n,ξ]"}}]&

Table[{n,HE1[n,ξ]}//FullSimplify,{n,0,4}]//TableForm[#, TableHeadings->{None,{"n" ,"HE1[n,ξ]"}}]&

Table[{n,(J[n,ξ]-HE1[n,ξ])}//Simplify,{n,0,10}]//TableForm[#, TableHeadings->{None,{"n" ,"J[n,ξ]-HE1[n,ξ]"}}]&

___________________________________________________
Mathematica Program-5
___________________________________________________

Simple Harmonics wave function
Fourier Transform of ψ[n,ξ]; ψ[n,κ]=Φ[n,κ]

Clear["Global`*"];conjugateRule={Complex[re_,im_]:>Complex[re,-im]};
Unprotect[SuperStar];SuperStar/:(SuperStar[exp_ ]):=exp/.conjugateRule;
Protect[SuperStar];ψ[n_,ξ_]:=(2^(-n/2))  (π^(-1/4)) (((n!))^(-1/2)) Exp[-((ξ^2)/2)]  HermiteH[n,ξ]

Φ[n_,κ_]:=FourierTransform[ψ[n,ξ],ξ,κ,FourierParameters->{0,-1}]

Table[{n,Φ[n,κ]},{n,0,5}]//TableForm[#, TableHeadings->{None, {"n","Φ[n,κ]"}}]&

χ[n_,κ_]:=(((-I))^n)(2^(-n/2))  (π^(-1/4)) (((n!))^(-1/2)) Exp[-((κ^2)/2)]  HermiteH[n,κ]

Table[{n,(Φ[n,κ]-χ[n,κ])//Simplify},{n,0,10,1}]//TableForm[#, TableHeadings->{None, {"n","χ[n,κ]"}}]&

pt1[n_]:=Module[{f1, f2},f1=Plot[Evaluate[(SuperStar[Φ[n,κ]]) Φ[n,κ]],{κ,-6,6},PlotPoints->100,PlotRange->All,PlotStyle->{Red, Thick},Frame->True];f2=Graphics[{Text[Style["|Φ[n,κ]((|)^2)", Black, Italic,12], {-4,0.5}],Text[Style["n=" <>ToString[n], Black, Italic,12], {4,0.5}],Text[Style["κ", Black, Italic,12], {5,0.04}], Blue, Thick, Arrowheads[0.05], Arrow[{{0,0},{0, 0.6}}], Arrow[{{-6,0},{6, 0}}]}];Show[f1, f2]];

pt2=Evaluate[Table[pt1[n],{n,0,8}]];Show[GraphicsGrid[Partition[pt2,1]]]

pt1[0]

pt1[1]

pt1[2]

pt1[3]

pt1[4]

pt1[5]

ξ

pt1[6]

pt1[7]

pt1[8]

pt1[9]

pt1[10]

pt1[15]

_______________________________________________
Mathematica Program - 6
_______________________________________________

Classical limit of the simple harmonics

Clear["Global`*"];ψ[n_,ξ_]:=(2^(-n/2))  (π^(-1/4)) (((n!))^(-1/2)) Exp[-((ξ^2)/2)]  HermiteH[n,ξ];plot1=Plot[Table[ψ[n,ξ]^2+n+0.5,{n,0,12}]//Evaluate,{ξ,-6,6},PlotStyle->Table[{Hue[0.07 i], Thick},{i,0,10}]];plot2=Plot[ (ξ^2)/2,{ξ,-6,6},PlotStyle->{Purple,Thickness[0.005]}];plot3=Plot[Table[n+0.5,{n,0,12}]//Evaluate,{ξ,-6,6},PlotStyle->Table[Hue[0.07 i],{i,0,10}],Prolog->AbsoluteThickness[2]];f1=Show[plot1,plot2,plot3,PlotRange->{{-6,6},{0,12}}]

f2=Show[plot1,plot2,plot3,PlotRange->{{-4,4},{0,4}}]

_______________________________________________
Mathematica Program -7
_______________________________________________

Classical limit->

Clear["Global`*"];wc[n_,ξ_]:=Which[-(Sqrt[2 n +1])⩽ξ⩽(Sqrt[2 n +1]),(1/( π ))(1/(Sqrt[2 n +1-(ξ^2)])),ξ⩾(Sqrt[2 n +1]),0,ξ⩽-(Sqrt[2 n +1]) ,0];ψ[n_,ξ_]:=(2^(-n/2))  (π^(-1/4)) (((n!))^(-1/2)) Exp[-((ξ^2)/2)]  HermiteH[n,ξ];

Class[n1_]:=Module[{f1, f2},f1=Plot[Evaluate[wc[n1,ξ]],{ξ,-12,12},PlotStyle->{Blue,Thickness[0.004]}, PlotPoints->200]; f2=Graphics[{Green, Thick, Line[{{-(Sqrt[2 n1 +1]),0},{-(Sqrt[2 n1 +1]),0.6} }],Line[{{(Sqrt[2 n1 +1]),0},{(Sqrt[2 n1 +1]),0.6} }] }]; Show[f1, f2]];Prob[n1_]:=Module[{g1,g2, g3}, g1=Plot[ψ[n1,ξ]^2,{ξ,-12,12},PlotStyle->{Red,Thickness[0.004]},PlotPoints->200,Ticks->{Range[-12,12,2], Range[0,0.40, 0.1 ]}, PlotRange->All];g2=Graphics[{Text[Style["ξ", Black, 12, Italic], {9,0.03}], Text[Style["ψ[n,ξ]^2", Black, 12, Italic], {-7,0.375}], Text[Style["n=" <>ToString[n1]], {3,0.375}], Black, Thick,Arrowheads[0.04], Arrow[{{0,0},{0,0.6}}], Arrow[{{-12,0},{12,0}}]}];g3=Show[g1, g2]];PC[n1_]:=Show[Prob[n1], Class[n1],PlotRange->{{-12,12},{0, 0.6}} ];

PC[0]

PC[1]

PC[2]

PC[3]

PC[4]

PC[5]

PC[6]

PC[7]

PC[8]

PC[9]

PC[10]

PC[11]

PC[12]

PC[15]

PC[20]

PC[30]

PC[50]

PC[60]

PC[70]

A1=(1/(π (Sqrt[2 n+1])))/.n->60//N

((φ[60,0])^2)//N

PC[65]

PC[70]

_______________________________________________
Mathematica Program - 8
_______________________________________________

Commutation relations for the creation and annihilation operators

Clear["Global`"];CR:=(1/(Sqrt[2]))(ξ # -D[#,ξ])&;AN:=(1/(Sqrt[2]))(ξ # +D[#,ξ])&;N1:=CR[AN[#]]&;N11:=(CR[AN[#]]+AN[CR[#]])&;

[a,  ((((SuperPlus[a])))^2)] = 2a+

f11=AN[CR[CR[ψ[ξ]]]]-CR[CR[AN[ψ[ξ]]]]//Simplify

f12=2CR[ψ[ξ]]//Simplify

f11-f12//Simplify

[(SuperPlus[a]), (a^2)] = - 2a

f21=CR[AN[AN[ψ[ξ]]]]-AN[AN[CR[ψ[ξ]]]]//Simplify

f22=-2AN[ψ[ξ]]//Simplify

f21-f22//Simplify

[(a^2),((((SuperPlus[a])))^2)] =2((SuperPlus[a]) a + a  (SuperPlus[a])]

f31=AN[AN[CR[CR[ψ[ξ]]]]]-CR[CR[AN[AN[ψ[ξ]]]]]//Simplify

f32=2N11[ψ[ξ]]//Simplify

f31-f32//Simplify

[a,((((SuperPlus[a])))^3)] = 3(((a+))^2)

f41=AN[CR[CR[CR[ψ[ξ]]]]]-CR[CR[CR[AN[ψ[ξ]]]]]//Simplify

f42=3CR[CR[ψ[ξ]]]//Simplify

f41-f42//Simplify

[ (a^2),((((SuperPlus[a])))^3)] = 2  ((((SuperPlus[a])))^2) a+2 a ((((SuperPlus[a])))^2) +2 (SuperPlus[a])a  (SuperPlus[a])

f51=2CR[CR[AN[ψ[ξ]]]]+2AN[CR[CR[ψ[ξ]]]]+2 CR[AN[CR[ψ[ξ]]]]//Simplify

f52=AN[AN[CR[CR[CR[ψ[ξ]]]]]]-CR[CR[CR[AN[AN[ψ[ξ]]]]]]//Simplify

f52-f51//FullSimplify

_______________________________________________
Mathematica Program - 9
_______________________________________________

Clear["Global`*"];

Derivative operator of Hamiltonian of SH

AN:=(1/(Sqrt[2]))(ξ # +D[#,ξ])&;
CR:=(1/(Sqrt[2]))(ξ # -D[#,ξ])&;H1:=ℏ ω (CR[AN[#]]+(1/2) #)&;

Schrodinger equation of SH;
(1+2 n-(ξ^2)) ψ[ξ]+ψ''[ξ])=0

g1=H1[ψ[ξ]]-ℏ ω(n+(1/2))ψ[ξ]//FullSimplify

We assume that ψ(ξ)=Exp[-((ξ^2)/2)]H0(ξ)

rule1={ψ->(Exp[-(((#)^2)/2)]  H0[#]&)};

g11=g1/.rule1//Simplify

H0(ξ) satisfies the Hermite differential equation.
H0(ξ)=HermiteH[n,ξ]

g2=2 n H0[ξ]-2 ξ H0'[ξ]+H0''[ξ];

DSolve[g2==0, H0[ξ], ξ]

Table[{n,HermiteH[n,ξ]},{n,0,10}]//TableForm[#, TableHeadings->{None,{"n" ,"HermiteH[n,ξ]"}}]&

_______________________________________________
Mathematica Program -10
_______________________________________________

Prove that 
(2ξ-(d/(dξ)))^n1=Hermite[n,ξ]

Clear["Global`*"];G1:=(2 ξ #-D[#,ξ])&;G[n_,ξ_]:=Nest[G1,1,n]//Expand;f[n_,ξ_]:=G[n+1,ξ]-2 ξ G[n,ξ]+ 2 n  G[n-1,ξ]//FullSimplify;Table[{n,f[n,ξ]},{n,1,6}]//TableForm[#, TableHeadings->{None,{"n","f(n,ξ)"}}]&

_______________________________________________
Mathematica Program -11
_______________________________________________

S[ξ,t]=Exp[-(t^2)+2t ξ]
((d^n)/(dt^n))S[ξ,t](Subscript[|, t=0])=Hermite[n,ξ]

Clear["Global`*"];D1:=D[#,t]&;S[ξ_,t_]:=Exp[-(t^2)+ 2 t ξ];G1[n_,ξ_]:=Nest[D1,S[ξ,t],n]/.t->0;Table[{n,G1[n,ξ]},{n,0,6}]//TableForm[#, TableHeadings->{None,{"n","(Subscript[H, n])(ξ)"}}]&