Spherical well

Quantum Mechanics
K. Konishi, G. Paffuti

Study numerically an attractive square well in radial coordinates, i.e.  the potential  -(Subscript[V, 0])θ[a-r]. Consider in particular the case

(Subscript[V, 0])= 300 MeV; a = 3 fm; μ=940 MeV/(c^2).

ClearAll["Global`*"];
Off[General::"spell",General::"spell1"];
Needs["PlotLegends`"]

Introduction and units

We want to study bound states for a radial potential of the form -(Subscript[V, 0]) θ[a-r], i.e. a potential well in radial coordinate r. Radial Schrödinger equation for reduced wave function is

-((ℏ^2)/(2m))(((d^2) χ)/(dr^2)) + (((ℏ^2)/(2m (r^2)))L(L+1) - (Subscript[V, 0])θ[a-r])χ = E χ

Bound states are expected only for E < 0.

It is convenient to change variable (units) with r = a x;

(((d^2) χ)/(dx^2)) + (-(1/(x^2))L(L+1) + ((2 m (a^2))/(ℏ^2))(Subscript[V, 0])θ[1-x] + ((2 m (a^2))/(ℏ^2)) E)χ = 0

The real (adimensional) parameters of the problem are

A = (Sqrt[((2 m (a^2))/(ℏ^2))(Subscript[V, 0])]);    ξ = (Sqrt[((2 m (a^2))/(ℏ^2))(E + (Subscript[V, 0]))]) ;  η = (Sqrt[- ((2 m (a^2))/(ℏ^2)) E]);

ξ is the adimensional wave number inside the well while η is the analogous quantity in the classically forbidden region (remember that we look for eigenstates with negative energy). These parameters satisfy the relation

(ξ^2)+(η^2) = (A^2).

Inside and outside the well the equation takes the form

(((d^2) χ)/(dx^2)) + (-(1/(x^2))L(L+1) + (ξ^2))χ = 0;      ξ < 1;  region  I
(((d^2) χ)/(dx^2)) + (-(1/(x^2))L(L+1) - (η^2))χ = 0;       ξ >1;   region  II

The technique of solution has been exposed in the text and the main points will be reviewed below.

Solution for L=0

The generic regular solutions in regions I and II are

χ1[x] = A Sin[ξ x];            χ2[x] = B Exp[-η x]

The continuity of logarithmic derivative,i.e.f'/f at the transition point x=1 gives, together with () the system of trascendental equations

ξ Cot[ξ] = - η;    (ξ^2)+(η^2) = (A^2)        ⇒       ξ Cot[ξ] = - (Sqrt[(A^2)-(ξ^2)])

It is clear that every solution is less than A and is obtained graphically from the intersection of the curve η = - ξ Cot[ξ], and the circle (ξ^2)+(η^2) = (A^2). This procedure has been folowed in the text. Here we will use Mathematica function FindRoot to locate the solutions. It is better to have an idea of the function before using FindRoot.

Example: A = 6

Plot[ξ Cot[ξ]+(Sqrt[(A^2)-(ξ^2)])/.A->6,{ξ,0.2`,6}]

We have clearly 2 roots, one before 3 and the second between 3 and 6. FindRoot admits as arguments a list of starting value, and produce a list of roots:

rootsξOld = ξ/.FindRoot[ξ Cot[ξ]+(Sqrt[(A^2)-(ξ^2)])/.A-> 6,{ξ,{2,5}}]

Solution for any L

Equation () in region I is identical to the (radial) equation for a free particle with momentum ξ, so we can write at once the regular solution (see also notebook[*])   χ1[x] = (Sqrt[x]) (Subscript[J, L+(1/2)])( ξ x)

χ1[L_,x_]:= ((( x))^(+1/2)) BesselJ[L+(1/2), x]

χ is the reduced function, χ = r R,  so the additional x with respect NB5-2.  Let us make a check

FullSimplify[(ξ^2)Derivative[0,2][χ1][L,ξ x] +( -(1/(x^2))L(L+1) + (ξ^2))χ1[L,ξ x]]

Let us make explicit a trivial point which can cause troubles sometimes. From the definition of derivative
(∂_x) f[ξ x] = ξ f'[ξ x] ⇒ (∂_x) f[ξ x] = Derivative[1][f][ξ x]
this is way there is a factor ξ in front of the Derivative operator in the above check.

The wave equation in region II becomes identical to the first one with x → i x so the general solution will be of the form

(((ξ x))^(+1/2)) ( a (Subscript[J, L+1/2])[i ξ x]+ b (Subscript[Y, L+1/2])[i ξ x])

From the asymptotic expansions

(Subscript[J, ν])[x]∼ (Sqrt[(2/(π x))])Cos[ x - (π/2)ν - (π/4)];                (Subscript[Y, ν])[x] ∼ (Sqrt[(2/(π x))])Sin[ x - (π/2)ν - (π/4)];

it follows easily that asymptotically  (Subscript[J, ν])[x] ± i (Subscript[Y, ν])[x]  ∼ Exp[ ± i x]. For immaginary argument only  the first has an exponential decay, and (Subscript[J, ν])[i x] + i (Subscript[Y, ν])[i x] is named (Subscript[K, ν])[x], in Mathematica BesselK[ν,x]. As a check for asymptotic behaviour

Off[Series::"esss"];
Normal[Series[BesselK[ν,x],{x,∞,0}]]
On[Series::"esss"];

So the solution in region II is of the form  χ2[ η x] with

χ2[L_,x_]:= ((( x))^(+1/2)) BesselK[L+(1/2),x]

as a check:

FullSimplify[(η^2)Derivative[0,2][χ2][L,η x] +( -(1/(x^2))L(L+1) - (η^2))χ2[L,η x]]

Continuity

The continuity of logarithmic derivative is

(Subscript[((((∂_x)χ2[L,η x])/(χ2[L,η x]))-(((∂_x)χ1[L,ξ x])/(χ1[L,ξ x]))), x=1])  = 0      i.e.      η ((χ2'[L,η])/(χ2[L,η ]))- ξ((χ1'[L,ξ])/(χ1[L,ξ ]))=0

using the relation between ξ and η, we have to find roots of the following function

diffLogDerivative[L_,ξ_,A_]:= ((η ((Derivative[0,1][ χ2][L,η])/(χ2[L,η])))/.η-> (Sqrt[(A^2)-(ξ^2)]) )- ξ ((Derivative[0,1][ χ1][L,ξ])/(χ1[L,ξ]))

As an example for A = 4 and L = 1

Plot[diffLogDerivative[1,x,4],{x,0,4}]

ξ/.FindRoot[diffLogDerivative[1,ξ,4],{ξ,2.2}]

Example A = 6

Consider again the case A = 6. Plots with L = 0, 1 .. 5 for the matching function are given below

Show[GraphicsGrid[Partition[Table[Plot[diffLogDerivative[L,x,6],{x,0,6},PlotRange->{{0,6},{-10,10}},
PlotLabel->"L = "<>ToString[L]],{L,0,5}],2]]]

It is apparent that only for L ≤ 3 there are roots. Let us check the method with the result found at the beginning of the notebook for S-waves:

{ξ/.FindRoot[diffLogDerivative[0,ξ,6],{ξ,{2,5},{0,0},{6,6}}],"old result:",
rootsξOld}

Physical parameters and the general case

The physical case proposed as a problem has paremeters given by

(Subscript[V, 0])= 300 MeV; a = 3 fm; μ=940 MeV/(c^2);    Pl= ℏ c= 200 MeV fm;
A=(Sqrt[ ((2μ (Subscript[V, 0])(a^2))/(ℏ^2))])=(Sqrt[ ((2μ  (c^2) (Subscript[V, 0])(a^2))/((ℏ^2)(c^2)))])

(Subscript[V, 0])= 300; PlanckC = 200; μNucleon = 940; aWell=3;
A = N[(Sqrt[ ((2 *μNucleon * (Subscript[V, 0]) * (aWell^2) )/(PlanckC^2))])]

Given a root ξ for the matching condition, the energy is given by

Ene[ξ_]:=  - (Subscript[V, 0]) + (((ξ^(2 )) (PlanckC^2))/(2 μNucleon *(aWell^2)))

Let us give a plot of the matching condition for different values of L

Show[GraphicsGrid[Partition[Table[Plot[diffLogDerivative[L,x,A],{x,0,A},PlotRange->{{0,A},{-10,10}},PlotLabel->"L = "<>ToString[L]],{L,0,8}],3]]]

As is apparent from the figures above and from the expression of matching condition the function diffLogDerivative is singular at the zero of the Bessel function J[L+1/2,x] and zero values are approached from below. All the roots must be less then A. The zeros of Bessel function, and the point r = 0 fix the intervals for search the roots. Last interval between last zero of BesselJ( less than A) and ξ = A must be taken into account only if the function is positive at the endpoint ξ = A.

Mathematica 6 provides directly the k zeo of Bessel functions, BesselJZero[ν,k]. Prior versions had a function BesselJZeros[ν,x] to compute these zeros.

To find all the solution of the matching conditions we select from the vector of zeros of Bessel function those which are less than A, prepend to the result r=0 and apply FindRoot to the resulting intervals. This is done in the folowing function

findMatchingPoints[L_,A_]:=Module[{zerosBessel,xEnd,xStart,xMed,tab,lastRoot},zerosBessel=Select[Table[N[BesselJZero[L+1/2,k]],{k,1,5}],#1<A&];
xEnd=Append[zerosBessel,A];xStart=Prepend[zerosBessel,0.01`];
xMed=((xEnd+xStart)/2);ε=(1/(10^6));tab=Table[ξ/.FindRoot[diffLogDerivative[L,ξ,A]==0,{ξ,xMed[[k]],xStart[[k]],xEnd[[k]]}],{k,1,Length[xStart]-1}];If[diffLogDerivative[L,A-ε,A]<0,tab,lastRoot=ξ/.FindRoot[diffLogDerivative[L,ξ,A]==0,{ξ,Last[xMed],Last[xStart],A-ε}];Append[tab,lastRoot]]]

The  If instruction in the function return the found root except when the function diffLogDerivative is negative near ξ=A, int his case a root is searched in the last interval, between the last zero of BesselJ and A.

In the present case the ξ roots are:

TableForm[tabRoots=Table[findMatchingPoints[L,A],{L,0,7}]]

and the corresponding energies

TableForm[tabEne= Ene[tabRoots]]

Let us make a graph for the results

dis[mytab_]:=Show[Table[
Show[Graphics[{
AbsoluteThickness[2],Hue[(L)/Length[mytab]],
Line[{{L-0.4,mytab[[L+1,s]]},{L+0.4,mytab[[L+1,s]]}}]}],AspectRatio-> 1/1],
{L,0,Length[mytab]-1},{s,1,Length[mytab[[L+1]]]}]];

Show[dis[tabEne],Axes-> True,PlotRange->{{-1,8},{-300,0}},AxesOrigin->{-1,0},
AxesLabel->{"L","E(MeV)"},Ticks->{Range[0,8],Automatic}]

In usual spectroscopic notation angular momentum states are denoted by s,p,d ... Let us order the various levels using this notation.

angMomentumSymbols = {"s","p","d","f","g","h","i","l","m","n"};
(levelSymbols=Table[ToString[s]<>angMomentumSymbols[[L+1]],
{L,0,7},{s,1,Length[tabEne[[L+1]]]}])//TableForm

order = Ordering[Flatten[tabEne]];
Flatten[levelSymbols][[order]]

Comparison with Hydrogen atom

The principal quantum number n is here written as nr + L;

tabEneHydrogen=Table[-(1/(2 (((L+nr))^2))),{L,0,4},{nr,1,8}];
Show[dis[tabEneHydrogen],Axes->True,PlotRange->{{-1,5},{-0.6`,0}},AxesOrigin->{-1,0},AxesLabel->{"L","E(Ry)"},Ticks->{Range[0,5],Automatic}]

Exercises

Plot wave functions for the eigenvalues found in this notebook.