HW #4

1. Hard Sphere Scattering

For the hard sphere scattering, the requirement on the radial wave functionn is (Subscript[R, l])(
          r) = (Subscript[j, l])(k r) cos (Subscript[δ, l]) + (Subscript[n, l])(k r) 
            sin (Subscript[δ, l])for r > a, and (Subscript[R, l])(a) = 0.  Therefore, tan (Subscript[δ, l]) = -(((Subscript[j, l])(
                k a))/((Subscript[n, l])(k a))).  Using the trigonometric identity \(TraditionalForm\`\(sin\^2\) θ = \(\(\(tan\^2\) θ\)\/\(1 + \(tan\^2\) θ\) = 1\/\(1 + \(cos\^2\) θ\)\)\), and for the case k a  =  100,

sin2delta = 
        Table[{l, 
            N[((((BesselJ[(((2  l + 1))/2), z]))^2)/((((BesselJ[(((2  l + 1))/2), z]))^2) + (((BesselY[(((2  l + 1))/2), z]))^2))) /. {z -> 
                    100}]}, {l, 0, 200}];

ListPlot[sin2delta]

It appears that (sin^2) (Subscript[δ, l])behaves more or less randomly between 0 and 1, and hence (1/2) on average.

Sum[
      4  π (((2  l + 1))/(100^2)) sin2delta[[l + 1, 2]], {l, 0, 
        200}]

We can analytically estimate the total cross section in the k a ≫ 1 limit, by regarding (sin^2) (Subscript[δ, l]) = (1/2) on average for 0 ≤ l ≤ k a.  Then \(TraditionalForm\`σ = \(∑\+\(l = 0\)\%\(k\ a\)\(\(4 \( π(2  l + 1)\)\)\/k\^2\) 
              1\/2 = \(\(\(2  π\)\/k\^2\) \((k\ a + 1)\)\^2 = 
            2  π\ a\^2\)\)\) up to corrections of order (((k a))^(-1)).  Indeed, the total cross section above for k a = 100 is quite close to

N[2  π (a^2) /. {a -> 1}]

2. Spherical Well Potential for l = 0. (Threshold bound states)

(a)

Bound state problem is given by r (Subscript[R, 0])(r) = 
        sin (Sqrt[(K^2) - (κ^2)]) r for r < a, and r (Subscript[R, 0])(r) = (e^(-κ  r)) for r > a.  The continuity of the logarithmic derivatives requires (Sqrt[(K^2) - (κ^2)]) 
          cot (Sqrt[(K^2) - (κ^2)]) a = -κ  or  (κ/(Sqrt[(K^2) - (κ^2)])) 
            tan (Sqrt[ (K^2) - (κ^2)]) a + 1 = 0.   To understand where the solutions appear, we plot the function z cot z

Plot[π x Cot[π x], {x, 0, 5}, 
      PlotRange -> {-2, 2}]

Note that the straight vertical lines are an artefact of the plotting by Mathematica.

For a bound state to form, we need (Sqrt[(K^2) - (κ^2)]) 
        cot (Sqrt[ (K^2) - (κ^2)]) a to be the same as -κ < 0.  Therefore, the solutions arise when (n + (1/2)) π < (Sqrt[(K^2) - (κ^2)]) 
          a < (n + 1) π, or \(TraditionalForm\`\@\(K\^2 - \((\((n + 1)\)\ π)\)\^2\) < κ < \@\(K\^2 - \((\((n + 1\/2)\) π)\)\^2\)\).
For the threshold bound states, κ = 0 and hence K cot K a = 0, and therefore K a = (n + (1/2)) π.  New solutions appear there, and the binding energies increase as K is increased.

The rest is not required, but this is how one can plot the binding energies as a function of the potential depth.

sol1 = 
      Table[{K, κ /. FindRoot[(Sqrt[(K^2) - (κ^2)]) 
                      Cot[(Sqrt[(K^2) - (κ^2)]) 
                          a] + κ /. {a -> 
                      1}, {κ, (Sqrt[(K^2) - ((((1/2) π))^2)]), (Sqrt[Max])[(K^2) - (π^2), 0], (Sqrt[(K^2) - ((((1/2) π))^2)])}][[1]]}, {K, (1/2) π, 
          5  π, (1/10) π}]

plot1 = ListPlot[sol1, PlotJoined -> True]

sol2 = 
      Table[{K, κ /. FindRoot[(Sqrt[(K^2) - (κ^2)]) 
                      Cot[(Sqrt[(K^2) - (κ^2)]) 
                          a] + κ /. {a -> 
                      1}, {κ, (Sqrt[(K^2) - ((((3/2) π))^2)]), (Sqrt[Max])[(K^2) - (((2  π))^2), 0], (Sqrt[(K^2) - ((((3/2) π))^2)])}][[1]]}, {K, (3/2) π, 5  π, (1/10) π}]

plot2 = ListPlot[sol2, PlotJoined -> True]

sol3 = 
      Table[{K, κ /. FindRoot[(Sqrt[(K^2) - (κ^2)]) 
                      Cot[(Sqrt[(K^2) - (κ^2)]) 
                          a] + κ /. {a -> 
                      1}, {κ, (Sqrt[(K^2) - ((((5/2) π))^2)]), (Sqrt[Max])[(K^2) - (((3  π))^2), 0], (Sqrt[(K^2) - ((((5/2) π))^2)])}][[1]]}, {K, (5/2) π, 5  π, (1/10) π}]

plot3 = ListPlot[sol3, PlotJoined -> True]

sol4 = 
      Table[{K, κ /. FindRoot[(Sqrt[(K^2) - (κ^2)]) 
                      Cot[(Sqrt[(K^2) - (κ^2)]) 
                          a] + κ /. {a -> 
                      1}, {κ, (Sqrt[(K^2) - ((((7/2) π))^2)]), (Sqrt[Max])[(K^2) - (((4  π))^2), 0], (Sqrt[(K^2) - ((((7/2) π))^2)])}][[1]]}, {K, (7/2) π, 5  π, (1/10) π}]

plot4 = ListPlot[sol4, PlotJoined -> True]

sol5 = 
      Table[{K, κ /. FindRoot[(Sqrt[(K^2) - (κ^2)]) 
                      Cot[(Sqrt[(K^2) - (κ^2)]) 
                          a] + κ /. {a -> 
                      1}, {κ, (Sqrt[(K^2) - ((((9/2) π))^2)]), (Sqrt[Max])[(K^2) - (((5  π))^2), 0], (Sqrt[(K^2) - ((((9/2) π))^2)])}][[1]]}, {K, (9/2) π, 5  π, (1/10) π}]

plot5 = ListPlot[sol5, PlotJoined -> True]

Show[plot1, plot2, plot3, plot4, plot5]

This plot shows the relationship between K on the horizontal axis and κ on the vertical axis.  New bound states indeed appear at K = (n + (1/2)) (π/a).  To plot the bound state energies E = -(ℏ^2) (κ^2)/2  m as a function of (Subscript[V, 0]) = 2  m (K^2)/(ℏ^2), we need to just square them in the unit of 2  m/(ℏ^2) = 1.

Dimensions[sol1]

E1 = 
      ListPlot[Table[{sol1[[i, 1]]^2, -sol1[[i, 2]]^2}, {i, 1, 
            46}], PlotJoined -> True, 
        PlotRange -> {{0, (((5  π))^2)}, {-(((5  π))^2), 
              0}}]

Dimensions[sol2]

E2 = 
      ListPlot[Table[{sol2[[i, 1]]^2, -sol2[[i, 2]]^2}, {i, 1, 
            36}], PlotJoined -> True, 
        PlotRange -> {{0, (((5  π))^2)}, {-(((5  π))^2), 
              0}}]

Dimensions[sol3]

E3 = 
      ListPlot[Table[{sol3[[i, 1]]^2, -sol3[[i, 2]]^2}, {i, 1, 
            26}], PlotJoined -> True, 
        PlotRange -> {{0, (((5  π))^2)}, {-(((5  π))^2), 
              0}}]

Dimensions[sol4]

E4 = 
      ListPlot[Table[{sol4[[i, 1]]^2, -sol4[[i, 2]]^2}, {i, 1, 
            16}], PlotJoined -> True, 
        PlotRange -> {{0, (((5  π))^2)}, {-(((5  π))^2), 
              0}}]

Dimensions[sol5]

E5 = 
      ListPlot[Table[{sol5[[i, 1]]^2, -sol5[[i, 2]]^2}, {i, 1, 
            6}], PlotJoined -> True, 
        PlotRange -> {{0, (((5  π))^2)}, {-(((5  π))^2), 
              0}}]

Show[E1, E2, E3, E4, E5]

This plot shows how new bound states appear at each (Subscript[V, 0]) = (1/(2  
                m (a^2))) (((ℏ^2)((n + (1/2)) π))^2) as the depth of the potential is increased.

(b)

Phase shift is given by matching r (Subscript[R, 0])(r) = 
        sin (Sqrt[(K^2) + (k^2)]) r for r < a and r (Subscript[R, 0])(r) = sin k r cos (Subscript[δ, 0]) + 
            cos k r sin (Subscript[δ, 0]) = sin(k r + (Subscript[δ, 0])).  The logarithmic derivatives are (Sqrt[(K^2) + (k^2)]) cot (Sqrt[ (K^2) + (k^2)]) a = 
        k cot(k a - (Subscript[δ, 0])), or (k/(Sqrt[(K^2) + (k^2)])) tan (Sqrt[ (K^2) + (k^2)]) 
          a = tan(k a + (Subscript[δ, 0])).  therefore, (Subscript[δ, 0]) = -k a + 
          arctan  (k/(Sqrt[(K^2) + (k^2)])) tan (Sqrt[ (K^2) + (k^2)]) a.  

We can rewrite it as
 (k/(Sqrt[(K^2) + (k^2)])) tan (Sqrt[ (K^2) + (k^2)]) 
          a = tan(k a + (Subscript[δ, 0]))
\(TraditionalForm\`\(\(=\)\(\(-i\) \(e\^\(i(k\ a + δ\_0)\) - e\^\(-\(i(k\ a + δ\_0)\)\)\)\/\(e\^\(i(k\ a + δ\_0)\) + e\^\(-\(i(k\ a + δ\_0)\)\)\) = \(-i\) \(\(e\^\(i\ k\ a\)\) e\^\(2  i\ δ\_0\) - e\^\(\(-i\)\ k\ a\)\)\/\(\(e\^\(i\ k\ a\)\) e\^\(2  i\ δ\_0\) + e\^\(\(-i\)\ k\ a\)\)\)\)\)
 we find
 i (k/(Sqrt[(K^2) + (k^2)])) 
          tan (Sqrt[ (K^2) + (k^2)]) a((e^(i k a)) 
                (e^(2  i (Subscript[δ, 0]))) + 
              (e^(-i k a))) = (e^(i k a)) 
            (e^(2  i (Subscript[δ, 0]))) - (e^(-i k a))
 and
 (e^(2 i (Subscript[δ, 0]))) (e^(i k a))(
            i (k/(Sqrt[(K^2) + (k^2)])) tan (Sqrt[ (K^2) + (k^2)]) a  - 
              1) = (-i (k/(Sqrt[(K^2) + (k^2)])) 
              tan (Sqrt[ (K^2) + (k^2)]) a  - 1)
 and hence
 \(TraditionalForm\`e\^\(2\ i\ δ\_0\) = \(e\^\(\(-i\)\ k\ a\)\) \(1 + \(\(i\)\(\ \)\(k\/\@\(K\^2 + k\^2\)\) \(tan\) \(\@\( K\^2 + k\^2\)\) \(a\)\(\ \)\)\)\/\(1 - \(\(i\)\(\ \)\(k\/\@\(K\^2 + k\^2\)\) \(tan\) \(\@\( K\^2 + k\^2\)\) \(a\)\(\ \)\)\)\)

(c)

The S-matrix element (e^(2  i (Subscript[δ, 0]))) has poles when
 1 - 
          i (k/(Sqrt[(K^2) + (k^2)])) tan (Sqrt[ (K^2) + (k^2)]) 
            a = 0 .  
 On the upper half plane, we can write k = i κ, where Re(κ)>0.  The poles appear where
 1 +  (κ/(Sqrt[(K^2) - (κ^2)])) 
            tan (Sqrt[ (K^2) - (κ^2)]) a = 0.  This is precisely the condition for the bound states studied in part (a).

(Subscript[δ, 0]) = 
      ArcTan[(k/(Sqrt[(k^2) + (K^2)])) Tan[(Sqrt[(k^2) + (K^2)]) a]] - 
        k a

Far away from a threshold bound state, for instance for K = (π/4),

Plot[((4  π)/(k^2)) Sin[(Subscript[δ, 0])]^2 /. {a -> 1, 
          K -> (π/4)}, {k, 0, 3}]

which is smaller than the geometric cross section.  On the other hand, exactly on the threshold resonance,

Plot[((4  π)/(k^2)) Sin[(Subscript[δ, 0])]^2 /. {a -> 1, 
          K -> (π/2)}, {k, 0, 1}]

and the cross section diverges at k = 0.  For K just above the threshold bound state, namely when the bound state exists just below E = 0,

Plot[((4  π)/(k^2)) Sin[(Subscript[δ, 0])]^2 /. {a -> 1, 
          K -> (π/2) + 0.01}, {k, 0, 0.1}]

The cross section is finite, but is much bigger than the geometric cross section.  Similarly, for K just below the threshold bound state, namely when the bound state has just disappeared above E = 0,

Plot[((4  π)/(k^2)) Sin[(Subscript[δ, 0])]^2 /. {a -> 1, 
          K -> (π/2) - 0.01}, {k, 0, 0.1}]

It still exhibits a cross section much larger than the geometric cross section.

3. Spherical Well Potential for l = 1. (Resonances)

(a)

Here is (Subscript[j, 1])(z):

PowerExpand[(Sqrt[(π/(2  z))]) BesselJ[(3/2), z]]

Here is (Subscript[n, 1])(z):

PowerExpand[(Sqrt[(π/(2  z))]) BesselY[(3/2), z]]

We use r (Subscript[R, 1])(r)

Rin = -Cos[(Sqrt[(k^2) + (K^2)]) r] + 
        Sin[(Sqrt[(k^2) + (K^2)]) r]/((Sqrt[(k^2) + (K^2)]) r)

Rout = (-Cos[k r] + Sin[k r]/(k r)) 
          Cos[(Subscript[δ, 1])] + (((-Cos[k r])/(k r)) - Sin[k r]) 
          Sin[(Subscript[δ, 1])]

Simplify[Rout]

Simplify[D[Rin, r]/Rin]

D[Rout, r]/Rout

sol = 
      Solve[-((Sqrt[(k^2) + (K^2)]) r Cos[(Sqrt[(k^2) + (K^2)]) r] + (-1 + (k^2) (r^2) + 
                          (K^2) (r^2)) Sin[(Sqrt[(k^2) + (K^2)]) r])/(r ((Sqrt[(k^2) + (K^2)]) r Cos[(Sqrt[(k^2) + (K^2)]) r] - 
                        Sin[(Sqrt[(k^2) + (K^2)]) r])) == (cot (Cos[k r]/r + k Sin[k r] - 
                        Sin[k r]/(k (r^2))) + (-k Cos[k r] + 
                      Cos[k r]/(k (r^2)) + 
                      Sin[k r]/r) )/(cot (-Cos[k r] + 
                        Sin[k r]/(k r)) + (-Cos[
                              k r]/(k r) - 
                      Sin[k r]) ) /. {r -> a}, cot]

Plot[(1/(1 + (cot^2))) /. 
            sol[[1]] /. {a -> 1} /. {k -> 0.1}, {K, 0, 10}, 
      PlotRange -> {0, 1}, PlotPoints -> 500]

Plot[(1/(1 + (cot^2))) /. 
            sol[[1]] /. {a -> 1} /. {k -> 0.1}, {K, 3.1, 
        3.2}, PlotRange -> {0, 1}]

Plot[(1/(1 + (cot^2))) /. 
            sol[[1]] /. {a -> 1} /. {k -> 0.1}, {K, 3.13, 
        3.14}, PlotRange -> {0, 1}]

firstpeak = 
      FindRoot[cot /. sol[[1]] /. {a -> 1} /. {k -> 
              0.1}, {K, 3.137}]

The corresponding value of (Subscript[V, 0]) at the peak is

(((ℏ^2) (K^2))/(2  m (a^2))) /. firstpeak

(b)

The cross section is dominated by the l = 1 partial wave around the peak.

Plot[((4  π)/(k^2)) 3  (1/(1 + (cot^2))) /. 
            sol[[1]] /. {a -> 1} /. firstpeak, {k, 0, 0.2}, 
      PlotRange -> {0, 4000}]

It shows a prominent peak at k a = 0.1, and the cross section is much larger than the geometric cross section 4  π (a^2).

(c)

We need a precise value for the phase shift at the peak

sol[[1]] /. {a -> 1, k -> 0.1} /. 
      firstpeak

(Subscript[δ, 1]) = ArcTan[(1/cot) /. %]

Rout /. {(Subscript[δ, 1]) -> -1.5707963267946148`} /. 
          firstpeak /. {a -> 1, k -> 0.1} /. {r -> 
          1}

Rin /. {(Subscript[δ, 1]) -> -1.5707963267946148`} /. 
          firstpeak /. {a -> 1, k -> 0.1} /. {r -> 
          1}

\(Plot[\(If[
            r > 1, \(Rout /. 
                sol[\([1]\)]\) /. {δ\_1 -> \(-1.5707963267946148`\)}, \(10.049875069427088`\/1.0010011782275658`\) Rin] /. 
          firstpeak\) /. {a -> 1, k -> 0.1}, {r, 0, 10}]\)

It appears more-or-less a bound-state solution at small r, because the classical turning point is ((l(l + 1))/(r^2)) = (k^2) and hence r = ((Sqrt[2])/k) ≃ 14 in our case.  However, it does oscillate far out in radius as expected for a non-bound state beyond the classical turning point,

\(Plot[\(If[
            r > 1, \(Rout /. 
                sol[\([1]\)]\) /. {δ\_1 -> \(-1.5707963267946148`\)}, \(10.049875069427088`\/1.0010011782275658`\) Rin] /. 
          firstpeak\) /. {a -> 1, k -> 0.1}, {r, 0, 100}, 
      PlotRange -> {\(-2\), 11}, PlotPoints -> 500]\)