Chapter 2 - Quantum Theory
Pt. 2

ECE 4339
Han Q. Le (c) U. of Houston

3. Modern quantum theory

3.1 Illustration of uncertainty principle

Let's consider a signal      Cos[2 π t]
We sample the signal from -T/2 to T/2. The Fourier transform of the signal is:

FullSimplify[(1/T)(∫_(-T/2)^(T/2))(E^(I 2 π f  t )) Cos[2 π  t]\[DifferentialD]t]

F[f_,T_]:=(((Cos[f π T] Sin[π T]-f Cos[π T] Sin[f π T])/(π T-(f^2) π T)))^2  ;

Source code: Fourier transfrom of a finite harmonic signal

Block[{F,xdat,y,Sp},
Manipulate[
If[INIT== 1,
F[f_,T_]:=(((Cos[f π T] Sin[π T]-f Cos[π T] Sin[f π T])/(π T-(f^2) π T)))^2  ;
INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]];
Manipulate[
If[RUN== 1,

Column[{Manipulate[
GraphicsRow[{Plot[Sin[2 π t] ,{t,-T/2,T/2},  PlotRange->{{-10,10}, {-1,1}}
,PlotStyle-> {Thick,Red}],
Plot[F[f,T],{f,0,2 }, PlotRange->{{0,2},{0,0.3}},
Ticks->{{0, 1, 2},Automatic} 
,PlotStyle-> {Thick,Blue}]}]
,{{T,5},1,20}]
,
Manipulate[xdat= Table[{x,Sin[2 π x]},{x,0, a,0.05} ];
y=Transpose[xdat][[2]]  ;
Sp=Take[Abs[Fourier[y]]^2,Ceiling[Length[xdat]/2] ];
GraphicsRow[{ListPlot[xdat, Joined-> True, PlotRange->{{0,20}, {-1,1}}
,PlotStyle-> {Thick,Red}],
ListPlot[Sp, Joined-> True, PlotRange-> {{0,0.5*Length[Sp]},All}
,PlotStyle-> {Thick,Blue}]}]
,{a,1,20}]} ]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
,Button["FOURIER TRANSFORM OF A FINITE HARMONIC SIGNAL",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]
]

App Demo: Fourier transfrom of a finite harmonic signal

3.2 Example of a free electron

We see that a free electron behaves like a plane wave!
                       Ψ[r,t]=(e^(i (k . r -ω t)))
 where:                             ω=(1/ℏ)(((ℏ^2) (k^2) )/(2 m))

Source code: WAVEFUNCTION OF A FREE ELECTRON

Block[{Ψ},
Manipulate[
If[INIT== 1,
Ψ[x_,t_,k_]:=Exp[  I 2 π ( x . k -(k . k)  t)] ;
INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]];
Manipulate[
If[RUN== 1,
Manipulate[
Manipulate[ 
Plot3D[Re[Ψ[{x,y},t,{k Cos[θ], k Sin[θ]} ] ],
{x,0, 3},{y,0,3}, PlotPoints->{30,30}, Mesh-> False]
,
Control[{ {t,0.,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0., 2 π /k^2}]
]
,
Control[{ {k,1.,Button["momentum",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.85,0.6,1,1]]},0.5, 2.}]
,Control[{{θ,0.,Button["angle (rad)",BaseStyle->{16, Bold,Black,FontFamily-> "Calibri"},Background->Hue[0.25,1,1,1]]},0.,Pi/2.}]
 ]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
,Button["WAVEFUNCTION OF A FREE ELECTRON",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]
]

App Demo: WAVEFUNCTION OF A FREE ELECTRON

3.3 Compare 2 free electrons

We will compare two free electrons

Source code: WAVEFUNCTION OF 2 FREE ELECTRONS

Block[{Ψ},
Manipulate[
If[INIT== 1,
Ψ[x_List,t_,k_List]:=Exp[  I 2 π ( x . k -(k . k)  t)] ;

INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]];
Manipulate[
If[RUN== 1,
Manipulate[
Manipulate[GraphicsRow[{Plot3D[Re[Ψ[{x,y},t,{k1 , 0. } ] ],
{x,0,3},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}] , Plot3D[Re[Ψ[{x,y},t,{k2 , 0. } ] ],
{x,0,3},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}]} ]
,Control[{ {t,0.,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0., 2 π /k1^2}]
] 
,Control[{ {k1,1.,Button["momentum",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.85,0.6,1,1]]},0.5, 2.}]
,Control[{ {k2,1.5,Button["momentum",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.6,0.6,1,1]]},0.5, 2.}]
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
,Button["WAVEFUNCTIONS OF 2 FREE ELECTRONS",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]
]

App Demo: Wavefunctions of 2 electrons

3.4 Compare electron wave with lightwave

Source code: Comparison of electron wave vs. lightwave

Block[{Ψ, ω1, ω2},
Manipulate[
If[INIT== 1,
Ψ[x_List,t_,k_List]:=Exp[  I 2 π ( x . k -(k . k)  t)] ;

INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]];
Manipulate[
If[RUN== 1,
Manipulate[
ω1=kp1 ; ω2=kp2 ;
Manipulate[

GraphicsGrid[{{Plot3D[Re[Ψ[{x,y},t,{k1 , 0. } ] ],
{x,0,3},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}] 
, Plot3D[Re[Ψ[{x,y},t,{k2 , 0. } ] ],
{x,0,3},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}]}
,{Plot3D[Cos[3 kp1 ( x -t) ],
{x,0,2. π},{y,0,1}, PlotPoints->{60,3}, Mesh-> False
, BoxRatios-> {3,1,1}] , Plot3D[Cos[3 kp2 (x - t) ],
{x,0,2. π},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}]} }]
,Control[{ {t,0.,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0., 2 π/kp1}]
] 
,Grid[{{Control[{ {k1,1.,Button["e-momentum 1",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.85,0.6,1,1]]},0.5, 2.}]
,Control[{ {k2,1.5,Button["e-momentum 2",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.6,0.6,1,1]]},0.5, 2.}]}
,{Control[{ {kp1,1.,Button["light k 1",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.85,0.6,1,1]]},0.5, 2.}]
,Control[{ {kp2,1.5,Button["light k 2",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.6,0.6,1,1]]},0.5, 2.}]}
,{TextCell["Time scale not the same for electron and photon"
,Black,17,Italic,FontFamily->"Calibri",Background->Hue[0.15,0.8,1]]}
}]
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]
]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
,Button["ELECTRON WAVE vs. LIGHTWAVE",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]
]

App Demo: Comparison of electron wave vs. lightwave

What is the difference?

4. Illustration: Bounded electron in an infinite quantum well

Source code: Electron in an infinite quantum well

Block[{k,a, En,Psi, Tmax},
Manipulate[
If[INIT== 1,

Psi[n_,a_,x_,t_]:=En[[n]]+
Exp[-2. I π En[[n]] t ]*Sqrt[2. k[[n]]/(a k[[n]] +Cos[(n-1.) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2.] ;

INIT=0 
,TextCell["Not ready"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]];
Manipulate[
If[RUN== 1,
Manipulate[
k=Pi*Range[4]/a ; En=k^2 /20. ;Tmax=1/Min[En] ;
Manipulate[
Show[Table[
Plot[Re[Psi[n,a,x,t]] ,{x,-a/2,a/2}
,PlotRange->{{-1.1,1.1},{-1,15}}
,PlotStyle-> {Thick,Hue[n*0.2,1,0.8]}
,Filling-> {1-> {En[[n]],Hue[n*0.2,1,1,0.5]}}
,Frame->True, AxesOrigin-> {0,0}
,AspectRatio-> 1.6] 
,{n,1,4} ]
,Graphics[{Hue[0.5,1,1,0.4]
,Rectangle[{-1.1,0.},{-0.5 a,15.}]
,Rectangle[{1.1,0.},{0.5 a,15.}]
,Rectangle[{-1.1,-1.},{1.1,0.}]
,Hue[0.6,1,1,1],Thickness[0.01]
,Line[{{-0.5 a,15.},{-0.5 a,0.},{0.5 a, 0.},{0.5 a,15.}}]}

] 
, ImageSize->{250,400}
]
,Control[{ {t,0.,Button["time",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Yellow]},0.,Tmax}]]

,
Control[{ {a,1.,Button["QW width",BaseStyle->{16, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.85,0.6,1,1]]},0.75, 2.}]
]
,TextCell["Not running"
,Black,20,Italic,FontFamily->"Calibri",Background->GrayLevel[0.5]]
]
,Row[{Button["Toggle box to start, and turn off when not running",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Black],Control[{{RUN,0,"->"},{0,1}}]
}]
]
,Button["ELECTRON WAVE IN AN INFINITE QUANTUM WELL",BaseStyle->{24,Black, Bold,FontFamily-> "Calibri"}
,Background->Hue[0.4,1,1]]
,Row[{Button["Click once to initialize",BaseStyle->{16,White, Bold,FontFamily-> "Calibri"}
,Background->Red],Control[{{INIT,0,"->"},{0,1}}]
}]
]
]

App Demo: Electron in an infinite quantum well

5. Hydrogen atom

5. 1 Spherical harmonics

Y[l_,m_,θ_,ϕ_]:=Sqrt[(l-m)!(2 l+1)/((l+m)!*4 π)] (E^(I m ϕ)) LegendreP[l, m, Cos[θ]] ;

YP[l_,m_,Cθ_,ϕ_]:=Sqrt[(l-m)!(2 l+1)/((l+m)!*4 π)] (E^(I m ϕ)) LegendreP[l, m, Cθ] ;

ParametricPlot3D[((\[LeftBracketingBar]Y(0,0,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π},Boxed->False, Axes-> False, PlotRange-> All,FaceGrids->None
, PlotPoints->{20,20}]

y10=ParametricPlot3D[((\[LeftBracketingBar]Y(1,0,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π} (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{20,20}, ViewPoint-> {1.,1.,1}]

y11=ParametricPlot3D[((\[LeftBracketingBar]Y(1,1,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π} (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{20,20}, ViewPoint-> {1.,1.,1}]

y22=ParametricPlot3D[((\[LeftBracketingBar]Y(2,2,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π} (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{50,20}]

l=2  ; m=1   ;
y21=ParametricPlot3D[((\[LeftBracketingBar]Y(l,m,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π} (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{50,20}]

l=2  ; m=0   ;
y20=ParametricPlot3D[((\[LeftBracketingBar]Y(l,m,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π} (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{50,20}]

l=3  ; m=0   ;
y50=ParametricPlot3D[((\[LeftBracketingBar]Y(l,m,θ,ϕ)\[RightBracketingBar])^2) {sin(θ) cos(ϕ),sin(θ) sin(ϕ),cos(θ)},{θ,0,π},{ϕ,0,2 π},  PlotRange-> All  (*  , Boxed->False, Axes-> False,  PlotRange-> All,FaceGrids->None   *)
, PlotPoints->{100,20}]

5.2 Radial function

Rad[n_,l_,r_,Z_]:=LaguerreL[n-l-1,2 l+1,2 Z r/n]  (E^(-Z r/n )) (Z r/n)^ l   ;
RadNorm={} ;
For[n=1, n≤ 10, n++,{lnorm={} ;
 For[l=0,l<n, l++,
AppendTo[lnorm,1/ Sqrt[
       Integrate[Rad[n,l,r,1]^2 r^2,{r,0,Infinity}]]
                         ]]    ;
AppendTo[RadNorm,lnorm]  }

]  ;
Radial[n_,l_,r_,Z_]:=RadNorm[[ n,l+1 ]]*LaguerreL[n-l-1,2 l+1,2 Z r/n]  (E^(-Z r/n )) (Z r/n)^ l   ;

n=1 ; l=0  ; Z=1  ;
Plot[Radial[n,l,r,Z]^2 r^2,{r,0,5}, PlotRange-> All]

n=2; l=0  ; Z=1  ;
Plot[{ Radial[n,0,r,Z]^2 r^2
,Radial[n,1,r,Z]^2 r^2}, {r,0,20}, PlotRange-> All
,PlotStyle-> { Thick, Hue[0.5*0/n] ,Thick,Hue[0.5*1/n] ,Hue[0.5*2/n] ,Hue[0.5*3/n] ,Hue[0.5*4/n] }]

n=3  ; Z=1  ;
Plot[{Radial[n,0,r,Z]^2 r^2
,Radial[n,1,r,Z]^2 r^2
,Radial[n,2,r,Z]^2 r^2} ,{r,0,40}, PlotRange-> All
,PlotStyle-> {{Hue[0.5*0/n] ,Thickness[0.005]}
, {Hue[0.5*1/n] ,Thickness[0.005]}
 ,{Hue[0.5*2/n]  ,Thickness[0.005]}
,{Hue[0.5*3/n]  ,Thickness[0.005]}
,{Hue[0.5*4/n] ,Thickness[0.005]} }]

5.3 Wavefunction contour

rho[x_,y_,z_]:=Sqrt[x^2+y^2+z^2]   ;
Cθ[x_,y_,z_]:=z/Sqrt[x^2+y^2+z^2]  ;
Φ[x_,y_]:=ArcCos[x+I y]  ;

5.3.1  n=1 ; l=0 (S state)

n=1;  l=0  ; m=0 ; Brange=4   ;
allc1={} ;
For[i=1,i≤ 10 , i++,
Print["i  =", i] ; 
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange},{y,-Brange,Brange},{z,-Brange,Brange}
,Contours->{0.002*(11-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -Brange,Brange}}
]  ]] ;
Np=Length[allc1]  ;

Manipulate[
Show[allc1[[i]], DisplayFunction-> $DisplayFunction] 
, {i,1,Np,1} ]

5.3.2 n=2 ; l=0 (S state)

n=2;  l=0  ; m=0 ; Brange=7.   ;
allc2={} ;
For[i=1,i≤ 10 , i++,
Print["i  =", i] ; 
AppendTo[allc2,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange},{y,-Brange,Brange},{z,-Brange,Brange}
,Contours->{0.0005*(11-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -Brange,Brange}}
 ]  ]] ;
Np=Length[allc2]  ;

Manipulate[
Show[allc2[[i]], DisplayFunction-> $DisplayFunction] 
, {i,1,Np,1} ]

5.3.3 n=2 ; l=1 (P state)

m=0

n=2;  l=1  ; m=0 ; Brange=7.   ;
allc3={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc3,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange},{y,-Brange,Brange},{z,-2 Brange,2 Brange}
,Contours->{0.001*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -2 Brange,2Brange}}
]  ]] ;
Np=Length[allc3]  ;

Manipulate[
Show[allc3[[i]], DisplayFunction-> $DisplayFunction] 
, {i,1,Np,1} ]

m=+ - 1

n=2;  l=1  ; m=-1 ; Brange=7.312   ;
allc1={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-2Brange,2Brange+0.1},{y,-2Brange,2Brange+0.1},{z,- 1.5 Brange,1.5 Brange}
,Contours->{0.01*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -2 Brange,2Brange},{  -2Brange,2Brange},{  -1.5 Brange,1.5 Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Print[Show[allc1[[i]]] ]
,Print[Show[allc1[[2*Np-i]]        ]]            ]             ]

5.3.4 n=3 ; l=2 (D state)

m=0

n=3;  l=2  ; m=0 ; Brange=20.   ;
allc1={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange+0.1},{y,-Brange,Brange+0.1},{z,- 1.5  Brange,1.5 Brange}
,Contours->{0.00005*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -1.5  Brange,1.5Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] ]]

m=+ - 2

n=3;  l=2  ; m=-2 ; Brange=50.   ;
allc1={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange+0.1},{y,-Brange,Brange+0.1},{z,- Brange, Brange}
,Contours->{0.0001*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -  Brange,Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction] ]]

m=+ - 1

n=3;  l=2  ; m=-1 ; Brange=40.   ;
allc1={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange+0.1},{y,-Brange,Brange+0.1},{z,- Brange, Brange}
,Contours->{0.0001*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -  Brange,Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction] ]]

n=4; l=3 (F state)

m=0

n=4;  l=3  ; m=0 ; Brange=35.   ;
allc1={} ;
For[i=1,i< 15 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange+0.1},{y,-Brange,Brange+0.1},{z,- 1.5  Brange,1.5 Brange}
,Contours->{0.000025*(15-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -1.5  Brange,1.5Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] ]]

n=4; l=0 (S state)

m=0

n=4;  l=0  ; m=0 ; Brange=25.   ;
allc1={} ;
For[i=1,i< 25 , i++,
Print["i  =", i] ;
AppendTo[allc1,
ContourPlot3D[Abs[Radial[n,l,rho[x,y,z],1]*YP[l,m,Cθ[x,y,z],Φ[x,y]]]^2 *(x^2+y^2+z^2)
,  {x,-Brange,Brange+0.1},{y,-Brange,Brange+0.1},{z,-  Brange, Brange}
,Contours->{0.0001*(25-i)}
,PlotPoints->{5,5,5}
,PlotRange->{{  -Brange,Brange},{  -Brange,Brange},{  -  Brange,Brange}}
,DisplayFunction-> Identity]  ]] ;
Np=Length[allc1]  ;

For[i=1,i≤ 2*Np-1,i++,
If[i≤ Np,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction, ViewPoint->{1,2.3,1}] ]]