Mathematica Program
Vector analysis (spherical coordinates)

Masatsugu Sei Suzuki and Itsuko S. Suzuki
Department of Physics, SUNY at Binghamton
(Date: February 11, 2023)

((Definition))

L->L=(r×p)
LP->(Lx+iLy),	LM->(Lx-iLy)
Leq->(L^2)
Peq->(p^2)
Lx1[n]->(((Lx))^n),		Ly1[n]->(((Ly))^n),		Lz1[n]->(((Lz))^n)		
ex->ux,	ey->uy,	ez->uz
er->ur
pr=pr->linear radial momentum (quantum mechanical)
H->H1=((1/(2 μ)) Peq+V[r]); Hamiltonian
ψ=ψ(r,θ,ϕ); wave function

Clear["Global`"];ux={Sin[θ] Cos[ϕ],Cos[θ] Cos[ϕ],-Sin[ϕ]};uy={Sin[θ] Sin[ϕ],Cos[θ] Sin[ϕ],Cos[ϕ]};uz={Cos[θ],-Sin[θ],0};ur={1,0,0};x=r ur.ux;y=r ur.uy;z=r ur.uz;Lap:=Laplacian[#,{ r,θ,ϕ},"Spherical"]&;Gra:=Grad[#,{ r,θ,ϕ},"Spherical"]&;Diva:=Div[#,{r,θ, ϕ},"Spherical"]&;Curla:=Curl[#,{r,θ, ϕ},"Spherical"]&;L:=(-I ℏ( Cross[(ur  r),Gra[#]])&)//Simplify;Lx:=(ux.L[#]&)//Simplify;Ly:=(uy.L[#]&)//Simplify;Lz:=(uz.L[#]&)//Simplify;LP:=(Lx[#]+I Ly[#])&//Simplify;LM=(Lx[#]-I Ly[#])&//Simplify;Leq:=(Nest[Lx,#,2]+Nest[Ly,#,2]+Nest[Lz,#,2])&;p:=(-I ℏ  Gra[#])&;px:=(ux.p[#])&;py:=(uy.p[#])&;pz:=(uz.p[#])&;Peq:=(Nest[px,#,2]+Nest[py,#,2]+Nest[pz,#,2])&;prq:=(((-I ℏ)/2) ur .Gra[#]+ ((-I ℏ)/2) Diva[# ur])&;pr:=(ℏ/(I r)) D[r #,r]&; H1=((1/(2 μ)) Peq[#]+V[r]#) &;Lx1[n_]:=Nest[Lx, #,n]&;Ly1[n_]:=Nest[Lx, #,n]&;Lz1[n_]:=Nest[Lz, #,n]&;

(a)	Angular momentum in spherical coordinates
Lx, Ly, Lz, 
LP=Lx+iLy,	LM=Lx-iLy
Leq=(Lx^2)+(Ly^2)+(Lz^2)
ψ=ψ(r,θ,ϕ)

a1=Lx[ψ[r,θ, ϕ]]//Simplify

a2=Ly[ψ[r,θ, ϕ]]//Simplify

a3=Lz[ψ[r,θ, ϕ]]//Simplify

a4=LP[ψ[r,θ, ϕ]]//Simplify

a5=LM[ψ[r,θ, ϕ]]//Simplify

a6=Leq[ψ[r,θ, ϕ]]//Simplify

(b)	Linear momentum in spherical coordinates
(p^2)=(px^2)+(py^2)+(pz^2)
pr=prq (linear radial momentum)
(p^2)=(pr^2)+(Leq/(r^2))
Leq=(Lx^2)+(Ly^2)+(Lz^2)

b1= Peq[ψ[r,θ,ϕ]]//Simplify

b2=Nest[prq,ψ[r,θ,ϕ],2]//Simplify

b3=(1/(r^2))Leq[ψ[r,θ,ϕ]]//FullSimplify

bs23=b2+b3//Simplify

b1-bs23//Simplify

(c)	Commutation relations 
Angular momentum, coordinates, and Hamiltonian

[Lx, y ]=iℏ z

c1=Lx[y ψ[r,θ,ϕ]]-y Lx[ ψ[r,θ,ϕ]]-I ℏ z ψ[r,θ,ϕ]//FullSimplify

[Ly, z ]=iℏ x

c2=Ly[z ψ[r,θ,ϕ]]-z Ly[ ψ[r,θ,ϕ]]-I ℏ x ψ[r,θ,ϕ]//FullSimplify

[Lz, x ]=iℏ y

c3=Lz[x ψ[r,θ,ϕ]]-x Lz[ ψ[r,θ,ϕ]]-I ℏ y ψ[r,θ,ϕ]//FullSimplify

Leq=(L^2)=(Lx^2)+(Ly^2)+(Lz^2)
(r^2)Peq=(r^2)((px^2)+(py^2)+(pz^2))

Proof
(L^2)=(r×p)·(r×p)=(r^2)(p^2)-(((r·p))^2)+iℏ(r·p)

c4=Leq[ψ[r,θ,ϕ]]//FullSimplify

c51=(r^2) Peq[ψ[r,θ,ϕ]]//Simplify

c52=-(r ur).p[(r ur).p[ψ[r,θ,ϕ]]]//Simplify

c53=I ℏ (r ur).p[ψ[r,θ,ϕ]]//Simplify

cs5=c51+c52+c53//Simplify

c4-cs5//Simplify

(d)	Commutation relation of Hamiltonian and angular momentum
H=(1/(2m)) (p^2)+V(r)
[H,(Lx^2)]=0,	[H,(Ly^2)]=0,	[H,(Lz^2)]=0,	
[H,(Lx^2)+(Ly^2)+(Lz^2)]=0,
[H,Lx]=0,	[(L^2),Lx]=0,
[H,Ly]=0,	[(L^2),Ly]=0,	
[H,Lz]=0,	[(L^2),Lz]=0,

d1=H1[ψ[r,θ, ϕ]]//Simplify

d2=H1[Lx1[2][ ψ[r,θ, ϕ]]]-Lx1[2][H1[ ψ[r,θ, ϕ]]]//FullSimplify

d3=H1[Ly1[2][ ψ[r,θ, ϕ]]]-Ly1[2][H1[ ψ[r,θ, ϕ]]]//FullSimplify

d4=H1[Ly1[2][ ψ[r,θ, ϕ]]]-Ly1[2][H1[ ψ[r,θ, ϕ]]]//FullSimplify

d5=H1[Lz1[2][ ψ[r,θ, ϕ]]]-Lz1[2][H1[ ψ[r,θ, ϕ]]]//FullSimplify

d6=H1[Lz[ ψ[r,θ, ϕ]]]-Lz[H1[ ψ[r,θ, ϕ]]]//FullSimplify

d7=H1[Leq[ ψ[r,θ, ϕ]]]-Leq[H1[ ψ[r,θ, ϕ]]]//FullSimplify

d8=H1[Leq[ ψ[r,θ, ϕ]]]-Leq[H1[ ψ[r,θ, ϕ]]]//FullSimplify

d9=Leq[Lz[ ψ[r,θ, ϕ]]]-Lz[Leq[ ψ[r,θ, ϕ]]]//FullSimplify

(e)	Commutation relations of quantum mechanical radial angular momentum and power of r
pr=prq: linear radial momentum
(pr^2), 
[pr,r], [pr,(r^2)], [pr,(r^n)],
[r,(pr^2)], [r,(pr^3)],
[pr,(1/r)], [pr,(1/(r^2))], 
[(pr^2),(1/r)], [(pr^2),(1/(r^2))],

e0=prq[ψ[r, θ,ϕ]]-pr[ψ[r, θ,ϕ]]//Simplify

e1=pr[pr[ψ[r, θ,ϕ]]]//Simplify

e2=pr[r ψ[r, θ,ϕ]]-r pr[ψ[r, θ,ϕ]]//Simplify

e3=pr[(r^2) ψ[r, θ,ϕ]]-(r^2) pr[ψ[r, θ,ϕ]]//Simplify

e4=pr[(r^n) ψ[r, θ,ϕ]]-(r^n) pr[ψ[r, θ,ϕ]]//Simplify

e5=r Nest[pr, ψ[r, θ,ϕ],2]-Nest[pr,r ψ[r, θ,ϕ],2]-2 I ℏ pr[ψ[r, θ,ϕ]]//Simplify

e6=r Nest[pr, ψ[r, θ,ϕ],3]-Nest[pr,r ψ[r, θ,ϕ],3]-3 I ℏ Nest[pr,ψ[r, θ,ϕ],2]//Simplify

e7=pr[(1/r) ψ[r, θ,ϕ]]-(1/r) pr[ψ[r, θ,ϕ]]//Simplify

e8=pr[(1/(r^2)) ψ[r, θ,ϕ]]-(1/(r^2)) pr[ψ[r, θ,ϕ]]//Simplify

e9=Nest[pr,(1/r) ψ[r, θ,ϕ],2]-(1/r) Nest[pr,ψ[r, θ,ϕ],2]//Simplify

e10=Nest[pr,(1/(r^2)) ψ[r, θ,ϕ],2]-(1/(r^2)) Nest[pr,ψ[r, θ,ϕ],2]//Simplify