Chapter 2 - Quantum Theory

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

0. Physical constants or frequently used formulas

Execute the package below will give you the full list of physical constants in Mathematica.

<<PhysicalConstants`

Execute the package below helps do unit conversion.

<<Units`

1. The black body radiation origin of quantum theory

1.1 Stefan's law

Use Stefan's law, what is the total radition power output from our body? calculate for normal temperature and at fever 104 F

We need to know the body surface area (BSA), take a typical value for an adult as 1.73 (m^2). We need to know human skin emissivity, which is practically close to water, 0.97. We look up Stefan Boltzmann's constant:

σ=StefanConstant

σ=((5.6704`*^-8 Watt)/((Kelvin^4) (Meter^2)))

We define a function according to Stefan's law:

Whuman[T_,ϵ_,σ_,A_]:=A ϵ σ (T^4) ;

The total power from human body at T=97.6 F (average for external skin, not at the body inner core) is:

ϵ=0.97  ;  A=1.73 ;
T=100*(97.6-32)/180 +273.15  ;
Whuman[T Kelvin,ϵ,σ, A  Meter^2]

We can plot to see power as a function of our fever:

Plot[Whuman[(100*(T-32)/180 +273.15 ) Kelvin,ϵ,σ, A  Meter^2]/Watt , {T,97.6, 105} , PlotRange-> {{95,105},{870,930}}
,PlotStyle->{ {Thick,Red}, {Thick,Black},{Thick,Blue}}
,Frame-> True, GridLines-> Automatic,FrameLabel-> {Style["Temperature (F)",20 ],Style["Power (W)",20]}  ,LabelStyle-> {20, FontFamily-> "Arial"}
]

Is it enough to detect?

What is the total BB power emitted by the entire population of the world?

Assume 2014 world population is 7 B

Total_power=7 (10^(9 )) Whuman[T Kelvin,ϵ,σ, A  Meter^2]

Or a bit over 6 Terawatt.

People have talked about thermophotovoltaic technology that can convert body thermal radiation into energy. Suppose you can wear a garment made up of this material, which has an energy conversion efficiency of 25%. How much electrical power you can get? Is it enough to power your laptop and phone?

Assume our garment area is ~ 1 (m^2). The converted power is:

T=97.6 ;
0.25 Whuman[(100*(T-32)/180 +273.15 ) Kelvin ,ϵ,σ, 1 Meter^2]

Yes.

1.2 Rayleigh-Jeans's law

Plot the Rayleigh-Jeans's energy spectral density for different temperature

We define a function:

uRJ[λ_,kB_,T_]:=((8 π)/(λ^4))kB T  ;

kB=BoltzmannConstant

kB=((1.3806504`*^-23 Joule)/Kelvin)

λ=10^-4 cm ;
T=1000. Kelvin ;
uRJ[λ,kB,T]

What we should plot is erg/(cm^3) nm, because this is energy per unit volume per unit of wavelength. We choose CGS unit: erg/(cm^3) for energy density, and nm for unit of wavelegnth. Also, for plotting, we must get rid of all units. A way to do this is:

UnitConvert[Quantity[1.,"Joules"],"Ergs"]

kBErg=kB * 10^7  *Kelvin/Joule

Plot3D[uRJ[λ*10^-4,kBErg,T]*10^-7, {λ,0.2, 0.7} ,  {T,500, 2000} , PlotPoints-> {20,5}, PlotRange ->  All, BoxRatios-> {1,1,0.67}
,AxesLabel->{"Wavelength (μm)","Temperature (K)","BB Spectral density (erg/(cm^3)nm)"}]

What is the total energy density between (Subscript[λ, 0]) and infinity?

We perform integration:

T= .  ;   λ=.  ;
8 π (Subscript[k, B]) T(∫_(Subscript[λ, 0])^∞)(1/(λ^4))\[DifferentialD]λ

Or:                ((8 π)/(3 (Subscript[λ, 0]^3)))(Subscript[k, B]) T
which is infinite as (Subscript[λ, 0]) -> 0. This obviously doesnot make physical sense

1.3 Planck's law

1.3.1 Quantization of photon

Calculate photon energy in unit of eV as a function of wavelength from 2 um to 0.2 um

Ephoton=h(c/λ)
We need various constant and proper unit conversion.
h is usually in Joule Second or Joule/Hertz
1 Joule= 1 Coulomb Volt = 1 (Coulomb/e) e Volt
Thus, we define:

heV=PlanckConstant q/ElectronCharge

heV=((4.1356673336325144`*^-15 Joule q Second)/Coulomb)

heV = heV/.  Joule-> Coulomb Volt

We want speed of light in unit of micron/second

c=SpeedOfLight

c=((299792458 Meter)/Second)

cμ=c*10.^6/Meter

Now we have the special product:

heV cμ

Pc=heV cμ/(q Volt)

Very useful formula:  (Subscript[E, ph])=(1.2398/(λ (μm))) eV

Plot[Pc/λ ,{λ,0.2,2},PlotRange->{0,6}
,PlotStyle->{ {Thick,Red}, {Thick,Black},{Thick,Blue}}
,Frame-> True, GridLines-> Automatic,FrameLabel-> {Style["Wavelength (μm)",20 ],Style["Photon Energy (eV)",20]}  ,LabelStyle-> {20, FontFamily-> "Arial"}]

In long-distance optical communication, the wavelength of light is ~ 1.55 μm. How many photons are generated per second in a 1-mW beam? If it takes 200 photons to register one bit, and 0 photons for bit-0, how many bits can 1-mW beam transmit per second?

The 1.55 μm photon energy is

λ=1.55 ;
E155=Pc/λ

The unit is eV. To convert to Joule, we just x by electron charge

2.3*SpeedOfLight

UnitConvert[Quantity["SpeedOfLight"]]

E155*QuantityMagnitude[UnitConvert[Quantity["ElectronCharge"]]]

The number of photons/sec in one 1 mW beam is

10^-3/(E155*QuantityMagnitude[UnitConvert[Quantity["ElectronCharge"]]] )

The number of bits, assuming approx equal bits 1 and bits 0 is:

7.80288208103913`*^15/100.

1.3.2 Planck's BB theory

Plot the Planck's BB radiation energy density

We define a function:

uPBB[λ_,T_,Planck_List]:=Module[ {ν, kB, c, h} ,
{kB, c, h} =Planck; ν=c/λ ;
((8 π)/(λ^4)) ((h ν )/((E^(h ν/(kB T)))-1))   ] ;

We will use CGS unit below:

Planck={10^7 BoltzmannConstant Kelvin/Joule
, 10.^2*SpeedOfLight Second/Meter
,10^7*PlanckConstant  /(Joule Second)}

Planck={1.3806504000000002`*^-16,2.99792458`*^10,6.62606896`*^-27}

λ=10^-4  ;
T=1000.   ;
uPBB[λ,T,Planck]*10^-7

We plot in erg/(cm^3) nm,

Plot[{uPBB[λ*10^-4,300,Planck] *10^-7  
,uPBB[λ*10^-4,1300,Planck] *10^-7 
,uPBB[λ*10^-4,2300,Planck] *10^-7
,uPBB[λ*10^-4,3300,Planck] *10^-7
,uPBB[λ*10^-4,4300,Planck] *10^-7
,uPBB[λ*10^-4,5300,Planck] *10^-7
,uPBB[λ*10^-4,6300,Planck] *10^-7}  ,{λ,0.2,5},PlotRange->All,Frame->True, GridLines->Automatic, PlotStyle->Table[{Hue[0.15*i,1,1], Thick},{i,0,6}]
,FrameLabel-> {Style["Wavelength (μm)",20 ],Style["BB Ener Spectral density",20]}  
,LabelStyle-> {20, FontFamily-> "Arial"}
]

We plot on Log scale

LogPlot[{uPBB[λ*10^-4,300,Planck] *10^-7  
,uPBB[λ*10^-4,1300,Planck] *10^-7 
,uPBB[λ*10^-4,2300,Planck] *10^-7
,uPBB[λ*10^-4,3300,Planck] *10^-7
,uPBB[λ*10^-4,4300,Planck] *10^-7
,uPBB[λ*10^-4,5300,Planck] *10^-7
,uPBB[λ*10^-4,6300,Planck] *10^-7}
, {λ,0.2,10}  , PlotRange -> {10^-10,10^-1} , Frame-> True, GridLines-> Automatic, PlotStyle->Table[{Hue[0.15*i,1,1], Thick},{i,0,6}] 
,FrameLabel-> {Style["Wavelength (μm)",20 ],Style["BB Ener Spectral density",20]}  
,LabelStyle-> {20, FontFamily-> "Arial"} ]

2. Classical atom theory (homework)

2.1 An exercise on "natural" units

Below is a list of defined quantities; we define them to be so, based on the fundamental quantities. We can use them as “natural” units. The fundamental quantities are electron charge e, electron mass m, speed of light in vacuum c, Planck’s constant h. In the questions below, use the Gaussian system of unit (CGS and electrostatic unit).

2.1.1 Find the Bohr’s radius in unit of nanometer

(Subscript[a, B])=((ℏ^2)/(m (e^2)))

ℏ= PlanckConstant / (2 Pi) *10^7 Erg/Joule
e=ElectronCharge*SpeedOfLight*10 Esu Second/(Coulomb Meter)
m=  Convert[ElectronMass,Gram]

ℏ=1.054571628251774`*^-27 Erg Second ;
e=4.80320427187535`*^-10 Esu ;
m=9.109382149999999`*^-28 Gram ;

(Subscript[a, B])

(Subscript[Na, B])=(Subscript[a, B]) /. Gram->(Erg^2) (Second^2)/(cm*Esu^2)

(Subscript[Na, B]) /.cm->10^8 Angstrom

2.1.2	Find the classical electron radius in unit of nanometer. (note: this is just a convenient number and has no real physical meaning; electron appears to be a point charge down to 10^-16 m, and there is no evidence that it has any internal structure - it is also very round)

e = .; ℏ=. ;c= .; m=. ;

(Subscript[r, cr])=((e^2)/(m (c^2)))

c=SpeedOfLight*100. cm/Meter 
e=ElectronCharge*SpeedOfLight*10 Esu Second/(Coulomb Meter)
m=  Convert[ElectronMass,Gram]

c=((2.99792458`*^10 cm)/Second)

e=4.80320427187535`*^-10 Esu

m=9.109382149999999`*^-28 Gram

(Subscript[r, cr])

(Subscript[Nr, cr])=(Subscript[r, cr]) /. Esu->(Sqrt[Gram*cm^3/(Second^2)])

Note: this is ONLY a number, it is NOT truly electron radius in the physical sense. 
It is likely to be <(10^(-16)) cm. 
It is also very round:
http://phys.org/news/2011-05-electron-surprisingly-scientists-year.html

2.1.3	Find the ratio of rc/aB in terms of fundamental quantities (not a numerical value). What is this ratio in terms of the fine structure constant ?

e = .; ℏ=. ;c= .; m=. ;
(Subscript[a, B])=((ℏ^2)/(m (e^2))) ;
(Subscript[r, cr])=((e^2)/(m (c^2)))  ;

((Subscript[r, cr])/(Subscript[a, B]))

We define the fine structure constant: (Subscript[α, F])=((e^2)/(ℏ c)) ;

((e^4)/((c^2) (ℏ^2))) /. ℏ->e^2/((Subscript[α, F]) c)

Thus, the ratio of classical electron radius to the Bohr radius is (Subscript[α, F]^2)

2.1.4	What is the unit and value of the fine structure constant ? (up to 4 significant digits is good enough). Redefine some unit system (e. g. change the Coulomb definition, meter definition, Joule definition…) so that the fine structure constant can exactly be equal to 1.

2.2 Classical quantum theory of the atom (Bohr's model)

Look at the Bohr’s model for the hydrogen atom in the figure. 


The electron has a charge -e, the proton charge is e, r is the radius of the electron orbit. v is the electron tangential velocity. The electron mass is m. For simplicity, let the proton be immovable. Use the de Broglie’s wavelength to derive the energy by doing the follow:

2.2.1 Express its momentum in terms of m and v; what is its de Broglie’s wavelength λ?

Momentum

p=m*v

De Broglie's wavelength

λ=(h/p)

2.2.2 In order for the orbit to repeat itself, it must be a multiple of the de Broglie’s wavelength, i. e. , where n is a positive integer. Use this quantization condition to express the kinetic energy as a function of r

Kinetic energy

p= .;λ = .;
KE=((p^2)/(2 m))

p=.  ;
KE=KE /. p->h/λ

Orbit quantization condition

quantize=Solve[n*λ==2*Pi*r,λ]

Substitute this in KE expression

KE=KE /.quantize[[1,1]]

Substitute h by 2π ℏ

(((n^2) (ℏ^2))/(2 m (r^2)))

This the is kinetic energy as a function of r:     KE=(((n^2) (ℏ^2))/(2 m (r^2)))

2.2.3 Express its electrostatic potential energy in terms of e and r.

Electrostatic potential

PE=-e^2/r

2.2.4 Express its total energy (kinetic energy + potential energy) as a function of r.

Total energy

TE=KE+PE

2.2.5 Qualitatively sketch the total energy as a function of the orbit radius; does the function have a minimum? Find the orbit radius where the total energy is a minimum.

We will qualitative plot the total energy:  (((n^2) (ℏ^2))/(2 m (r^2)))-((e^2)/r)
We need to simplify the expression:  (((n^2) (ℏ^2))/(2 m (r^2)))-((e^2)/r), using the various natural unit above. (See 2.1). We will use the Bohr radius:

(Subscript[a, B])=  . ; 
(((n^2) (ℏ^2))/(2 m (r^2)))-((e^2)/r)  /. ℏ->(Sqrt[(e^2) m (Subscript[a, B])])

We factor out a unitless energy term:      ((e^2)/(2 (Subscript[a, B])))((( (n^2) ((Subscript[a, B])^2))/(r^2))-((2 (Subscript[a, B]))/r)) 
Now we can define a unit less variable ρ=(r/(Subscript[a, B]))

((e^2)/(2 (Subscript[a, B])))((( (n^2) ((Subscript[a, B])^2))/(r^2))-((2 (Subscript[a, B]))/r))   /. r-> ρ  (Subscript[a, B])

We now can sketch the energy as a function of ρ

TEf[ρ_,n_]:=(((n^2)/(ρ^2))-(2/ρ))  ;

Plot[Evaluate@Table[TEf[ρ,n],{n,1,3}],{ρ,0.5,20}, PlotRange-> {-1,0.1}, PlotStyle-> Thick]

To see where the minimum comes from:

kplot=Plot[1/r^2,{r,0.3,5}];
pplot=Plot[-1/r,{r,0.3,5}];
Show[pplot,kplot, PlotRange-> All, AxesOrigin-> {1,0}]

To find the minimum, we take the derivative of the total energy function

(∂_ρ)TEf[ρ,n]

Solve[(2/(ρ^2))-((2 (n^2))/(ρ^3))==0,ρ]

2.2.6 Find the quantized total energy.

The quantized total energy is  (((e^2) )/(2 (Subscript[a, B])))(((n^2)/(ρ^2))-(2/ρ)) for ρ->(n^2)

(((e^2) )/(2 (Subscript[a, B])))(((n^2)/(ρ^2))-(2/ρ)) /. ρ->(n^2)

Or we can write:   Quantized energy =  -((e^2)/(2 (n^2) (Subscript[a, B])))=-((e^2)/(2 (Subscript[a, B])))(1/(n^2))

2.2.7  Use the above definition, express the total energy find in 2.2.6 in terms of Rydberg.

Since: (Subscript[R, ∞])=((e^2)/(2 (Subscript[a, B]))) 
Quantized energy =  -((e^2)/(2 (n^2) (Subscript[a, B])))=-((Subscript[R, ∞])/(n^2))

2.2.8 Find the quantized angular momentum

Angular momentum:  L= p r
The quantized angular momentum is obtained by substitute p and r at the quantized condition

L= p r  /. p-> 2 π ℏ/λ

Substitute the relation between r and λ:   λ->((2 π r)/n)

L= L /. quantize[[1,1]]

Hence, angular momentum is quantized    L = n ℏ

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:

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

F[w_,T_]:=(((4 π sin(π T) cos(((T w)/2))-2 w cos(π T) sin(((T w)/2)))/(T (4 (π^2)-(w^2)))))^2  ;

We plot to see how the square of this function varies as a function of T. Below is DFT plot

F[w_,T_]:=(((4 π sin(π T) cos(((T w)/2))-2 w cos(π T) sin(((T w)/2)))/(T (4 (π^2)-(w^2)))))^2  ;
Manipulate[
GraphicsRow[{Plot[Sin[2 π t] ,{t,-T/2,T/2},  PlotRange->{{-10,10}, {-1,1}}],
Plot[F[w,T],{w,0,4 π }, PlotRange->{{0,4 π },{0,0.3}},
Ticks->{{0, 2 Pi, 4 Pi},Automatic} ]}]
,{T,1,20}]

We notice that as we increase the sampling time, in other words, larger Δt, then the frequency uncertainty Δν is less, and vice versa.
Below is DFT plot

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}}],
ListPlot[Sp, Joined-> True, PlotRange-> {{0,0.5*Length[Sp]},All}]}]
,{a,1,20}]

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))

Ψ[x_List,t_,k_List]:=Exp[  I 2 π ( x . k -(k . k)  t)] ;

k=1.  ;
Manipulate[Plot3D[Re[Ψ[{x,y},t,{k Cos[θ], k Sin[θ]} ] ],
{x,0, 3},{y,0,3}, PlotPoints->{30,30}, Mesh-> False],{t,0, 2 π /k^2} 
,{ θ,0., Pi/2} ]

3.3 Compare 2 free electrons

We will compare two free electrons

Ψ[x_List,t_,k_List]:=Exp[  I 2 π ( x . k -(k . k)  t)] ;

k1=1.  ; k2=1.7  ;
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}]} ]
,{t,0, 2 π /k1^2}]

Compare with light wave:

ω1=1.  ; ω2=2  ; 
k1=ω1;  k2=ω2;
Manipulate[GraphicsRow[{Plot3D[Cos[k1 x -ω1 t ],
{x,0,8 π/k2},{y,0,1}, PlotPoints->{60,3}, Mesh-> False
, BoxRatios-> {3,1,1}] , Plot3D[Cos[k2 x -ω2 t ],
{x,0,8 π/k2},{y,0,1}, PlotPoints->{40,3}, Mesh-> False
, BoxRatios-> {3,1,1}]} ]
,{t,0, 8π /ω2}]

What is the difference?

3.5 Quantum well

Below is the plot of electron wavefunctions in 3 coupled quantum wells.

4. Illustration: infinite quantum well

4.1 Eigenstates

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;

Tmax=1/Min[En]   ;   Np=100  ;
disp1={}  ;
For[i=0,i<Np,i++,{
  cycle={} ;
For[j=1,j≤ Length[k],j++,
AppendTo[cycle,
Plot[{En[[j]]+Psi[j,x]*Re[Exp[I*(-2*Pi*En[[j]]*i*Tmax/Np)]],En[[j]]},
{x,-a/2,a/2}, Filling->{1-> {2}}, FillingStyle-> {Hue[(j-1)*0.2,1,1]}]  ] ]  ;
  AppendTo[disp1, cycle]
}
]

Manipulate[Show[disp1[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{-1,10}}, ImageSize->{100,250}
,Frame->True, AspectRatio-> 2.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp1],1} ]

4.2 Non steady-state 1+2

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
disp2={}  ;
For[i=0,i<Np,i++,
AppendTo[disp2, 
Plot[{Re[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ],
Im[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ]}
,{x,-a/2,a/2},Filling-> {1-> {2}}, FillingStyle->{Hue[0.7,1,1] } 
, PlotRange-> {{-a/2,a/2},{-2.5,2.5}} , ImageSize->{400,400}]
]

 ]

Manipulate[Show[disp2[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{-3,10}}, ImageSize->{100,150}
,Frame->True, AspectRatio-> 1.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp2],1} ]

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
disp3={}  ;
For[i=0,i<Np,i++,
AppendTo[disp3, Plot[Abs[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ]^2,
{x,-a/2,a/2}, Filling-> Axis, FillingStyle->{Hue[0.9,0.3,0.5] } 
, PlotRange-> {{-a/2,a/2},{0,8}} , ImageSize->{400,400}]
]

 ]

Manipulate[Show[disp3[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{-3,10}}, ImageSize->{100,150}
,Frame->True, AspectRatio-> 1.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp3],1} ]

4.3 Non steady-state 1+3

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
disp4={}  ;
For[i=0,i<Np,i++,

AppendTo[disp4,Plot[Re[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[3,x]*(E^(-I 2 π  En[[3]]*i*Tmax/Np)) ],
{x,-a/2,a/2}, Filling-> Axis, FillingStyle->{Hue[0.4,1,1] } 
, PlotRange-> {{-a/2,a/2},{-3,3}} , ImageSize->{400,400}]
]

 ]

Manipulate[Show[disp4[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{-3,10}}, ImageSize->{100,150}
,Frame->True, AspectRatio-> 1.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp4],1} ]

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
disp5={}  ;
For[i=0,i<Np,i++,
AppendTo[disp5, 
Plot[Abs[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) - Psi[3,x]*(E^(-I 2 π  En[[3]]*i*Tmax/Np)) ]^2,
{x,-a/2,a/2}, Filling-> Axis, FillingStyle->{Hue[0.05,1,1] } 
, PlotRange-> {{-a/2,a/2},{0,8}} , ImageSize->{400,400}]
]

 ]

Manipulate[Show[disp5[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{-3,10}}, ImageSize->{100,150}
,Frame->True, AspectRatio-> 1.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp5],1} ]

4.4 Non steady-state 1+2+3

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
disp6={}  ;
For[i=0,i<Np,i++,
AppendTo[disp6, Plot[Abs[2 Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  ))
 - Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np  )) 
+ Psi[3,x]*(E^(-I 2 π  En[[3]]*i*Tmax/Np)) ]^2,
{x,-a/2,a/2}, Filling-> Axis, FillingStyle-> {Hue[0.1,1,1] } 
, PlotRange-> {{-a/2,a/2},{0,25}} , ImageSize->{400,400}]
 ]
]

Manipulate[Show[disp6[[i]], DisplayFunction-> $DisplayFunction,
PlotRange->{{-a/2,a/2},{0,25}}, ImageSize->{100,150}
,Frame->True, AspectRatio-> 1.5, AxesOrigin-> {0,0}] 
,{i,1, Length[disp6],1} ]

4.5 Re and Im

a=1.  ;
k=Pi*Range[4]/a  ;  
En=k^2 /20. ;
Psi[n_,x_]:=Sqrt[2 k[[n]]/(a k[[n]] +Cos[(n-1) Pi]*Sin[a k[[n]] ])] *Cos[k[[n]] x+(n-1)*Pi/2]   ;
Tmax=1/Min[En]   ;   Np=100  ;
dispRI={}  ;
For[i=0,i<Np,i++, {alline={}  ;
For[j=1,j≤ 50,j++, xj=-a/2+(j-1)*a/49.  ;
AppendTo[alline,
Line[{{xj, 0,0}, 
{xj, 
 Re[Psi[1,xj]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,xj]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ],
Im[Psi[1,xj]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,xj]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ]}}
            ]  ]  ];
Psicurve=ParametricPlot3D[{x, 
 Re[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ],
Im[Psi[1,x]*(E^(-I 2 π  En[[1]]*i*Tmax/Np  )) + Psi[2,x]*(E^(-I 2 π  En[[2]]*i*Tmax/Np)) ]}
, {x,-a/2,a/2} , DisplayFunction-> Identity]  ;
   AppendTo[dispRI,                
Show[Psicurve, Graphics3D[alline], PlotRange->{{-a/2,a/2},{-3,3},{-3,3}}
, BoxRatios-> {5,3,3}, ImageSize->{400,400}
, DisplayFunction -> $DisplayFunction ]
]  }

 ]

Manipulate[Show[dispRI[[i]]] 
,{i,1, Length[dispRI],1} ]

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++,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++,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,
Show[allc1[[i]], DisplayFunction-> $DisplayFunction] 
,Show[allc1[[2*Np-i]], DisplayFunction-> $DisplayFunction] ]]

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++,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++,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++,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++,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++,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}] ]]