卒業論文発表会
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 方程式

ddr(GF)(=− κrµ−εµ+εκr)(GF)

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 0i}+gρφρ (r) ρ3 (r) + eA0 (r) ρ p (r) 4πr2dr

ni：i 番目の核子軌道を占拠する核子の個数
φσ 、φω 、φρ：中間子場
ρ s 、ρv 、ρ3：各種の核子密度

核の全エネルギーの誤差∆E208 Pb

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
0r000√r+Δr(波動関数の誤差) ＝ (∆F)2 + (∆G)2∆r ＝ 0.005fm で得た解を正確とみなした。

r束縛エネルギーの大きい軌道と小さい軌道の比較

Double magic nucleus N=126 Z=82 , neutron orbitals
0.00010s1/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

0510

152025r (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

0510152025

r (fm)

・r が小さいところで誤差が大きいという傾向は変わらない。
・誤差の絶対値は、大きい r で不安定核の方がやや大きいが、傾向は同じである。

角運動量の大きい軌道と小さい軌道の比較（ 束縛エネルギーが同程度のもの）

Double magic nucleus N=126 Z=82 , neutron orbitals
10.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-160510152025r (fm)

角運動量の大きい軌道では、r が小さいところで誤差は小さくなるが、これは遠心力ポテンシャルのせいで波動関数が小さいためであり、相対誤差として見れば、やはり小さい r で大きな誤差が発生するといえる。

非等間隔グリッド の導入rξ= f( )
RsRs∫ xcoshtsinhx
dt=carcsinh√√0
c2 + 1
1 + ( cosht )2cf’(x)f(x)

cosh x~等間隔cc

~等比数列
1~等間隔

01.501
等間隔グリッド

Rs=0.01fm(固定テーラー展開 Runge-Kutta法 漸近解

c=50,100,200,∞
積分は台形公式（＝この場合の最高精度公式）

全エネルギーの誤差∆E0.010.001c=∞c=200ΔE(MeV)

0.0001

c=100

1e-05
1e-06
1e-07
1e-08
1e-09
1001000

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-1482980

3 S 1/2
1 p 3/2-1.985985990
iteration number

9951000-1.99980985990995

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.9451 p 3/2
3 s 1/2-1.9491 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.957950960970

980

iteration number

990

1000

-1.948950960970980

iteration number

振動の振幅を減少させるが消すことはできない

990

1000

2. 有限温度化：理論（ １）フェルミ分布関数kB ： ボルツマン定数
T ： 温度 : 一粒子準位
F ： フェルミ準位 (化学ポテンシャル )11, β=f () =1 + exp (β( − F ))k BT

軌道 i の縮退度を gi 、エネルギー準位を i 、占拠核子数を ni とおき、ni = gi f (i)にとると、温度 T のときの分布になる。
フェルミ準位 F は粒子の個数を N として、∑N=niiを満たすように決定する。
1 は、反復で逆転のおきる２準位の間隔程度以上の大きさにとる必要がある。
β

占拠粒子個数

占拠粒子個数2. 有限温度化：理論（２）1エネルギー準位f (i) =11 + 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

100150
iteration number

200

250300-1.966

50100150200

iteration number

全エネルギーの差は k BT = 0.1(MeV) の時 ∆E = 70(KeV)
振動を止めることができるが、これは温度が０の場合の解ではない

2503003. 占拠数の固定：方法

振動現象がおきたらその後の配位を固定させて、自己無撞着解に収束させれば振動現象を抑制できると期待できる。

自動でさせるようにプログラムを修正するのがベスト ですが、今回は N = 176 Z = 82の場合のみについて手動で行うことにする。
１００回目以降で固定した場合と、１０１回目以降で固定した場合について考察する。

3. 占拠数の固定：結果

拡大図（ 配位固定後の部分）

N=176 Z=82-1466
BA-1468

N=176 Z=82-1473.75
BA-1473.8-1472
total energy(MeV)

total energy(MeV)-1470-1474-1476-1478-1473.85-1473.9-1480-1473.95-1482
5060708090100110

iteration number

120130140150-1474
110115120125130

iteration number(A):100 回目の反復以降、配位を凍結させる(B):101 回目の反復以降、配位を凍結させる

配位を固定することで振動を抑制することができた

135140145

3. 占拠数の固定：考察

(A)

(B)

(A):100 回目の反復以降、配位を凍結させる(B):101 回目の反復以降、配位を凍結させる

A の場合はエネルギーが低い軌道から順に詰めた配位

B の場合は間に空いた軌道の挟まった不安定な解下の軌道は MeV で変化し ているのに対して、フェルミ面を挟む２軌道のエネルギー準位の差は 20KeV ぐらいになっています。

まとめ＜本研究の目的＞
自己無撞着解を求める計算のループの中でおこることがある振動現象をとめて、計算を収束させる

1. 減衰因子の導入では、振動の振幅を抑えることができるが、消すことはできない。

2. フェルミ分布関数を導入して有限温度化すると振動を止めることができる。しかしこれは温度がゼロの場合の解ではない。

3. 振動が始まった後、配位を固定すると振動を止め自己無撞着解に収束させることができる。
全エネルギーが低い方の解はエネルギーの低い軌道から順に詰めた配位になっている。