球対称井戸型ポテンシャル/メモ の変更点

更新 
&katex();
* 球ベッセル関数の導出 [#q4e20ed3]

　&math(
R''+\frac{2}{\rho}R'+\left\{1-\frac{l(l+1)}{\rho^2}\right\}R=0
);

は &math(\rho\to\infty); にて &math(R''=-R); となるから、
&math(R(\rho)\propto\sin \rho); または &math(R(\rho)\propto\cos \rho); となる。
そこで、

　&math(
R(\rho)=\sum_{k=0}^\infty \frac{s_k\sin \rho+c_k\cos \rho}{\rho^k}
);

と置いて代入すれば、

　&math(
R''&=\sum_{k=0}^\infty \left[
k(k+1)\frac{s_k\sin \rho+c_k\cos \rho}{\rho^{k+2}}
-2k\frac{s_k\cos\rho-c_k\sin\rho}{\rho^{k+1}}
-\frac{s_k\sin\rho+c_k\cos\rho}{\rho^{k}}\right]\\
&=\sum_{k=2}^\infty (k-2)(k-1)\frac{s_{k-2}\sin \rho+c_{k-2}\cos \rho}{\rho^k}\\
&\phantom{=}-\sum_{k=1}^\infty 2(k-1)\frac{s_{k-1}\cos\rho-c_{k-1}\sin\rho}{\rho^k}\\
&\phantom{=}-\sum_{k=0}^\infty\frac{s_k\sin\rho+c_k\cos\rho}{\rho^{k}}\\
);

　&math(
\frac{2}{\rho}R'&=\sum_{k=0}^\infty \left[
-2k\frac{s_k\sin \rho+c_k\cos \rho}{\rho^{k+2}}
+2\frac{s_k\cos\rho-c_k\sin\rho}{\rho^{k+1}}\right]\\
&=-\sum_{k=2}^\infty 2(k-2)\frac{s_{k-2}\sin \rho+c_{k-2}\cos \rho}{\rho^k}\\
&\phantom{=}+\sum_{k=1}^\infty 2\frac{s_{k-1}\cos\rho-c_{k-1}\sin\rho}{\rho^k}\\
);

　&math(
\left\{1-\frac{l(l+1)}{\rho^2}\right\}R&=\sum_{k=0}^\infty \frac{s_k\sin \rho+c_k\cos \rho}{\rho^k}\\
&\phantom{=}-l(l+1)\sum_{k=2}^\infty \frac{s_{k-2}\sin \rho+c_{k-2}\cos \rho}{\rho^k}
);

&math(\rho^0); については自動的に満たされる。

&math(\rho^{-1}); については、

　&math(2(k-1)c_{k-1}-2c_{k-1}-l(l+1)s_{k-2}=0);

　&math(-2(k-1)s_{k-1}+2s_{k-1}-l(l+1)c_{k-2}=0);

すなわち、

　&math(c_0=\frac{l(l+1)}{-2}s_{-1}=0);

　&math(s_0=\frac{l(l+1)}{2}c_{-1}=0);

&math(k\ge 2); については、

&math(\frac{\sin\rho}{\rho^k}); の係数より、

　&math(
&(k-2)(k-1)s_{k-2}+2(k-1)c_{k-1}-\cancel{s_k}\\
&\hspace{1cm} -2(k-2)s_{k-2}-2c_{k-1}+\cancel{s_k}-l(l+1)s_{k-2}\\
&=\{(k-3)(k-2)-l(l+1)\}s_{k-2}+2(k-2)c_{k-1}\\
&=0
);

すなわち、

　&math(
2(k-2)c_{k-1}=\{l(l+1)-(k-3)(k-2)\}s_{k-2}
);　（&math(k\ge 2);）

あるいは、

　&math(
2kc_{k+1}=\{l(l+1)-(k-1)k\}s_k
);　（&math(k\ge 0);） 

&math(\frac{\cos\rho}{\rho^k}); の係数より、

　&math(
&(k-2)(k-1)c_{k-2}-2(k-1)s_{k-1}-\cancel{c_k}\\
&\hspace{1cm}-2(k-2)c_{k-2}+2s_{k-1}+\cancel{c_k}-l(l+1)c_{k-2}\\
&=\{(k-3)(k-2)-l(l+1)\}c_{k-2}-2(k-2)s_{k-1}\\
&=0
);

すなわち、

　&math(
2(k-2)s_{k-1}=\{(k-3)(k-2)-l(l+1)\}c_{k-2}
);　（&math(k\ge 2);）

あるいは、

　&math(
2ks_{k+1}=\{(k-1)k-l(l+1)\}c_k
);　（&math(k\ge 0);）

となる。

得られた２つの漸化式は &math(k=0); で意味をなさないから、
&math(s_1,c_1); は自由に選べて、&math(k\ge 1); において、

　&math(
c_{k+1}=\frac{l(l+1)-(k-1)k}{2k}s_k
);

　&math(
s_{k+1}=\frac{(k-1)k-l(l+1)}{2k}c_k
);

となる。&math(k\to\infty); にて &math(s_k\ne 0,c_k\ne 0); であれば、

　&math(c_{k+1}\sim-\frac{k}{2}s_k);

　&math(s_{k+1}\sim \frac{k}{2}c_k);

となって明らかに発散するから、この漸化式は &math(k=l+1); で打ち切られる必要がある。

&math(l=0); のとき &math(c_2=s_2=0); より、

　&math(R=\frac{1}{\rho}(s_1\sin\rho+c_1\cos\rho));

であるが、&math(c_1\ne 0); では &math(\rho=0); で発散してしまうため、&math(c_1=0); であり、

　&math(j_0(\rho)\propto\frac{\sin\rho}{\rho});

&math(l=1); のとき、&math(c_2=s_1);、&math(s_2=-c_1); より、

　&math(R=s_1\left(\frac{\sin\rho}{\rho}+\frac{\cos\rho}{\rho^2}\right)+
c_1\left(\frac{\cos\rho}{\rho}-\frac{\sin\rho}{\rho^2}\right)
);

であるが、&math(s_1\ne 0); では &math(\rho=0); で発散してしまうため、&math(s_1=0); であり、

　&math(j_1(\rho)\propto \frac{\cos\rho}{\rho}-\frac{\sin\rho}{\rho^2}
);

上記で発散する側を選んだのが球ノイマン関数になる。

* 球ベッセル関数 [#p507b042]

 LANG:mathematica
 MySphericalBesselJ[l_, x_] := 
   Nest[D[#, x]/x &, Sin[x]/x, l] x^l // FullSimplify
 Table[MySphericalBesselJ[l, x], {l, 0, 4}]

~
&math(
&\Bigg\{\frac{\sin (x)}{x},\frac{x \cos (x)-\sin (x)}{x^2},
-\frac{\left(x^2-3\right) \sin (x)+3 x \cos (x)}{x^3},\\
&\hspace{5mm}\frac{3 \left(2 x^2-5\right) \sin (x)-x \left(x^2-15\right) \cos (x)}{x^4},\frac{5 x \left(2 x^2-21\right) \cos (x)+\left(x^4-45 x^2+105\right) \sin (x)}{x^5}\Bigg\}
);

* 球ベッセル関数のグラフ [#a9909050]

 LANG:mathematica
 Plot[
  Join[
    SphericalBesselJ[{0, 1, 2, 3}, r], 
    {1/r, -1/r}
  ] // Evaluate, 
  {r, 0, 40}, PlotRange -> {-0.3, 1.05}, ImageSize -> 800, BaseStyle -> 20, 
  PlotStyle -> {Thick, Thick, Thick, Thick, {Thick, Dotted, Gray}, {Thick, Dotted, Gray}}, 
 ]
  
 Plot[Join[
   SphericalBesselJ[{1, 5, 9}, (Pi r)], {1/(Pi r), -1/(Pi r)}] // 
   Evaluate, {r, 0, 16}, ImageSize -> 800, BaseStyle -> 20, 
   PlotStyle -> {Thick, Thick, 
   Thick, {Thick, Dotted, Gray}, {Thick, Dotted, Gray}}, 
   PlotRange -> {-0.3, 0.5}, AspectRatio -> 0.4]
 
 Plot[
  r^2 SphericalBesselJ[{0, 1, 2, 3}, r]^2 // Evaluate, 
  {r, 0, 40}, PlotRange -> Full, ImageSize -> 800, 
  BaseStyle -> 20, PlotStyle -> Thick, Filling->Axis]

* エネルギー [#s488a51e]

 LANG:mathematica
 RootsOfSphericalBesselJ[l_, xmax_] := 
   Map[ Round[#[[1]][[2]], 0.00001]&, 
     Table[
       FindRoot[SphericalBesselJ[l, x], {x, s}], 
       {s, 1, xmax, 0.1}]] // Sort // Union // 
         Select[#, Function[x, 0.1 <= x <= xmax]] &
 
 energies = 
   Table[ 
     MapIndexed[{#1^2, #2[[1]], l}&, 
       RootsOfSphericalBesselJ[l, 40]], 
    {l, 0, 10}] // Flatten[#, 1] & // 
       Sort[#, (#1[[1]] < #2[[1]]) &] &
 
 ListPlot[{#[[3]], #[[1]]} & /@ energies, PlotStyle -> PointSize[Large], 
 PlotRange -> {{-0.2, 10.2}, {0, 800}}, 
 AxesLabel -> {l, "(\!\(\*SuperscriptBox[SubscriptBox[\(\[Rho]\), \(n\)], \(l\)]\)\!\ \(\*SuperscriptBox[\()\), \(2\)]\)"}, LabelStyle -> 16, 
 GridLines -> {{}, Range[0, 800, 50]}]

* 境界条件 [#vcd00b8a]

 LANG:mathematica
 RootsOfSphericalBesselJ[l_, xmax_] :=
   Map[Round[#[[1]][[2]], 0.00001] &,
     Table[FindRoot[SphericalBesselJ[l, x], {x, s}], {s, 1, xmax, 0.1}]] //
       Sort // Union // Select[#, Function[x, 0.1 <= x <= xmax]] &
 
 roots = Table[RootsOfSphericalBesselJ[l, 40], {l, 0, 4}]
 
 ScaledSphericalBesselJ[l_, n_, x_, xs_] := 
   SphericalBesselJ[l, x roots[[l + 1]][[n]]]/
     (FindMaximum[
       SphericalBesselJ[l, xx roots[[l + 1]][[n]]], {xx, xs}][[1]]) // FullSimplify
 
 Table[Plot[
   Table[ScaledSphericalBesselJ[l, n, x, 0.3], {l, 0, 3}] // 
     Evaluate, {x, 0, 1}], {n, 1, 4}] // 
       GraphicsRow[#, ImageSize -> 1336] &
 
 Table[Plot[
   Table[ScaledSphericalBesselJ[l, n, x, 0.3]^2, {l, 0, 3}] // 
     Evaluate, {x, 0, 1}], {n, 1, 4}] // 
       GraphicsRow[#, ImageSize -> 1336] &
 
 ScaledSphericalBesselJ[l_, n_, x_, xs_] := 
   SphericalBesselJ[l, x roots[[l + 1]][[n]]]/
     (FindMaximum[
       xx SphericalBesselJ[l, xx roots[[l + 1]][[n]]], {xx, xs}][[1]]) // FullSimplify
 
 Table[Plot[
   Table[x^2 ScaledSphericalBesselJ[l, n, x, 0.3]^2, {l, 0, 3}] // 
     Evaluate, {x, 0, 1}], {n, 1, 4}] // 
       GraphicsRow[#, ImageSize -> 1336] &

* 有限エネルギー障壁 [#kcddb8f1]

波動関数は

$$
\begin{cases}
A j_l\Big(\sqrt{\frac{2m}{\hbar^2}}\sqrt\varepsilon r\Big)&(r<a)\\
B h_l^{(1)}\Big(i\sqrt{\frac{2m}{\hbar^2}}\sqrt{V-\varepsilon} r\Big)&(a<r)\\
\end{cases}
$$

と表せるから、$r=a$ において

$$
A j_l\Big(\sqrt{\frac{2m}{\hbar^2}}\sqrt\varepsilon r\Big)=
B h_l^{(1)}\Big(i\sqrt{\frac{2m}{\hbar^2}}\sqrt{V-\varepsilon} r\Big)
$$

$$
A \frac{d}{dr}j_l\Big(\sqrt{\frac{2m}{\hbar^2}}\sqrt\varepsilon r\Big)=
B \frac{d}{dr}h_l^{(1)}\Big(i\sqrt{\frac{2m}{\hbar^2}}\sqrt{V-\varepsilon} r\Big)
$$

が成り立つ必要がある。$\rho_a=\sqrt{\frac{2m}{\hbar^2}}\sqrt\varepsilon a$, $\alpha=B/A$, $\beta=\sqrt{V-\varepsilon}/\sqrt{\varepsilon}$ と置けば、$d/dr=\sqrt{\varepsilon}d/d\rho$ などに注意して、

$$
j_l(\rho_a)=\alpha h_l^{(1)}(i\beta\rho_a)
$$

$$
j_l'(\rho_a)=
i\alpha \beta h_l^{(1)\prime}(i\beta\rho_a)
$$

両辺を割り算すると $\alpha$ を消去できて、

$$
\frac{j_l'(\rho_a)}{j_l(\rho_a)}=i\beta\frac{h_l^{(1)\prime}(i\beta\rho_a)}{h_l^{(1)}(i\beta\rho_a)}
$$

左辺は $-\tan r$ のように周期的に発散を繰り返す関数、右辺は $r$ に対して $-\infty$ から $0$ まで単調に増加する関数となる。

一方、右辺を $\beta$ の関数として見ると、$\beta=0$ で最大値 $-(l+1)/\rho_a$ を取る単調減少関数で、$\beta$ の大きいところで直線 $-b-1/r$ に漸近する。

したがって、左辺が $\frac{j_l'(\rho_a)}{j_l(\rho_a)}<-(l+1)/\rho_a$ を満たす時に限って上式を満たす $\beta$ が存在することになる。

距離を $a$、エネルギーを $\frac{\hbar^2}{2m}$ の単位に計ることにすると $\rho_a=\sqrt{\varepsilon}$ となり、

$$
\sqrt\varepsilon\,\frac{j_l'(\sqrt\varepsilon)}{j_l(\sqrt\varepsilon)}=
i\sqrt{V-\varepsilon}\,\frac{h_l^{(1)\prime}(i\sqrt{V-\varepsilon})}{h_l^{(1)}(i\sqrt{V-\varepsilon})}
$$

を得る。左辺と右辺とをそれぞれ $\varepsilon$ の関数としてプロットし、両者の交点を求めることでエネルギー固有値を計算できる。

&attachref(find_energy.svg,,40%);

ここでは左辺が青、右辺が黄色で表されている。

 LANG:mathematica
 djljl[l_, r_] = 
 D[SphericalBesselJ[l, r], r]/SphericalBesselJ[l, r] // FullSimplify
 dhlhl[l_, r_] = 
 D[SphericalHankelH1[l, I r], r]/SphericalHankelH1[l, I r] // 
  FullSimplify
 With[{l = 3, V = 250},
  Plot[{Sqrt[e] djljl[l, Sqrt[e]],
    Sqrt[V - e] dhlhl[l, Sqrt[V - e]]},
   {e, 0, V}]]
 
 With[{l = 3, V = 250},
  Table[
       e /. FindRoot[
          Sqrt[e] djljl[l, Sqrt[e]] -
           Sqrt[V - e] dhlhl[l, Sqrt[V - e]],
          {e, e0}] // Chop ,
       {e0, V/100, V - V/100, V/100}
       ] // Round[#, 10^(-9)] & // DeleteDuplicates // Sort // N
  ]
 
 PlotRadiusDistribution[l_, V_, ee_] := Table[
    Module[{a, n},
     a = SphericalBesselJ[l, Sqrt[e]]/
       SphericalHankelH1[l, I Sqrt[V - e]];
     n = NIntegrate[
       r^2 If[r < 1, SphericalBesselJ[l, Sqrt[e] r], 
          a SphericalHankelH1[l, I Sqrt[V - e] r]]^2, {r, 0, 10}];
     {e, a, n} // Chop
     ], {e, ee}] //
   Plot[
     {Table[e, {e, ee}],
      Table[e[[1]] + 10
         r^2 If[r < 1, SphericalBesselJ[l, Sqrt[e[[1]]] r], 
            e[[2]] SphericalHankelH1[l, I Sqrt[V - e[[1]]] r]]^2/e[[3]],
       {e, #}
       ],
      If[r < 1, 0, V]},
     {r, 0, 2}, Exclusions -> None, 
     PlotStyle -> {Dotted, Thick, Thick}, PlotRange -> {0, V + 50}] &
 RootsOfSphericalBesselJ[l_, xmax_] := 
  Map[Round[#[[1]][[2]], 0.00001] &, 
      Table[FindRoot[SphericalBesselJ[l, x], {x, s}], {s, 1, xmax, 
        0.1}]] // Sort // Union // 
   Select[#, Function[x, 0.1 <= x <= xmax]] &
 roots = Table[RootsOfSphericalBesselJ[l, 40], {l, 0, 4}]
 PlotRadiusDistribution2[l_, V_, m_] := 
  Table[
    {roots[[l + 1]][[k]], 
     NIntegrate[
      r^2 SphericalBesselJ[l, r roots[[l + 1]][[k]]]^2, {r, 0, 
       1}]}, {k, 1, m}] //
   Plot[
     {Table[e[[1]]^2, {e, #}],
      Table[e[[1]]^2 + 10
          If[r < 1, r^2 SphericalBesselJ[l, r e[[1]]]^2/e[[2]], 0],
      {e, #}
       ],
      If[r < 1, 0, V + 80]},
     {r, 0, 2}, Exclusions -> None, 
     PlotStyle -> {Dotted, Thick, Thick}, PlotRange -> {0, V + 50}] &
 GraphicsGrid[{{PlotRadiusDistribution[0, 250, {8.724234448`, 34.81866899`, 78.005855882`, 137.603946206`, 211.434813541`}],
    PlotRadiusDistribution2[0, 250, 5]},
   {PlotRadiusDistribution[1, 250, {17.836719112`, 52.556701517`, 104.118224984`, 171.426106714`, 248.24309799`}],
    PlotRadiusDistribution2[1, 250, 5]},
   {PlotRadiusDistribution[2, 250, {29.324751452`, 72.724643339`, 132.529167804`, 206.897892775`}],
    PlotRadiusDistribution2[2, 250, 4]},
   {PlotRadiusDistribution[3, 250, {43.076560923`, 95.225384656`, 163.090056824`, 242.876487902`}],
    PlotRadiusDistribution2[3, 250, 4]}
   }]
Counter: 4763 (from 2010/06/03), today: 3, yesterday: 3