原子核の相対論的 平均場模型のプログラムの改良

卒業論文発表会
2 月 1 日, 2011, 福井大学工学部物理工学科
原子核の相対論的
平均場模型のプログラムの改良
1. 動径グリッド の改良による誤差の低減
数理・量子科学講座 原子核理論研究室
物理工学科 西川恵子
原子核の構造
R = r0 A1/3
r0 ' 1.2 (fm)
A:質量数＝陽子の個数
＋中性子の個数
原子核の相対論的平均場模型
Dirac 方程式：スピン 12 のフェルミ粒子の量子状態を相対論的に記述する。
研究の目的
原子核の相対論的平均場模型プログラム（ 三和之浩修士論文、２００８年２月）
の改良
三和は動径波動関数を等間隔に配置したグリッド 点上の値で表現し 、
それらを Runge-Kutta 法により決定した。
しかし 、十分な精度を得るにはグリッド 間隔を非常に密にとらなければならず、
長い時間がかかった。
本研究では、
1. グリッド 間隔と誤差の関係
2. r 軸上のどこで大きな誤差が発生しているのか
3. 非等間隔グリッド の導入の効果
を調べた。
原子核の記述に置ける相対論的理論型式の利点
シュレディンガー方程式と違い、
1. スピン・軌道結合力が自然に導かれる。
2. 運動エネルギーが非常に大きくなる
高温・高圧下の原子核への外挿の信頼性が高い。
3. より基本的な理論（ 相対論的な場の理論）との関係がつく。
1. グリッド の間隔と原子核の全エネルギーの誤差の関係
核子の波動関数
テーラー展開解
Runge-Kutta法による解
漸近解
２５fm
Dirac 方程式
d
dr
(
G
F
)
(
=
− κr
µ−ε
µ+ε
κ
r
)(
G
F
)
m：核子の質量
V s：スカラーポテンシャル
Vv：ベクターポテンシャル
Ei：核子のエネルギー固有値
κ：角運動量と偶奇性に対応する
量子数で決まる定数
{
r
µ = m + Vs
ε = Ei − VV
全エネルギー E
∫
∑
1 ∞
E =
ni E i −
{−gσφσ (r) ρ s (r) + gωφω (r) ρv (r)
2 0
i
}
+gρφρ (r) ρ3 (r) + eA0 (r) ρ p (r) 4πr2dr
ni：i 番目の核子軌道を占拠する核子の個数
φσ 、φω 、φρ：中間子場
ρ s 、ρv 、ρ3：各種の核子密度
核の全エネルギーの誤差 ∆ E 208Pb N=126 0.1
0.01
|ΔE|(MeV)
0.001
0.0001
1e-05
1e-06
1e-07
1e-08
1e-09
0.001
0.01
Δr (fm)
∆r:グリッド 間隔（ fm ）
∆r=0.002fm での計算値を正確な値と見なした。
0.1
1
r 軸上のどこで大きな誤差が発生するか
Runge-Kutta 法で ∆r=0.05fm 進む際に生じる動径波動関数の誤差の大きさを
r の関数としてプロット
正確な値
数値解
F(r)
ΔF
0
r
000
√
r+Δr
(波動関数の誤差) ＝ (∆F)2 + (∆G)2
∆r ＝ 0.005fm で得た解を正確とみなした。
r
束縛エネルギーの大きい軌道と小さい軌道の比較
Double magic nucleus N=126 Z=82 , neutron orbitals
0.0001
0s1/2 (-64.80822MeV)
2d5/2 (-0.00800MeV)
1e-06
|error of wabefunction|
1e-08
1e-10
1e-12
1e-14
1e-16
1e-18
1e-20
1e-22
0
5
10
15
20
25
r (fm)
・束縛エネルギーによらず、r が小さいところで誤差は大きい。
・束縛エネルギーの小さい軌道では r が大きくなっても誤差は小さく、大きい軌道では
r が大きくなっても誤差は大きい。
安定核（ N= １２６、Z=８２）と、不安定核（ N= １７０、Z=８２）の比較
neutron 0s1/2 orbital (-64.80822MeV) (-62.53229MeV)
1e-06
N=126
N=170
|error of wabefunction|
1e-08
1e-10
1e-12
1e-14
1e-16
1e-18
1e-20
1e-22
0
5
10
15
20
25
r (fm)
・r が小さいところで誤差が大きいという傾向は変わらない。
・誤差の絶対値は、大きい r で不安定核の方がやや大きいが、傾向は同じである。
角運動量の大きい軌道と小さい軌道の比較（ 束縛エネルギーが同程度のもの）
Double magic nucleus N=126 Z=82 , neutron orbitals
1
0.01
|error of function|
0.0001
1e-06
0j15/2 (-0.40251MeV)
0j15/2 (-0.40251MeV)|wave func|
1e-08
3s1/2 (-0.09134MeV)
3s1/2 (-0.09134MeV)|wave func|
1e-10
1e-12
1e-14
1e-16
0
5
10
15
20
25
r (fm)
角運動量の大きい軌道では、r が小さいところで誤差は小さくなるが、これは遠心力ポ
テンシャルのせいで波動関数が小さいためであり、相対誤差として見れば、やはり小さ
い r で大きな誤差が発生するといえる。
非等間隔グリッド の導入
r
ξ
= f( )
Rs
Rs
∫ x
cosht
sinhx
dt
=
carcsinh
√
√
0
c2 + 1
1 + ( cosht )2
c
f’(x)
f(x)
cosh x
~等間隔
c
c
~等比数列
1
~等間隔
0
1.5
0
1
等間隔グリッド
Rs=0.01fm(固定
テーラー展開 Runge-Kutta法 漸近解
c=50,100,200,∞
積分は台形公式（＝この場合の最高精度公式）
全エネルギーの誤差 ∆ E
0.01
0.001
c=∞
c=200
ΔE(MeV)
0.0001
c=100
1e-05
1e-06
1e-07
1e-08
1e-09
100
1000
10000
Npt
N pt：Runge-Kutta 法で扱うグリッド 点の個数
100000
まとめ
１．208Pb 核のエネルギーの精度 ∆E は、それぞれ等間隔グリッド の間隔が
∆ r=0.10fm のとき、∆E=7.74keV
∆ r=0.005fm のとき、∆E=0.26keV
となり、(∆r)4.7 に比例すると分かった。
２．等間隔のグリッド では、r の小さいところで大きな誤差が発生することを見出した。
３．r の小さいところで密にグリッド 点をとれば 、少ない個数のグリッド 点で高精度の
結果を得られると期待して、最適なパラメータを探索中である。
卒業論文発表会
2 月 1 日, 2011, 福井大学工学部物理工学科
原子核の相対論的
平均場模型プログラムの改良
２自己無撞着解への収束を妨げる振動現象の解決策
物理工学科 酒井優
自己無撞着解とは
核子の波動関数とエネルギー
配位
1. 仮定し たポテンシャルをもとに核子の
エネルギーと波動関数を求める。
2. 求めたエネルギーを基に各軌道に占拠
核子数を割り当て配位を決定する。
3. 核子密度を得る。
核子密度
中間子場
核子の感じるポテンシャル
4. 中間子場の古典解を求める。
5. 中間子場を組み合わせてポテンシャル
を作る。
１と５のポテンシャルの差が非常に小さく
なれば、自己無撞着解が得られたことになる
占拠個数の振動現象とは
Ei
フェルミ面
• 配位を決定するときに、計算をさせるたびに 2 軌道の順番が入れ替わることがある。
• 入れ替わりがおこることで、計算を収束させることが出来ず自己無撞着解が求めら
れない。
振動の具体例
中性子 N=176 陽子 Z=82 の鉛の時に、軌道が入れ替わる事によって振動が
起こる（ 他にもおこる場合は多々ある）
N=174 Z=82
N=176 Z=82
-1466
-1.945
-1468
-1.95
single-particl levels(MeV)
total energy(MeV)
-1470
-1472
-1474
-1476
-1478
-1.955
-1.96
-1.965
-1.97
-1.975
-1.98
-1480
-1482
980
3 S 1/2
1 p 3/2
-1.985
985
990
iteration number
995
1000
-1.99
980
985
990
995
iteration number
iteration number：自己無撞着解を求めるための計算の繰り返しの回数
(iterationnumber ≤ 1000)
1000
振動を抑制する方法
1. 減衰因子の導入
2. 有限温度化（フェルミ準位を滑らかにする）
3. 占拠数の固定
＊現象論的に対相関力を導入して、BCS 近似解を作っても振動を止めることができる
が、本研究では、模型を変更せずにできる方法のみを調べた。
1. 減衰因子の導入：方法
自己無撞着解を求めるための計算を繰り返すときに、次の計算に使うポテンシャルを
Vnext として、今のポテンシャルを Vold 、計算で求めたポテンシャルを Vnew とおく。単
純に、
Vnext = Vnew
するのではなく、減衰因子 p を導入して、
Vnext = (1 − p) × Vold + p × Vnew (0 ≤ p ≤ 1)
とおきかえてみる。
p を小さくしていくことで振動を抑えることができると期待されるが、、、
1. 減衰因子の導入：結果
p = 0.1 の時
p = 0.03 の時
N=176 Z=82
N=176 Z=82
-1.948
-1.945
1 p 3/2
3 s 1/2
-1.949
1 p 3/2
3 s 1/2
single-particl levels(MeV)
single-particl levels(MeV)
-1.9455
-1.95
-1.951
-1.952
-1.953
-1.954
-1.946
-1.9465
-1.947
-1.955
-1.9475
-1.956
-1.957
950
960
970
980
iteration number
990
1000
-1.948
950
960
970
980
iteration number
振動の振幅を減少させるが消すことはできない
990
1000
2. 有限温度化：理論（ １）
フェルミ分布関数
kB ： ボルツマン定数
T ： 温度
: 一粒子準位
F ： フェルミ準位 (化学ポテンシャル )
1
1
, β=
f () =
1 + exp (β( − F ))
k BT
軌道 i の縮退度を gi 、エネルギー準位を i 、占拠核子数を ni とおき、
ni = gi f (i)
にとると、温度 T のときの分布になる。
フェルミ準位 F は粒子の個数を N として、
∑
N=
ni
i
を満たすように決定する。
1 は、反復で逆転のおきる２準位の間隔程度以上の大きさにとる必要がある。
β
占拠粒子個数
占拠粒子個数
2. 有限温度化：理論（２）
1
エネルギー準位
f (i) =
1
1 + exp (β(i − F ))
ni = gi f (i)
1
エネルギー準位
ni：占拠粒子個数
i：エネルギー準位
gi：縮退度
f (i)：フェルミ分布関数
フェルミ分布関数を乗じることで、準位逆転による振動を抑制できると期待される。
2. 有限温度化：結果
1 = 0.1(MeV) の時
β
3s 12 および 1p 23 準位
N=176 Z=82
N=176 Z=82
-1.955
-1473.4
3s1/2
1p3/2
-1.956
-1.957
single-particl levels(MeV)
total-energy(MeV)
-1473.5
-1473.6
-1473.7
-1473.8
-1.958
-1.959
-1.96
-1.961
-1.962
-1.963
-1.964
-1473.9
-1.965
-1474
50
100
150
iteration number
200
250
300
-1.966
50
100
150
200
iteration number
全エネルギーの差は k BT = 0.1(MeV) の時 ∆E = 70(KeV)
振動を止めることができるが、これは温度が０の場合の解ではない
250
300
3. 占拠数の固定：方法
振動現象がおきたらその後の配位を固定させて、自己無撞着解に収束させれば振動現象
を抑制できると期待できる。
自動でさせるようにプログラムを修正するのがベスト ですが、今回は N = 176 Z = 82
の場合のみについて手動で行うことにする。
１００回目以降で固定した場合と、１０１回目以降で固定した場合について考察する。
3. 占拠数の固定：結果
拡大図（ 配位固定後の部分）
N=176 Z=82
-1466
B
A
-1468
N=176 Z=82
-1473.75
B
A
-1473.8
-1472
total energy(MeV)
total energy(MeV)
-1470
-1474
-1476
-1478
-1473.85
-1473.9
-1480
-1473.95
-1482
50
60
70
80
90
100
110
iteration number
120
130
140
150
-1474
110
115
120
125
130
iteration number
(A):100 回目の反復以降、配位を凍結させる
(B):101 回目の反復以降、配位を凍結させる
配位を固定することで振動を抑制することができた
135
140
145
3. 占拠数の固定：考察
(A)
(B)
(A):100 回目の反復以降、配位を凍結させる
(B):101 回目の反復以降、配位を凍結させる
A の場合はエネルギーが低い軌道から
順に詰めた配位
B の場合は間に空いた軌道の挟まった
不安定な解
下の軌道は MeV で変化し ているのに対し
て、フェルミ面を挟む２軌道のエネルギー
準位の差は 20KeV ぐらいになっています。
まとめ
＜本研究の目的＞
自己無撞着解を求める計算のループの中でおこることがある振動現象をと
めて、計算を収束させる
1. 減衰因子の導入では、振動の振幅を抑えることができるが、消すことはできない。
2. フェルミ分布関数を導入して有限温度化すると振動を止めることができる。しかし
これは温度がゼロの場合の解ではない。
3. 振動が始まった後、配位を固定すると振動を止め自己無撞着解に収束させることが
できる。
全エネルギーが低い方の解はエネルギーの低い軌道から順に詰めた配位になって
いる。