6. Quantum Mechanics

Legal Notice, Disclaimer, and Book Ordering Information:

Null

Setup

\(Clear["Global`*"]; 
Off[General::spell]; 
Off[General::spell1]; 
    \)

Needs["Calculus`VectorAnalysis`"]

Needs["Graphics`ParametricPlot3D`"]

<< "D:\Mma\Inits\ch6init.m"

Clear["Global`*"];

6.2 Problems

6.2.3  Three Dimensional Schrödinger Equation

Problem 4:  Particle with negative energy in a finite spherical box.

Consider a three-dimensional square well of infinite depth and radius a defined by V[r] = -V0 for r<a and V[r] = 0 for r>a..  Assume the wave function separates and is of the form ψ[r,θ,ϕ] = R[r] Θ[θ] Φ[ϕ] .  The radial wave function satisfies the equation  ((k r)^2 - ℓ(ℓ+1)) R[r] + 2r R'[r] + r^2 R''[r] == 0 , where En = (ℏ k)^2/(2 m) + V .

Part (a)   Use the conditions that the wave function must be finite at the origin and vanish as r->∞.  Write the general form of the wave function inside and outside the well.

rEq = ((k r)^2 - ℓ(ℓ+1)) R[r] + 2r R'[r] + r^2 R''[r] == 0//ExpandAll

rRule = DSolve[rEq,{R[r]},r]//Flatten//PowerExpand

kRule = {k1 -> Sqrt[2 m (V0 - Wn)]/ℏ
		,k2 -> Sqrt[2 m (   - Wn)]/ℏ}//PowerExpand

wave1[k_,el_,r_] = A BesselJ[(1/2)+el, k r]/Sqrt[r];

symBes[k_,el_,r_] :=
	B (BesselJ[(1/2)+el, k r] + I BesselY[(1/2)+el, k r])/Sqrt[r]

asymBes[k_,el_,r_] :=
	B (BesselJ[(1/2)+el, k r] - I BesselY[(1/2)+el, k r])/Sqrt[r]

((symBes[I q,2,r])//TrigToExp)//ExpandAll//Factor

((asymBes[I q,2,r])//TrigToExp)//ExpandAll//Factor

wave2[k_,el_,r_] = symBes[k,el,r]

Part (b)  Use the boundary condition at the wall to derive an expression for the eigenvalues.

bc1 =(wave1[k1,ℓ,r] == wave2[k2,ℓ,r])//.{r->a};

bc2 =(D[wave1[k1,ℓ,r],r] == D[wave2[k2,ℓ,r],r])//.{r->a};

bc1Sol = Solve[bc1,B]//Flatten

bc2Sol = (bc2//.bc1Sol/.A->1)//Simplify

Part (c)  Find the energy eigenvalues for ℓ=0 and plot the radial solution.

allRules = Union[kRule,{ℏ->1,a->1,m->1,V0->100}];

eq3 = bc2Sol//.{ℓ->0}//PowerExpand//Simplify

Plot[eq3[[1]]//.allRules//Evaluate,{Wn,0,99}];

eigen[0] = {33,62,83,95}//Map[FindRoot[eq3[[1]]//.allRules//Evaluate,{Wn,#}]&,#]&

eq3[[1]]//.allRules//.eigen[0]//Chop[#,10^-6]&

region1s = (wave1[k1,0,r]/.A->1//.allRules//.eigen[0]);

region2s = (wave2[k2,0,r]//.bc1Sol/.ℓ->0/.A->1//.allRules//.eigen[0]);

plotOp[region_,range_List] :=
	Plot[region//Evaluate,range
		,PlotStyle->{Dashing[{.05,.00}]
					,Dashing[{.04,.02}]
					,Dashing[{.03,.03}]
					,Dashing[{.02,.04}]	}
		,DisplayFunction->Identity
	]

Show[ plotOp[region1s,{r,0.001,1}]
	 ,plotOp[region2s,{r,1,1.3}]
	 ,Graphics[Line[{{1,-1},{1,3}}]]
	 ,Graphics[
	 	Text[FontForm["Boundary",{"Courier-Bold",10}],{1,2}]]
	 ,PlotRange->{-0.6,2.8}
	 ,DisplayFunction->$DisplayFunction
	];

Part (d)  Find the energy eigenvalues for ℓ=1 and plot the radial solution.

Part (c)  Find the energy eigenvalues for ℓ=0 and plot the radial solution.

eq4 = bc2Sol//.{ℓ->1}//PowerExpand//Simplify

Plot[eq4[[1]]//.allRules//Evaluate,{Wn,0,99}];

eigen[1] = {17,49,74,91}//Map[FindRoot[eq4[[1]]//.allRules//Evaluate,{Wn,#}]&,#]&//Chop[#,10^-6]&

eq4[[1]]//.allRules//.eigen[1]//Chop[#,10^-6]&

region1p = (wave1[k1,1,r]/.A->1//.allRules//.eigen[1]);

region2p = (wave2[k2,1,r]//.bc1Sol/.ℓ->1/.A->1//.allRules//.eigen[1]);

Show[ plotOp[region1p,{r,0.001,1}]
	 ,plotOp[region2p,{r,1,1.3}]
	 ,Graphics[Line[{{1,-1},{1,3}}]]
	 ,Graphics[
	 	Text[FontForm["Boundary",{"Courier-Bold",10}],{1,2}]]
	 ,PlotRange->{-0.5,1.3}
	 ,DisplayFunction->$DisplayFunction
	];