Eigenvalues of a Spherical Square Well

In this notebook, we compute the eigenvalues En of the potential with V = -V0 for r<R and V=0 for r>R.  The derivation can be found in any graduate and many undergraduate quantum mechanics texts.

l=0

The l=0 case is quite simple.  Inside the well, the radial function has the form u(r) = A Sin[kp*r], with kp given by

kp:=((Sqrt[2 mu (V0+En)])/hbar)

while outside the well the radial function has the form u(r) = B Exp[-kappa*r], with kappa given by

kappa:=((Sqrt[-2 mu En])/hbar)

In these definitions, V0 > 0 and En < 0, mu is the reduced mass and hbar is Planck's constant divided by 2π.   The eigenvalue condition is found by matching the radial wave function and its derivative (or, equivalently, its logarithmic derivative) at r = R.  This implies that the function f(En):

f[En_]:=(kappa/kp)+(1/(Tan[kp R]))

is zero at eigenvalues En.  We can plot this function to see graphically where the eigenvalues lie and then use FindRoot to get numerical values.

A specific case

Set the mass and hbar to simple values:

hbar=1; mu=1;

Choose a particular square well:

R=1; V0=20;

Now plot f[En] in the acceptable range for bound states (-V0 < En < 0):

Plot[f[En],{En,-V0,0}]

We can see  bound state eigenvalues, near -16 and -6.  (Don't confuse the straight lines Mathematica draws to connect -∞ and +∞ with places where f[En] has a zero.)  Look up FindRoot in the Help Browser to see the possible options.  Here we apply Findroot with a guess AND bounds on where the root is expected (this is one reason to ALWAYS make a plot first):

E1 = En /. FindRoot[f[En],{En,-16,-20,-15}]

E2=En /. FindRoot[f[En],{En,-6,-15,0}]

If we want a lot of digits, we could try using N[expression, digits]:

N[E1,16]

but that doesn't work.  We can verify that these quantities have 16 digits of precision:

Precision[E1]

Instead of N, use NumberForm:

NumberForm[{E1,E2},16]

or ScientificForm:

ScientificForm[{E1,E2},16]