2019 年度計算物理学２第 5 回レポート課題
2019/7/10 更新
1 次元調和振動子ポテンシャル内の粒子の運動を量子力学で考える。Schrödinger 方程式は()h̄2 d212 2−+mωxϕn (x) = En ϕn (x)2m dx22(1)で与えられる。En は固有エネルギー、ϕn (x) は固有波動関数である。この微分方程式の数値近似解を考える。遠方では波動関数はゼロに近づくため、有限の領域 |x| ≤ xmax のみで波動関数は値を持つとする。領域|x| ≤ xmax を N 点で等分割し、分割点を xi (i = 1, 2, · · · , N )xi = −xmax + (i − 1)∆x(2)とする。ただし ∆x = 2xmax /(N − 1) である。これらの分割点上での波動関数 ϕi,n = ϕn (xi ) を考える。(1)
式の Schrödinger 方程式は x = xi で−h̄2 d2 ϕi,n1+ mω 2 x2i ϕi,n = En ϕi,n2m dx22(3)となる。

1. 上の Schrödinger 方程式の二階微分に三点公式 (第 5 回配布資料、式 (11) を参照) を用いることによってこの方程式が En を固有値、(ϕn )t = (ϕ1,n , ϕ2,n , · · · , ϕN,n ) を固有ベクトルとする、N × N 行列の固有値問題となること N × N 行列を書いて示せ。ただし |x| ≤ xmax よりも外側での波動関数 ϕ0,n や ϕN +1,n などはゼロと近似してよい。

2. LAPACK のライブラリを用いて この固有値問題を解き、n = 0, 1, · · · , 4 の固有エネルギーおよび波動関数を求めよ。ただし簡単のため h̄ = 1, m = 1, ω = 1 とし、N = 1000、xmax はこの条件で x = xmax での波動関数の値が十分小さくなる適当な値を選ぶこと。(しかし xmax を大きくしすぎると分割の幅 ∆x が大きくなって近似が悪くなる) 解析解と比較して計算の近似の精度が気になる場合は二階微分の公式を点数の多い精度の良いものに置き換えても構わない。波動関数は規格化してグラフで表示せよ。波動関数の規格化は∫ ∞−∞dx|ϕk (x)|2 = 1(4)で与えられる (固有ベクトルの規格化とは異なることに注意) がこの積分は単純な和で置き換えてよい。

1. の式 (手書きで構いません。TeX が使える人は TeX で)、2. で作成したプログラムと固有エネルギーと波動関数のグラフを提出してください。
■ヒント• 三点公式を用いた場合は三重対角行列の固有値問題になり、LAPACK だと DSYEV や DGEEV を用いて対角化する。実対称行列用の DSYEV を用いた場合は固有値は昇順に配列に格納されるためこちらを使うことをすすめるが、実一般行列用の DGEEV を用いた場合は固有値が任意の順番で出てくるため、並べ替えを行う手間が発生する。

• DGEEV で並べ替えをしたい場合は、Fortran にはソートの組み込み関数はないので、並べ替えのアルゴリズムを考えて自分で行う。組み込み関数で使えそうなものには MINVAL, MAXVAL(配列の指定された領域での最小値、最大値を返す) や、MINLOC, MAXLOC(配列の指定された領域での最小値、最大値を与えるインデックスを返す) などを使ってもよい。

1

• 1 次元調和振動子の問題は解析解があるので計算結果がある程度の精度で一致しているかチェックする。

(以下は発展問題。余力のある方だけどうぞ)
3. 次に 3 次元調和振動子を考える。3 次元調和振動子ポテンシャル VHO (r) = 21 mω 2 r2 は中心力であるため波動関数 ϕn (r, ϕ, φ) は動径部分 Rnl (r) と角度部分 Ylm (θ, ϕ) の積に分けることができ、角度部分の解は球面調和関数 Ylm (θ, ϕ) で与えられる。動径部分は(−)
1 1 d2h̄2 l(l + 1)
r++V(r)Rnl (r) = Enl Rnl (r)HO2m r dr22mr2(5)で与えられる。l は軌道角運動量量子数で l = 0, 1, 2, · · · である。r ≤ rmax を N 点で分割して固有値問題としてこの微分方程式の解を求めてみる。rmax では動径波動関数はゼロと近似し、原点付近では Rnl (r) ∼ r lと振る舞うことを使って良い。ただし原点では l ̸= 0 のとき遠心力ポテンシャル項 (∼ 1/r2 ) が発散するため、原点を含まない離散点を取るようにする。

h̄ = 1, m = 1, ω = 1 の場合に動径部分の方程式を解き、l = 0, 1, 2, 3 それぞれの場合の固有エネルギー Enlを n = 0, 1, 2, · · · 4 を求めよ。どの n, l に対応する解が縮退するか。

4. 原子核内の陽子や中性子は核力によって複雑に相互作用しているにもかかわらず、中心力の WoodsSaxon 型のポテンシャル中を独立に運動している粒子としてよく近似できることが知られている。WoodsSaxon ポテンシャルは
VWS (r) =V(0)r − R01 + expa(6)で与えられる。rmax = 20 fm, h̄c = 197 MeV fm2 , mc2 = 939 MeV, V0 = −45.5 MeV, R0 = 6.26 fm,

a = 0.67 fm の場合 (c は光速だが値は使わない)、(5) 式の調和振動子ポテンシャル VHO (r) を Woods-Saxon
ポテンシャル VWS (r) に置き換えた動径方向の Schrödinger 方程式を解き、負の固有エネルギーをすべて求めよ。3. の 3 次元調和振動子の解と対応させたとき、調和振動子で縮退していた解はどうなっているか。

5. 陽子や中性子はスピン s = 1/2 のフェルミ粒子であり、パウリの排他律により複数の同種の粒子が一つの固有状態にいることが許されないため、エネルギーの低い順に固有エネルギー状態を占有する。軌道角運動量 l を持つ準位の縮退度は m = l, l − 1, · · · , −l + 1, −l の 2l + 1 で、さらに各準位はスピンの自由度に関して二重に縮退 (sz = ±1/2) しているとすると、3. のポテンシャル問題ではちょうど粒子が 2, 8, 20, 40, 70 個詰まった時に次に空いているエネルギー準位との間にエネルギーギャップができる閉殻構造となり安定化する。
その一方で実験的には 2, 8, 20, 28, 50, 82,... において閉殻構造になることが知られており (魔法数と呼ばれる)、4. の Woods-Saxon ポテンシャルでも 28, 50, 82 は再現できない。これを再現するためには中心力に加えてスピン軌道力が必要となる。

VSO (r) = ηSO

1 dVWS (r)l̂ · ŝrdr2(7)

スピン軌道力を導入して原子核の魔法数を説明した Mayer と Jensen は 1963 年にノーベル賞を受賞した。4.の Woods-Saxon ポテンシャルにスピン軌道ポテンシャルを加えて魔法数の再現を試みてみよう。
スピン軌道ポテンシャルの l̂ · ŝ 項は全角運動量演算子 ĵ = l̂ + ŝ を導入すると

l̂ · ŝ =1 2(ĵ − l̂2 − ŝ2 )2(8)と書けるため、角運動量を合成して固有波動関数を ϕn,l,m,s,sz から ϕn,j,jz ,l,s に変更すると、この波動関数に対して[()]
1
1 1(l̂ · ŝ)ϕn,j,jz ,l,s =j(j + 1) − l(l + 1) −+ 1 ϕn,j,jz ,l,s22 2(9)と j と l のみで表せる。j は l ± 12 である (ただし l = 0 の場合は j = 1/2 のみ)。波動関数 ϕn,j,jz ,l,s (r) の動径波動関数 Rnjl (r) に関する Schrödinger 方程式()
h̄2 1 d2h̄2 l(l + 1)−++V(r)+V(r)Rnjl (r) = Enjl Rnjl (r)WSSO2m r dr22mr2(10)を解いて負の固有エネルギー En,j,l をすべて計算し、閉殻構造の現れる領域を議論せよ。ただし、ηSO = −0.7とし、状態 ϕn,j,jz ,l,s の縮退度は jz = j, j − 1, · · · , −j + 1, −j の 2j + 1 である。
■ヒント

• R(r) =χ(r)と置くと微分方程式は簡単になる。
r• 3 次元調和振動子も解析解があるので比較して計算のチェックをする。
解いた方は 3.-5. で作成したプログラムと求めた固有エネルギーのリストなどを提出してください。

3