球対称平均場模型と近似的な変形の扱いによる原子核質量公式

2009 年 2 月
福井大学 大学院 工学研究科 物理工学専攻
平成 19 年度入学 13 番 山田昌平

目次
序章

3

第 1 章 原子と原子核の質量1.1 原子質量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 原子質量超過 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 原子と原子核の束縛エネルギー . . . . . . . . . . . . . . . . . . . . . .
1.4 原子核の質量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 評価済みの原子質量データ . . . . . . . . . . . . . . . . . . . . . . . .
1.6 質量が既知の核種 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.7 導入した略号のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
5
5
6
7
9
10

第 2 章 原子核の構造
2.1 液滴模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Weizsäcker-Bethe の質量公式 . . . . . . . . . . . . . . . . . .
2.1.2 質量公式のパラメータの決定 . . . . . . . . . . . . . . . . . .
2.1.3 液滴模型から得た値と実験値との比較 . . . . . . . . . . . . .
2.2 殻模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 魔法数の存在 . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 スピン軌道結合力の導入 . . . . . . . . . . . . . . . . . . . . .

11
11
11
14
15
17
17
17

第 3 章 変形を考慮した質量公式の作成
3.1 変形核と四極子モーメント . . . . . . . . . . . . . . . . . . . . . . . .
3.2 変形核の殻エネルギー . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.1 球形核の重ね合わせとしての変形核 . . . . . . . . . . . . . .
3.2.2 平均変形エネルギー . . . . . . . . . . . . . . . . . . . . . . .
3.2.3 変形を含めた原子核の質量の計算 . . . . . . . . . . . . . . . .

22
22
23
24
37
38

第 4 章 平均場模型
4.1 Skyrme-Hartree-Fock 法 . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Hartree-Fock 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Skyrme 相互作用 . . . . . . . . . . . . . . . . . . . . . . . .
4.3 相対論的平均場模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 核子の媒体としての有効中間子場 . . . . . . . . . . . . . . . .
4.3.2 中間子の古典場近似 . . . . . . . . . . . . . . . . . . . . . . .
4.3.3 負エネルギー状態の核子の無視 . . . . . . . . . . . . . . . . .

40
40
40
43
45
45
46
47

1

第 5 章 新しい質量公式49

第 6 章 まとめ

53

参考文献

55

謝辞

57

付録 Program List
6.0.4 変形核を取り扱う大切な方法 . . . . . . . . . . . . . . . . . .
6.0.5 球形核の重ね合わせとしての変形核 . . . . . . . . . . . . . .

58
58
60

付録

64

Program List

2

序章
我々の身の回りの物質、又我々自身も含め、全て元素からできている。その元素は宇宙の誕生のビッグバン以来の約 137 億年という長い歴史の中で，原子核反応により合成されたと考えられている。この原子核反応は恒星のエネルギー源であるともともに超新星爆発など、他の天体現象でも主要な役割を果しており、宇宙の歴史をたどり
物質の起源を探るためにも原子核の性質を詳し く知る必要がある。
原子核の持つ最も基本的な量に原子質量がある。質量とエネルギーは同等であり、原子質量は静止状態の原子の全エネルギーである。それ故に、原子、原子核の質量は原子核の安定性や、崩壊や、反応の起こりやすさに決定的な影響を持っている。原子核の質量を知ることは、α 線や β 線の放出エネルギーや核子分離エネルギーなどを知る
ことに直接つながる。
現在、実験的に質量が確認されている原子核は約 2,200 種である。ところが理論的には、それ以外に何千種もの原子核が存在しうると予想されている。こうした原子核の特徴を知ろうとする場合、まず問題にしなければならないのが原子核の質量であり、原子核質量公式はそれを理論的に予想する手段である。質量公式は、原子核の質量を中性子の個数と陽子の個数の関数として与えるものである。最も簡単な質量公式のひとつとして、広く知られている Weizsäcker-Bethe の質量公式がある。今日では、巨視
的 (液滴模型) および微視的 (殻模型) 手法を組み合わせた様々な質量公式が提案されていて、質量の実験値がない原子核の質量推定に使用されている。[15]
本論文では、近似的な変形の扱いによる原子核質量公式を作成するにあたり、核変形を考慮しないで済めば 、計算量の激減により、非常に広範なパラメータの最適化を行うことが可能になる。そこで本研究では、2000 年に小浦寛之氏・橘孝博氏・宇野正宏
氏・山田勝美氏によって発表された論文 [19] に基づく KUTY 公式の変形核を球形核の重畳と見る近似を採用する。そして、様々な球対称平均場模型と組み合わせて新たな質量公式の改良を試みる。

3

第1章

原子と原子核の質量

本章に記した物理定数は全て参考文献 [2] によった。また、本論文では質量は光速
c の２乗を乗じて得られるエネルギー単位の値で表す。

1.1

原子質量

原子質量 (atomic mass) とは、原子核のまわりに電子がついて電気的に中性となった原子の質量のことを指し 、原子核の質量 (nuclear mass) とは、電子が全くついていない裸の原子核質量のことを指す。[3][4]
一般に、原子番号 (atomic number) を Z, 中性子数 (neutron number) を N で表す。Z は陽子の個数 (proton number) に等しい。また A = N + Z を質量数 (massnumber) と言う。この論文では、原子の質量に AM(N , Z ) という記号を用い、原子核の質量を表すのに NM(N , Z ) という記号を用いることにする。
AM(N , Z ) と NM(N , Z ) との間には、下記の関係式が成り立つ。

AM(N , Z ) = NM(N , Z ) + Zme − kel Z 2.39

(1.1)

これは、原子核質量に陽子と同じ個数の電子の静止質量を加え、電子の束縛エネルギー(binding energy) を引いたものである。ここで、電子の束縛エネルギーは kel × Z 2.39でよく近似され、比例定数の値は、

kel = 1.433 × 10−5 MeV(1.2)

である。そして、電子の静止質量 (electron mass) は me と表してあり、

me = 5.485799110 × 10−4 u

(1.3)

である。
例として (1.1) 式を最も軽い原子である水素原子に適用すると、この原子の原子核は1 個の陽子 (proton) であるから、水素原子の質量 mH は陽子質量 (proton mass)mpと下式のように関係付けられる。

mH = mp + me − kel 12.39= mp + me − kel ,(1.4)

mp = 1.00727646688u.

(1.5)

ただし、水素原子の束縛エネルギーは正確には 13.6eV であり、kel = 14.33eV とは 0.7keVの誤差があるが、展型的には数百 keV、将来的な目標として 100keV と言われているので、この論文で議論する原子核質量公式の誤差は無視できる値である。
41.2

原子質量超過原子質量はほとんど質量数に比例するので、その比例関係からのズレ の部分をデータとし て与える習慣がある。これを質量超過 (mass excess) と呼び 、本論文では AME(N , Z ) で表し 、AME(N , Z ) = AM(N , Z ) − (N + Z )u(1.6)

で定義される。ここで、u は原子質量単位 (atomic mass unit) と呼ばれる。現行の定義では「質量数 12 の炭素原子 1 個の質量の 12 分の 1」であり、

u = 931.494013 MeV

(1.7)

である。

1.3

原子と原子核の束縛エネルギー

まず束縛エネルギーとは、粒子が集まって束縛系をつくっているとき、これらの粒子を引き離してばらばらにしてしまうために必要なエネルギーを言う。束縛系の全エ
ネルギーを M 、それを構成する各粒子のエネルギーを mi とすると、束縛エネルギー
BE は、

BE =

n
∑

mi − M(1.8)

i=1で与えれる。そして、原子の束縛エネルギー (atomic binding energy) とは、中性原子を N 個の中性子と Z 個の中性な水素原子にバラバラに分けるのに必要なエネルギー
であり、この論文では ABE(N , Z ) で表す。式で表すと、ABE(N , Z ) = Nmn + ZmH − AM(N , Z )(1.9)と表される。ここで、mn は中性子質量 (neutron mass) で、

mn = 1.00866491578u(1.10)

であり、mH は水素原子質量 (hydrogen atom mass) で、

mH = 1.007825032u(1.11)である。
一方、原子核の束縛エネルギー (nuclear binding energy) とは、電子の付いていない裸の原子核を N 個の中性子と Z 個の陽子に分けるのに必要なエネルギーのことであり、この論文では NBE(N , Z ) で表す。式で表すと、NBE(N , Z ) = Nmn + Zmp − NM(N , Z )(1.12)となる。(1.1) 式よりNM(N , Z ) = AM(N , Z ) − Zme + kel Z 2.39
5

(1.13)

となる。また、(1.6) 式よりAM(N , Z ) = AME(N , Z ) + (N + Z )u

(1.14)となり、(1.14) 式を (1.1) 式に代入するとNM(N , Z ) = AME(N , Z ) + (N + Z )u − Zme + kel Z 2.39

(1.15)

ここで、上式を (1.12) 式に代入するとNBE(N , Z ) = N mn + Zmp − AME(N , Z ) − (N + Z )u + Zme + kel Z 2.39
= (mn − u)N + (mp − u)Z − AME(N , Z ) + Zme + kel Z 2.39
= (mn − u)N + (mp + me − u)Z + kel Z 2.39 − AME(N , Z ).
(1.16)(1.4) 式より
mp + me = mH + kel(1.17)

上式を (1.16) 式に代入するとNBE(N , Z ) = (mn − u)N + (mH + kel − u)Z − kel Z 2.39 − AME(N , Z )= (mn − u)N + (mH − u)Z − kel (Z 2.39 − Z ) − AME(N , Z )
となる。

1.4

原子核の質量本研究で必要とする量は原子核の質量 NM(N , Z ) であるが、実験値は原子に関係する量 AME(N , Z ) で与えられることが多い。そこで、NM(N , Z ) を AME(N , Z ) を使って表す式を求めと、

NM(N , Z ) = Nmn + Zmp − NBE(N , Z )となる。次に (1.18) 式を代入すると、NM(N , Z ) = N mn + Zmp − (mn − u)N − (mH − u)Z+ kel (Z 2.39 − Z ) + AME(N , Z )= uN + (mp − mH + u)Z + kel (Z 2.39 − Z ) + AME(N , Z )= uN + (kel − me + u)Z + kel (Z 2.39 − Z ) + AME(N , Z )= uN + (u − me )Z + kel Z 2.39 + AME(N , Z )(1.18)となる。
この論文では、等式 (1.18) を使用して実験データから原子核の質量を求める。
61.5

評価済みの原子質量データ

原子質量の実験値としては G.Audi と A.H.Wapstra と C.Thibault によりまとめられた評価済み原子質量データ表 [13] の値が広く用いられている。
この論文では彼らが公開しているデータファイルのうち『 mass.mas03 』(2003 年版
原子質量推奨値表) を用いる。図 1.1 にデータファイルの一部を示す。その各コラムの表す量について、左端から順番に以下で説明する。

• 第 1 コラムは、印刷時のページおくりの指示 (フォート ラン書式) であり、データではない。
• 第 2 コラムは、『 N-Z 』で中性子数から陽子数を引いたものである。
• 第 3 コラムは、『 N 』で中性子数 (neutron number) という。
• 第 4 コラムは、『 Z 』で陽子数 (proton number) という。
• 第 5 コラムは、『 A 』で原子核を構成する核子の数であり、質量数 (mass number)という。
• 第 6 コラムは、『 EL 』で元素記号である。
• 第 7 コラムは、『 O 』で二次核種の値の起源を表す。
• 第 8 コラムは、『 MASS EXCESS 』で原子質量超過であり、この論文では AM E(N, Z)と表すことにする。
• 第 9 コラムは、『 MASS EXCESS 』の誤差である。
• 第 10 コラムは、『 BINDING ENERGY 』で 1 核子当たりの原子の束縛エネル
ギーである。
• 第 11 コラムは、『 BINDING ENERGY 』の誤差である。
• 第 12 コラムは、β 崩壊のことで β − 崩壊と β + 崩壊がある。
• 第 13 コラムは、『 BETA-DECAY ENERGY 』で β 崩壊エネルギーである。
• 第 14 コラムは、『 ATOMIC MASS 』で原子質量である。
• 第 15 コラムは、『 ATOMIC MASS 』の誤差である。

7

1
0
0

a0peysza********************** file : mass.mas03 **********************A T O M I C
M A S SA D J U S T M E N T
DATE 18 Nov 2003 TIME 14:12
A=
0 TO 293

This is one file out of a series of 3 files published in:"The Ame2003 atomic mass evaluation (II)" by G.Audi, A.H.Wapstra and C.Thibault
Nuclear Physics A729 p. 337-676, December 22, 2003.
for files : mass.mas03 atomic masses
rct1.mas03 react and sep energies, part 1
rct2.mas03 react and sep energies, part 2
A fourth fileis the "Rounded" version of the atomic mass table (the first file)mass.mas03round atomic masses "Rounded" version

All files are 3179 lines long with 124 character per line.
Headers are 39 lines long.
Values in the files 1, 2 and 3 are exact (unrounded) copy of the published ones
Values in the files 4
are exact
copy of the published ones
col 1
format::Warnings:

Fortran character control: 1 = page feed 0 = line feed
a1,i3,i5,i5,i5,1x,a3,a4,1x,f13.5,f11.5,f11.3,f9.3,1x,a2,f11.3,f9.3,1x,i3,1x,f12.5,f11.3,1x
cc NZ N Z A
el omass unc binding unc
B beta unc
atomic_mass
unc
this format is not identical to the one used in Ame1993 and Ame1995
in particular "Mass Excess" and "Atomic Mass" values are given now, when necessary,with 5 digits after decimal point.
also, the "Binding Energy" column is replaced by a "Binding Energy per A" one.
decimal point is replaced by # for (non-experimental) estimated values.
* in place of value : not calculable....+....1....+....2....+....3....+....4....+....5....+....6....+....7....+....8....+....9....+...10....+...11....+...12....
MASS LIST
for analysis
1N-Z

N

Z

A

0

1
0
1
2
1
0
3
2
1
4
3
2
1
5
4
3
2
1
6
5
4
3
2
6
5
4
3
2
7
6
5
4
3
8
7
6
5
4
3

0
1
1
1
2
3
1
2
3
1
2
3
4
1
2
3
4
5
1
2
3
4
5
2
3
4
5
6
2
3
4
5
6
2
3
4
5
6
7

1 n
1 H
2 H
3 H
3 He
3 Li
4 H
4 He
4 Li
5 H
5 He
5 Li
5 Be
6 H
6 He
6 Li
6 Be
6 B
7 H
7 He
7 Li
7 Be
7 B
8 He
8 Li
8 Be
8 B
8 C
9 He
9 Li
9 Be
9 B
9 C
10 He
10 Li
10 Be
10 B
10 C
10 N

1
-1
0 0
0 1
-1
-3
0 2
0
-2
0 3
1
-1
-3
0 4
2
0
-2
-4
0 5
3
1
-1
-3
0 4
2
0
-2
-4
0 5
3
1
-1
-3
0 6
4
2
0
-2
-4

EL

O

-pp
-n
-p
-nn
-n
-p
x
-3n

x
-nn
-n

+3n

4n

-pp
++
-n

--

MASS EXCESS(keV)
8071.31710
0.000537288.97050
0.0001113135.72158
0.00035
14949.80600
0.0023114931.21475
0.00242
28667#
2000#
25901.518
103.2862424.91565
0.00006
25323.185
212.132
32892.440
100.000
11386.233
50.00011678.886
50.00037996#
3996#
41863.757
264.90617595.106
0.755
14086.793
0.015
18374.947
5.44843603#
699#
49135#
1005#
26101.038
16.65814908.141
0.079
15770.034
0.106
27868.346
70.71231598.044
6.868
20946.844
0.095
4941.672
0.035
22921.490
1.00035094.060
23.06840939.429
29.41824954.264
1.935
11347.648
0.39812415.681
0.983
28910.491
2.13848809.203
70.00133050.581
15.12412606.670
0.40112050.731
0.386
15698.682
0.403
38800.148
400.000

BINDING ENERGY/A
BETA-DECAY ENERGY
(keV)
(keV)
0.0
0.0
B782.347
0.001
0.0
0.0
B*
1112.283
0.000 B*
2827.266
0.001 B18.591
0.001
2572.681
0.001 B- -13736#
2000#-2267#
667#
B*
1400.351
25.821 B- 23476.602 103.286
7073.915
0.000 B- -22898.270 212.132
1153.761
53.033 B*
1336.360
20.000 B- 21506.207 111.803
5481.132
10.000 B-292.653
70.711
5266.132
10.000 B- -26317#
3996#
-154#
799#
B*
963.633
44.151 B- 24268.651 264.908
4878.017
0.126 B3508.313
0.755
5332.345
0.002 B- -4288.154
5.448
4487.262
0.908 B- -25228#
699#
152#
116#
B*
940#
144#
B- 23034#
1005#
4119.070
2.380 B- 11192.898
16.6585606.291
0.011 B-861.893
0.071
5371.400
0.015 B- -12098.312
70.7123531.306
10.102 B*
3925.975
0.859 B- 10651.200
6.869
5159.582
0.012 B- 16005.172
0.101
7062.435
0.004 B- -17979.819
1.001
4717.164
0.125 B- -12172.570
23.0893097.800
2.883 B*
3348.637
3.269 B- 15985.165
29.4815037.839
0.215 B- 13606.616
1.928
6462.758
0.044 B- -1068.034
0.899
6257.160
0.109 B- -16494.809
2.353
4337.476
0.238 B*
3033.927
7.000 B- 15758.623
71.6154531.555
1.512 B- 20443.910
15.1236497.711
0.040 B555.939
0.557
6475.071
0.039 B- -3647.951
0.118
6032.041
0.040 B- -23101.466 400.000
3643.660
40.000 B*

図 1.1:ATOMIC MASS
(micro-u)
1 008664.91574
0.00056
1 007825.03207
0.00010
2 014101.77785
0.00036
3 016049.27767
0.00247
3 016029.31914
0.00260
3 030775#
2147#
4 027806.424
110.881
4 002603.25415
0.00006
4 027185.558
227.733
5 035311.488
107.354
5 012223.624
53.677
5 012537.800
53.6775 040790#
4290#
6 044942.594
284.388
6 018889.124
0.810
6 015122.794
0.016
6 019726.317
5.848
6 046810#
750#
7 052749#
1079#
7 028020.618
17.883
7 016004.548
0.084
7 016929.828
0.113
7 029917.901
75.9128 033921.897
7.373
8 022487.362
0.101
8 005305.103
0.037
8 024607.233
1.073
8 037675.025
24.764
9 043950.286
31.581
9 026789.505
2.077
9 012182.201
0.426
9 013328.782
1.055
9 031036.689
2.295
10 052398.837
75.14910 035481.259
16.236
10 013533.818
0.430
10 012936.992
0.413
10 016853.228
0.432
10 041653.674
429.417

G.Audi と A.H.Wapstra と C.Thibault に よ り まと められ た 原 子 質 量デ ー タ ファイル
『 mass mas03 』の 1 部分

8

1.6

質量が既知の核種

原子質量表データの『 MASS EXCESS 』、『 BINDING ENERGY 』、『 ATOMICMASS 』など中にある＃が付うているデータは実験データではなく、一部の系統的な傾向から出されたものである。今回の研究ではより信頼の置ける値を用いるためそのデータは飛ばすことにする。
この条件を付けて『 mass mas03 』を読み込むと、2226 個の (N 、Z) のデータが得られた。図 1.2 に示したものは横軸を中性子数、縦軸を陽子数にとった実験データの存在する核種である。

120
known nuclei
Proton Number Z

100
80
60
40
20
0
0

20

40 60 80 100 120 140 160
Neutron Number N

図 1.2: 2226 個の質量が既知の原子核 横軸:中性子数 縦軸:陽子数

9

1.7

導入した略号のまとめ

この章で導入した略号のま記号や略称を下記にまとめておく。

A ： mass number
ABE ： atomicbinding energy
AM ： atomic mass
AME ： atomic massexcess
kel ： electric binding energy
me ： electron massmH ： hydrogen atom mass
mn ： neutron mass
mp ： proton mass
N ： neutron number
NBE ： nuclear bindingenergy
NM ： nuclear massu ： atomic massunit
Z ： proton number10

質量数原子の束縛エネルギー
原子の質量
原子質量超過電子の束縛エネルギー
電子質量
水素原子質量
中性子質量
陽子質量
中性子数
原子核の束縛エネルギー
原子核の質量
原子質量単位
陽子数

第2章

原子核の構造液滴模型

2.1

前章で求めた原子核の質量を用いて 1 核子当たりの束縛エネルギー NBE(N , Z )/(N +Z ) を求めると質量数との関係は図 2.1 のようになる。横軸が質量数 A = N + Z 、縦軸が 1 核子当たりの束縛エネルギーである。
9
Binding Energy BE/A(MeV)

8
7
6
5
4
3
2
1
0
-1
0

50

100

150
Mass Number A

200

250

300

図 2.1: 質量数 A が原子核について核子当たりの束縛エネルギーを示したもの

軽い原子核を別とすれば 核子当たりの束縛エネルギーはほぼ 8MeV である。このことは 1 核子当たりの束縛エネルギーは質量にほとんどよらないことを意味している。
この現象を結合エネルギーの飽和性と呼ぶ [1]。
自然界で飽和性を示す最も身近な系は液滴である。このことから原子核を考える時には液滴模型 (liquid-drop model) が用いられ、この様な模型を原子核の液滴模型と呼ぶ。この模型を用いて原子核の質量を表す公式を考えてみる。

2.1.1

Weizsäcker-Bethe の質量公式

質量公式 (mass formula) とは原子核の質量を陽子数 Z 個と中性子数 N の関数として表したものである。この章では、最も簡単で基本的とされる Weizsäcker-Bethe の
質量公式 [8][9] を取り扱う。それを用いると、Z 個の陽子と N 個の中性子を持つ原子
11核の質量は次のような現象論的な式で与えられる。

M (N, Z) = Mvol + Msurf + MC + Msym + Mpair(2.1)

この式の 5 つの項はそれぞれ体積項 Mvol 、表面項 Msurf クーロン項 MC 、対称項 Msym 、
偶奇項 Meo と呼ばれる。

(1) 体積項
液滴中の粒子間に働く力は到達距離が短いため、各粒子はそれぞれの近くにいる粒子とだけ引っ張り合って結合しているので、表面を除外すれば 、結合の強さは粒子数に比例する。そこで、結合エネルギーを表す公式の最大項として、核子数すなわち質量数に比例する項を導入する。密度一定の液体では粒子数は液滴の体積に比例するので、この項は体積項と呼ばれる。原子核でも飽和性の故に、わずかな例外を除けばどの原子核の密度も中央部において同じである。このことにより、原子核の体積 V は、
V = (4π/3)R3 、R = r0 A1/3 で表せる。こうすれば V は A に比例し 、密度 V /A は一定であるからである。avol を体積に比例するパラメータとするとMvol = avol A

(2.2)

になる。

(2) 表面項
液滴の表面にある粒子は引っ張り合う相手の粒子が少ないので、結合エネルギーは少さい。したがって、液滴全体のエネルギーは表面積に比例して減少する負の項をもつ。この項を表面エネルギーと呼ぶ。これは表面張力と密接に関係する。表面上に仮想的に境界線を描いて、表面をその右側と左側にわけて考えるとき、右側と左側の引き合う力の大きさを境界線の単位長さ当たりで表せば 表面エネルギーに等し くなるのである。[10] 液滴全体の表面エネルギーは表面積に比例する。表面積 S は S = 4πR2で表せる。さらに R = r0 A1/3 なので、これを代入すると表面積 S は A2/3 に比例する
ことがわかる。ここで asuf を表面積に比例するパラメータとすると、Msurf = asuf A2/3

(2.3)

が得られる。

(3) クーロン項
原子核内の陽子間に働くクーロン力は、陽子同士が互いに反発し合うように働くので、結合エネルギーを更に減少させる。そこで、この減少分を全電荷 Ze が半径 RC = rC A1/3の球に一様に分布している時の静電エネルギーとして評価する。
電荷 Q が一様な体積密度 ρ で半径 RC の球面内に集まっているとき、静電エネルギー12

は次のように計算できる。徐々に集めてきて半径 r になったところを考え、これを drだけ増すときの仕事を dU とする。このときの表面における電位は ρ を電荷密度として 43 πρr 3 /4πε0 r で、球殻 dr の電荷は 4πρr2 dr であるから、

dU =4πρ2 4r dr
3ε0がよって全エネルギー U は
∫ RC

U=
0

4πρ2 RC54πρ2 4r dr =
3ε0
15ε0

(2.4)

4
3
Ze = πρRC3(2.5)となる。陽子の全電荷 Ze は、なので、(2.4) 式の両辺を (2.5) 式の両辺の２乗でそれぞれ割ってから両辺に (Ze)2 をかけると

1
3 Z 2 e2· ·U=
4πε0 5 RC(2.6)となる。RC = rC A1/3 なので、クーロンエネルギーの最終的な表式として

Z2
,
A1/3
3e2
aC =
20πε0 rC

MC = aC(2.7)

を得る。通常は (rC , RC , ρ でなく)aC をクーロン項のパラメータとする。[12](4) 対称項
質量数が小さな場合は原子核は同じ数の陽子と中性子を持つ傾向がある。しかし、原子核が重くなればクーロン斥力を核力で部分的に帳消しにするためにだんだん多くの中性子を含む。それにより中性子数と陽子数の非対称が生じる。しかし 、核子間の力は同種核子間よりも中性子-陽子間で π 中間子による相互作用が働くために強い力となる。また、同種の核子が多くなるとパウリの禁止則によって、後から核に取り込まれる核子の状態は他のものより高いエネルギー状態に入らなけらばならない。そこで、なるべく中性子数と陽子数の差を小さくした方がエネルギーは低い。これらの効果によるエネルギーを対称エネルギーと呼ぶ。ここで asym をこの項に比例するパラメータとすると、

Msym = asym

(N − Z)2.
A(5) 偶奇項13(2.8)

~ 、スピン角運動量の総和を S~ とする
原子核内の核子がもつ軌道角運動量の総和を L~ +S~ なる全角運動量を原子核のスピンという。実験によれば 、Z も N もとき、I~ = L
偶数の核のスピンは例外なく 0 である。この事実を説明するには、原子核の基底状態
では 2 陽子は必ず対になって合成角運動量が 0 になり、2 中性子についても同じことがいえると考えればよい。すなわち、同核子間では 2 個ずつが 0 の合成角運動量に組んでエネルギーを下げる傾向がある。それを対エネルギーまたは対ギャップと呼ぶ。これも核力の性質に起因するものである。質量数 A のべき乗 A−α に比例するとして実験値にフィット すると α ≈ 21 が得られる。

1/2


 aeo /A

Meo =

2.1.2





0
−aeo /A1/2

Z, N ともに偶数のとき (偶偶核)
A 奇数のとき (奇 A 核, 奇核)
Z, N ともに奇数のとき (奇奇核).

(2.9)

質量公式のパラメータの決定

Weizsäcker-Bethe の質量公式を陽子数 Z, 中性子数 N , 質量数 A の関数として表すと、
M (N, Z) = avol A + asurf A2/3 + aC


1/2


 aeo /AZ2
(N − Z)2
+
a
+
sym

A1/3
A



0
−aeo /A1/2(2.10)である。右辺の最後の項の場合分けは後述する ((2.11) 式を参照)。次にこの公式で、未知数のパラメータを求めるにあたって最小二乗法 (least squares method) を使い、前章で読み込んだ信頼できる原子質量実験データのある 2226 個全ての核種を用いて求める。また平均二乗誤差 (root-mean-square) からどの程度合わせれるかを見ることにする。
データのある核に、通し 番号をふり i 番目 (0 ≤ i < I = 2226) の核質量を f (i), 中性
子数を Ni , 陽子数を Zi , 質量数を Ai = Ni + Zi と表すものとする。

2/3

1/3

f0 (i) = Ai , f1 (i) = Ai , f2 (i) = Zi2 /Ai , f3 (i) = (Ni − Zi )2 /Ai ,

−1/2


 A

f4 (i) =





−1/2

−A
0

(Ni , Zi 共に偶数のとき)(Ni , Zi 共に奇数のとき)(上記以外のとき)(2.11)と定義すると質量公式は、

f (i) =K−1∑

ck fk (i)(2.12)k=0となる。ここで、ck は各項の未知定数であり (c0 = av , c1 = av , c2 = aC , c3 = asym , c4 =aeo )、K = 5 である。
i 番目 (0 ≤ i < I) の核質量の計算値を f (i), 実験値を yi とし f (i) と yi の誤差の平均

14

が最小となるパラメータ ck を知りたい。このために最小二乗法を用いる。最適の f (i)は、

E(c0 , c1 , c2 , c3 , c4 ) =I−1∑(yi − fk (i)) =
2

i=0

I−1
∑

(

i=0

yi −

K−1
∑

)2

ck fk (i)(2.13)(l = 1, 2, · · · , K − 1)

(2.14)

k=0を最小にする。したがって、
{

(

)}

I−1
K−1
∑
∑
∂E
= −2
fl (i) yi −ck fk (i)
∂cl
i=0
k=0

=0

を満足する ck の組を求めればよい。これから、最適値 c̃k に対して
K−1
∑
i=0

c̃k

(I−1
∑

)

fl (i)fk (i) =

i=0

I−1
∑

(l = 1, 2, · · · , K − 1)

yi fl (i)(2.15)

i=0が導かれる。

∑I−1

A
=

l,k
i=0 fl (i)fk (i)

Xk = c k∑
Bl = I−1
i=0 yi fl (i)(2.16)とおくと、Al,k 、Bl は与えられた I −1 個の値から直ちに計算できる。したがって (2.15)
式は
K−1
∑

Al,k Xk = Bl(l = 1, 2, · · · , K − 1)

(2.17)k=0となる。この連立 1 次方程式を数値計算によって解く方法として、ガウス法 (Gaussianelimination) を利用する [21]。この結果、求められたパラメータは (MeV 単位で )、av = −15.308
as =16.480

aC =

0.693

asym =

22.527(2.18)(2.19)

aeo = −11.245
となる。

2.1.3

液滴模型から得た値と実験値との比較

各パラメータの値を (2.10) 式の右辺の対応するそれぞれの項に代入する。そしてM (N, Z) を求め、計算値と実験値との誤差、つまり原子核質量の誤差を縦軸、中性子数と陽子数を横軸にとったものが図 2.2、2.3 である。2226 個のデータでの平均二乗
誤差は 3.41MeV となる。

15

10

Mexp-MWB

5

Mass Error (MeV)

0
-5
-10
-15
-20
-25
0

20

40

60

80
100
Neutron Number N

120

140

図 2.2: 原子質量の実験値と Weizsäcker-Bethe の質量公式の差と中性子数

10Mexp-MWB

5

Mass Error (MeV)

0
-5
-10
-15
-20
-25
0

20

40

60
Proton Number Z

80

図 2.3: 原子質量の実験値と Weizsäcker-Bethe の質量公式の差と陽子数

16

100

160

2.2

殻模型

図 2.2、2.3 を見るといくつかのところで急激な変化を見せ、誤差が負の方向へ落ち込む点が存在する。これは実験値が計算値より実際は低くなっているからであるが、なぜこういうことが起こるのか。この問題を解決したのが 1949 年、Mayar と Jensen
によって提唱された核の殻模型 (shell model) である。

2.2.1

魔法数の存在極めて安定な原子核の陽子数あるいは中性子数は特定の数になることは、経験的に知られている。このような数 (2、8、20、28、50、82、126) は魔法数 (magic number)と呼ばれる。陽子あるいは中性子が魔法数の原子核には、安定な原子核が非常に沢山ある。中性子数が魔法数であると、中性子 1 個を原子核から分離するのに大きなエネル
ギーを要する。これに対し、魔法数より 1 個だけ多いと、分離に要するエネルギーは小さい。これは陽子数についても同様である。つまり、これらの核が閉核 (closed shell)を形成し 、これを壊して励起状態をつくるには高いエネルギーを要する。
特に安定なのは二重魔法数 (double magic) の核、つまり陽子数も中性子数も魔法数の原子核である [5]。

2.2.2スピン軌道結合力の導入

では、魔法数がなぜ存在するのか。それは殻模型 (shell model) によってうまく説明ができる。そのためにまず殻模型に適したポテンシャルを導入しなければならない。
強い相互作用、つまり核力の到達距離は非常に短いので、ポテンシャルの形は原子核の中の核子の密度分布に比例するはずである。このポテンシャルを近似して、数学的に取り扱いやすい調和振動子 (harmonic oscillator) で近似する。調和振動子ポテン
シャルを

1
V (r) = M ω 2 r2
2

(2.20)

ここで、M は核子の質量、 ω は調和振動子の角速度である、r は核子の半径である。1つの核子が , 等式 (2.20) という平均ポテンシャル内を運動しているとすると、この核子に対する Schrödinger 方程式は
)

(

1
h̄
∇ + M ω 2 r2 φ(r) = Eφ(r)−2m
2(2.21)

この上式を解くと、エネルギー固有値 E として()

3
E = h̄ω N +,2
N = nx + ny + nz

(2.22)(2.23)

ここで nx 、ny 、nz は x 軸、y 軸、z 軸についてのそれぞれの量子数であり、N はそれらを総和した量子数である。
17

一番エネルギーの低い状態は、量子数が N = 0 つまり (nx , ny , nz ) = (0, 0, 0) で E =(3/2)h̄ω の準位エネルギーをもつ。次に N = 1 では、(nx , ny , nz ) = (1, 0, 0), (0, 1, 0), (0, 0, 1)という 3 つの状態が縮退している。そして、パウリの原理にしたがって、1 つの準位に
2 つの中性子 (または陽子) までとし 下の準位から準に詰めていくと、閉殻数は表 2.1のようになる。 表 2.1 の閉殻数を見ると、N = 3 までは魔法数と同じであるがそれ
N

0

1

2

3

4

5

6

···

縮退数

2

6

12

20

30

42

56

···

閉殻数

2

8

20

40

70

112

168

···

表 2.1: 調和振動子型での量子数 N に対する縮退数と閉殻数

以降は魔法数と異っている。そこで、核子密度分布をより考慮したポテンシャルとして Woods-Saxon ポテンシャルを用いる。

V (r) =

−V0
1 + exp (r − R)/a(2.24)

ここで、R = r0 A1/3 は核半径、a は核表面のぼやけを表すパラメータである。このポテンシャルを加えた Schrödinger 方程式を解くと、調和振動子の場合は主量子数 n と軌道角運動量量子数 l が異なっていても N が同じであれば縮退している。しかし、WoodsSaxon ポテンシャルの場合では N が同じでも nl が異なれば 縮退していない。それではなぜ違いが生じるのか。
そこで N = 2 のときを見てる。調和振動子では 2s と 1d が縮退していることが、各量子数によって明らかになる。しかし、Woods-Saxon ポテンシャルの場合では図 2.5(b)のようにわかれている。それは図??を見ると半径 r が 0 に近い場合には Woods-Saxon
ポテンシャルエネルギーと調和振動子のポテンシャルエネルギーはさほど変わりは無
いが、調和振動子のポテンシャルエネルギーに比べて半径 r が核半径付近では低くなっている。例えば、2s は原子核の中心近くで存在する確率が高く、一方 1d は核表面付近に存在する確率が高いことになる。つまり、角運動量が小さい 2s は核中心部での存在確率が高いために、エネルギーの変化はそんなに変わらない。一方、角運動量が大き
い 1d は核表面付近に高く存在するためエネルギーは下がる。その結果、図 2.5(b) のように 2s と 1d は分かれるのである。

18

Potential Energy V

Harmonic Oscillator
Woods-Saxon

Radius r
図 2.4: 調和振動子ポテンシャルと Wood-Saxon ポテンシャル

これにより、Woods-Saxon ポテンシャル閉殻数を見てみると表 (2.2) のようになる。 この場合でも魔法数 2、8、20 は説明できても、それ以降の魔法数について説明
N

0

1

2

2

3

3

4

4

4

···

nl

1s

1p

1d

2s

1f

2p

1g

2d

3s

···

縮退数

2

6

10

2

14

6

18

10

2

···

閉殻数

2

8

18

20

34

40

58

68

70

···

表 2.2: Wood-Saxon 型での量子数 N に対する軌道名 nl と量子数縮退数と閉殻数できない。そこで、さらに核子の内部角運動量であるスピンと軌道角運動量を結合させるスピン軌道力を考慮に入れることによって全ての魔法数を説明できることが発見された。
核子の軌道角運動量 l とスピン s の結合による合成角運動量 j は
 j =l± 1
 j = 122

(l 6= 0)(l = 0)(2.25)となる。スピン軌道力の κ を正の定数として、−κ~l · ~sという形を仮定すると}{~l · ~s = 1 (~l + ~s)2 − ~l2 − ~s2219(2.26)

{

(

)}

1
1 1
(j 2 − l(l + 1) −
+1
2
2 2


l
(j = l + 12 のとき)=
 −(l + 1) (j = l − 1 のとき)
2=

(2.27)なので、合成角運動量が大きいとき、すなわち j = l + 21 の場合、エネルギーは低くなる。合成角運動量が小さいとき、すなわち j = l − 12 の場合は、エネルギーは高くなる。
つまり、j = l + 21 のレベルは常に j = l − 12 のレベルより下にある。このスピン -軌道結合を加えたことにより、閉殻数と魔法数との一致する。そのエネルギー準位を表したものが図 2.5(c) である。

20hω/2π
5

N=5
p

126

f
h

2p1/2
1f5/2
2p3/2
0i13/2
0h9/2
1f7/2

82
N=4

s
d

4

1d3/2
2s1/2
0h11/2
0g7/2
1d5/2

g
50
N=3

p

3
f

28

0g9/2
1p1/2
0f5/2
1p3/2
0f7/2

20
s

N=2

0d3/2
1s1/2

2
d

0d5/2

8
N=1

0p1/2

p

1

0p3/2

2
N=0
0

s

(a)

(b)

0s1/2

(c)

図 2.5: (a) 調和振動子ポテンシャルによるエネルギー準位、(b) Wood-Saxon ポテンシャルによるエネ
ルギー準位、(c) スピン -軌道結合を加えた場合のエネルギー準位

21

第3章

3.1

変形を考慮した質量公式の作成

変形核と四極子モーメント

原子核を球対称として記述する殻模型近似は、スピン -軌道相互作用を付け加えれば 、閉殻の二重魔法数の近傍での原子核についてはよく成り立つ。一方、殻の途中まで埋まっている原子核、つまり合成軌道角運動量 j に核子が全て詰まっていない原子
核は変形しておりポテンシャルは球対称ではないので球対称の殻模型はよい近似ではない。では、なぜ原子核は変形するのか。
球状が球形からずれて歪んでいる原子核を変形核と呼ぶ。最も一般的な変形は四極型変形である。特に希土類元素の原子核は非常に大きな電気四極子モーメント の値をもつ。そこで、核の変形と四極子モーメント Q との関係について明らかにする。
殻模型では閉殻では球対称のため四極子モーメント Q をもたず、Q = 0 となる。次に殻模型の単一陽子状態における電気四極子モーメント Qj を求めると

Qj = −

2j − 1< r2 >
2j + 2(3.1)

である。ここで < r2 > は r2 の平均値である。球対称の閉殻の外を 1 個の陽子が回っている四極子モーメント Qj は、(3.1) 式により Qj < 0 となる。これは m = j の状態における期待値であるから、陽子は核の赤道面を回っており、したがって原子核全体が赤道面方向にふくれてみえるわけである (oblate 変形)。一方、もし 1 個の陽子が付け加わると球対称な閉殻になる場合 (単一陽子空孔状態) は、赤道面を回る陽子が 1 個取り除かれたわけであるから、核は両極の方向に延びたようにと見える (prolate 変形)。
したがってめ四極子モーメント は Q−1 = −Qj > 0 となる。(図 3.1 参照)。
< r2 > のだいたいの値を知るために、半径 R = r0 A1/3 の球内で一様で、外では 0と仮定すると3< r2 >= r02 A2/3
5

(3.2)

となるから、

Qj = −

2j − 1 3 2 2/3
r A
2j + 2 5 0(3.3)

を得る。例えば 、17571 Lu では、j = 7/2 であり、r0 = 1.2fm とすると
−28 2
Q−1
m
j = 0.19 × 10

(3.4)

となる。一方、実験値は Qobs = 5.68 × 10−28 m2 であるから、Qobs /Q−1 = 30 となり、陽子の単一粒子模型ではこの大きな Qobs の値を全く説明できない。
22

図 3.1: (a)prolate 変形, (b)oblate 変形

このことから、J.Rainwater は一般の原子核の安定な形は必ずしも球形ではなく、球形からわずかにずれた回転楕円体の形になっている、すなわち、変形していると考え、この矛盾を解決した [11]。

3.2

変形核の殻エネルギー

変形した原子核を論じるにあたっての正攻法では、変形したポテンシャル中の一粒子の固有状態を求めることが必要である。しかし 、それを数値計算で求めるのに要す
る計算機資源の分量（ 計算時間、記憶容量など ）は、変形していない球対称なポテン
シャル中の固有状態を求める場合に較べて桁違いに多くなる。個々の原子核を論じるのでなく、原子核の質量公式を論じる場合には、計算量がさらに２〜３桁大きくなる。
なぜなら質量公式がすでに決定済みの場合でさえ、公式と実験データとの平均誤差を計算するには核図表一面の計算（ 即ち、数千個の原子核についての計算）を実行しなければならないからである。ポテンシャルが変形している場合はこの核図表一面の計算だけでも既に（ 最新の計算機をもってしても ）重い計算である。そして、本論の目標である質量公式を（即ち、その 15〜30 個程度のパラメータの値を）決定し 、改良するという目的のためには、この計算を何千回、何万回、と繰り返さなければならない。
最小化の対象とする関数（ 質量の平均２乗誤差）の高次元のパラメータ空間における挙動は単純なものではないと予想されるので、おそらくは最小値への接近は遅々としており、計算機資源の利用できる限りの計算を続行したいという状況になるであろう。
そのような状況下においては、個々の原子核を計算する速度をできる限り短く抑える
ことが決定的に重要となる。
そこで、もし、球対称なポテンシャル中の固有状態に関する情報のみに基づいて、原子核が変形した場合のエネルギーの変化を求める近似法が存在するなら、大変有益であると言えよう。そのような方法のひとつが、2000 年に小浦らによって KUTY 質量公式の重要な要素として発表された [19]。小浦らの KUTY 質量公式は、核質量の平23均二乗誤差が 0.715MeV というまずまずの成果を挙げていることから判断して、彼らによる変形の近似的取り扱い法を、KUTY 公式から切り離し、他の手法との組み合わせで試してみる価値があると我々は考えた。
そこで、まず、この節では、小浦らが導入した「変形核を球形核の重畳と見る近似」
を解説する。
なお、小浦らの方法では、「変形度」と「核図表上で近傍にある原子核の殻エネルギーの混合の仕方」とを密接に関連づけるのであるが、その根拠として小浦らが主張することは難解であり万人が理解し 認めるものとは言えないと思われる。しかし 、この方法には、平均場の変形を正確に考慮するという正攻法に優る現象論としての強みも逆に期待できるであろう。例えば 静的変形を超えた動的揺動の効果や複雑な相関の効果も同時に取り入れ得る可能性が期待される。そうであればこそ、核質量の平均二乗誤差を 0.715MeV にまで低減することができたのであろうと我々は推測する。

3.2.1

球形核の重ね合わせとしての変形核

変形核の固有殻エネルギーを得るため、幾つかの仮定を用いる。まず、変形核は適当な球形核の重ね合わせとして表せるものとする。その重ね合わせについて、干渉を伴わないものと仮定する。変形核状態は原則として、様々な陽子数と中性子数を持つ適当な球形核の最低エネルギーでの占有率の重ね合わせによって近似させる。次に、変形核の固有殻エネルギーは占有率と同じ重ね合わせの重みで球形殻エネルギーの重ね合わせとして表されることができるとする。
変形核状態を球対称の重畳として表すときの占有率は、変形核の (３次元空間における) 密度分布を、球形核の密度分布の重畳で表すための占有率として求めることにする。後者を決定するため、以下のような仮定をおく。変形核に対する表面の内部を一様な密度と仮定する。核の変形を比較的簡易なものに制限する。つまり、軸対称かつ鏡映対称な球体をとり、中性子と陽子のグループについて同様な形状であるとする。核の表面は sharp-cut であるものとし 、中性子と陽子の分布は一様であるとする。したがって、核形状は以下に述べる極角 (polar angle)θ の関数として記述される。

r(θ) =

R0[1 + α2 P2 (cos θ) + α4 P4 (cos θ) + α6 P6 (cos θ) + · · ·]
λ

(3.5)

P2i (cos θ)(i = 1, 2, ...) は 2i 次の Legendre 多項式であり、α2i は形状を指定するためのパラメータである。θ の範囲は 0 ≤ θ ≤ π/2 とする。R0 は球形核の半径は、
1

R0 (Z, N ) = r0 (N + Z) 3

(3.6)

と表せるとし 、r0 = 1.2(fm) とする。体積が考慮中の核の球形での体積と等し くなり、
[

λ= 1+

3
4π

(

4π 2 4π 2 4π 2α +α +α + ···5 2
9 4 13 6

)]1/3

(3.7)

となり、P6 (cos θ) 項まで等式 (3.5) 内の展開をとる。
重ね合わせの重み W (N, Z) の意味であるが、球形核の密度分布から余分な所を削り

24

W(N1)

図 3.2: 重み W のイメージ図

落として各準位の所がむきだしになった立体角を表し 、その面積の総和を 1 とするものと仮定する (図 3.2 参照)。
変形核の表面の内部密度は一様であると仮定するために、球形核の内部密度も一様である必要がある。そのため、r から rim に半径座標を変換する。そして、あたかも全内部密度が一様に成るように問題を扱うことの出来る、中間形状 (imtermediary shape)
へ核形状を対応させる。
変形核が与えられるとき、占有立体角 (occupied solid angle)Ωoc を下式定義する。
∫ πΩoc (r0 ) =0dθ · sin θ∫ 2π0dϕΘ(r (θ, ϕ) − r0 )(3.8)

ここで θ は極角、ϕ は方位角、Θ は階段関数で
 1 (x > 0)Θ=
 0 (x < 0)(3.9)

である。つまり、半径 r の球体の外側に広がる核表面部分が占有する立体角が Ωoc である。軸対称形について、占有立体角 Ωoc を計算する際には、prolate 変形 (葉巻型)と oblate 変形 (パンケーキ型) の２種類の変形に分けて考える必要がある。下記は、Prolate 形については至る所で dr< 0 である場合のみを考え、Oblate 形についてはdθ
至る所で dr> 0 である場合のみを考えて得られた結果である。なお、本論文の執筆時
dθ
点の研究段階では、形状パラメータのうち α2 のみをゼロでないとして考慮するため、
上記のどちらかの場合しか起き得ず、形状が 0 < θ < π2 で極値を持つ場合を考慮する必要はない。
∫ θΩoc (r(θ)) = 4π

sin θdθ0= 4π(1 − cos θ)∫ πΩoc (r(θ)) = 4π

2(prolate 変形)(3.10)(oblate 変形)(3.11)sin θdθθ= 4π cos θ

さて、以降のシナリオを簡単化してに述べておく。数式による正確な定義は後述する。
25

Ω oc(r)

r(θ)

r(θ)
r
r(π−θ)(0<θ<π/2)= =図 3.3: 占有立体角 Ωoc (r )1. r = r0 A1/3 として、r を核子の個数に対応させる。
2. 中性子と陽子の個数比は考察下の原子核と同一とする。
3. 変形により r が一定値をとるのでなく有限の幅に分布することを、A が一定値でなく有限の幅に分布しているのだと解釈する。
4. このときの r の分布の（ 即ち A の分布の) 密度関数として表面の立体角を使う。理由は殻エネルギーが表面効果であるからであろう。
実は、小浦らは、「中間形状」を導入することで、彼らの処方をもう一段複雑化させた。
「中間形状」の導入が必要な理由は、恐らく、実際の形状をそのまま使用したのでは質量公式の精度を上げ ることに困難があったのであろうと我々は想像している。中間形状の定義には、物理的意味の不明瞭な調節可能なパラメータ (後述の h) が １個含まれており、このパラメータの自由度を利用することで、精度向上に立ちはだかった
困難を克服することができたのであろう。
まず、中間形状を表すため、半径 r から新たな半径変数 rim への写像を下式で定義
する。
{
[
]}
Ωoc (r )drim = 1 − h 1 −dr(3.12)4π
ここで、h は正符号で１以下の調節可能なパラメータである (1 ≤ h ≤ 1). また、r = 0
において rim = 0 とする。
この新変数の意味を考察しておこう。変形核の最小半径を Rmin と表せば 、r ≤ Rmin
では、Ωoc (r ) = 4π ゆえに drim = dr であるから、rim = r である。また、r > Rmin においては rim ≤ r である。従って、変形に伴う半径の分布の幅を小さくするために rimを導入したのだと推測される。
次に、この新変数を使って核の中間形状（核の重畳の確率を計算する用途にのみ使用
26

r

rim

図 3.4: rim > r1 での方向の占有立体角 Ωim (r )される仮想的な形状）を下式で定義する (図 3.4 参照)。
2
Ωim (rim )rim
drim = Ωoc (r)r2 dr(3.13)

ここで、Ωim は中間形状において、半径が引数（ 上式では rim ）より大きい方向の占有
立体角を表す (図)。定義 (3.13) の意味は、半径が [r, r + dr] の範囲に含まれる物質は、中間形状の核では新変数が [rim , rim + drim ] の範囲にちょうど移動できるということである。また、その結果として、核の中間形状とオリジナルの変形核とは密度の積分値（ 核子の総数）が同一になる。
(3.12) 式と (3.13) 式により、与えられた形状 r = r(θ) の r ないし θ の各値に対応する中間形状の rim と θim の２値が定まるのである。例えば prolate 形の場合は、下記のような関係が成り立つ。
• (r = Rmin , θ = π2 ) に (rim = Rmin , θim = π2 ) が対応する。
• (r = Rmax , θ = 0) に (rim < Rmax , θim = 0) が対応する。
• θ = θim = 0, π2 以外では、rim ≥ r, θim ≥ θ である。
以下の式変形を追うにあたって、重要な留意事項は、４変数 r, θ, rim , θim は従属関係にあり、一つの変数を決めると、残りの３変数の値が決まってしまうこと、また、一般には r 6= rim 、θim 6= θ ということである。
なお θim は以下の議論では使用しないが、形状を図示したいときには必要である。表√面の座標 (x, y, z) は、本当の形状では、z = r cos θ 、 x2 + y 2 = r sin θ で与えられる√が、中間形状の表面の座標は、z = rim cos θim 、 x2 + y 2 = rim sin θim で与えられるの
27

である。
(3.12) 式に (3.10)、(3.11) 式を代入するとdrim =

 (1 − h cos θ)dr

(prolate 変形)
 1 − h(1 − cos θ)dr (oblate 変形)(3.14)となる。r : rim = Rmin に対応する θ = θim の値を θmin と書くことにする。すなわち、

Rmin = rim (θmin )


 θ

π
min = 2(prolate 変形) θmin = 0 (oblate 変形)(3.15)である。(3.15) 式の両辺を θmin から θ までの範囲で変数 θ で積分すると
∫ θ

∫ θ
θmin

drim
dθ
θmin dθ∫ θdrim dr=dθθmin dr dθdrim =(3.16)となる。右辺に現れる drim /dr は (3.15) 式より、dr /dθ は (3.5) 式より求まる。また、左辺は下記のように書くことも出来る。
∫ θ
θmin

drim = rim (θ) − rim (θmin ).

(3.17)(3.17) 式の rim (θmin ) は Rmin に等しい。(3.16)、(3.17)、(3.15) 式をまとめると∫ θ

rim (θ) = Rmin +drim drdθθmin dr dθ(3.18)(im)となる。この (3.18) 式から rim の最大値 Rmaxが求まる。
中間形状の半径変数 rim に対応する原子核の中性子数 N1 、陽子数 Z1 は(rim
N1 (rim ) =
R0

)3

(

rim
N, Z1 (rim ) =
R0

)3

Z

(3.19)

で与えられるものとする。逆変換は下記のとおりである。
(

N1
rim (N1 ) =
N

)1/3

(

Z1
R0 , rim (Z1 ) =
Z

)1/3

R0

(3.20)

この等式 (3.19) によって、rim の最大値と最小値から重ね合わせで使われる原子核の N1 と Z1 の範囲が決まる。ここで最小の N1 を N1min とし 、最大の N1 を N1max と表すことにする。道様に最小の Z1 と Z1min 、最大の Z1 と Z1max と表すことにする。θ のもつ範囲からそれそれ変数の最大値と最小値が決まり、それらをまとめたのが表 (3.1)
である。また、図 (3.5) と図 (3.5) は , 中性子数 N = 98 と陽子数 Z = 68 のエルビウム（ Er ）で α2 = 0.2 と設定した場合の r と rim の関係表したものであり、r 、rim と θ の極座標を表したものである。

(3.18) 式で計算できる θ と rim の関数関係の逆関数、
θ = θ(rim )28(3.21)

prolate 変形

oblate 変形

θ

0

π /2

0

π /2

r

Rmin

Rmin

Rmax

rim

Rmax
(im)
Rmax (< Rmax )

Rmin

Rmin

Rmax

Ωoc

4π

0

0

4π

Ωim

4π

0

0

4π

N1

N1max

N1min

N1min

N1max

(im)

表 3.1: 2 種の変形の、θ が最小と最大のときの各関数の値域

7.5

Radius rim (fm)

(im)

Rmax

6.5

Rmin

5.5
5.5

Rmin

6.5

Radius r (fm)

7.5

Rmax

図 3.5: エルビウム (N = 98, Z = 68, α2 = 0.2) での、半径座標 r 、rim の最小値 Rmin と最大値 Rmax

29

は Lagrange 補間多項式により導き出す。
ここで Lagrange 補間多項式について軽くふれておく [20]。
n 次式の Lagrange の補間の場合、相異なる n + 1 個の点
x 0 , x1 , x2 , . . . , x n
でそれぞれの値
f0 , f1 , f2 , . . . , fn
をとる n 次多項式 f (x) はf (x) =

n∑Ln,k (x)fk ,(3.22)

k=0(x − x0 )(x − x1 ) . . . (x − xk−1 )(x − xk+1 ) . . . (x − xn )(xk − x0 )(xk − x1 ) . . . (xk − xk−1 )(xk − xk+1 ) . . . (xk − xn )n∏x − xi.
(3.23)=
i=0 xk − xi

Ln,k (x) =i6=k

この Lagrange 補間多項式を用いることにより、rim から θ を導き出すことができる。
以上の準備により N1 → rim → θ → r の順に求める事が可能となる。
次に中間形状の占有立体角 Ωim (rim (N1 )) を等式 (3.10)、(3.11)、(3.13)、(3.14) から導き出すと( )2
 4π(1−cos θ) r
1−h cos θ(rim )2Ωim (rim (N1 )) =
4π cos θr
1−h[1−cos θ]

rim(prolate 変形)(3.24)(oblate 変形)となる。核 (N1 , Z1 = ZN1 /N ) の重み W (N1 ; N, Z) はW (N1 ; N, Z) = −1 dΩim (rim (N1 ))4πdN1と表される。N1min から N1max までの重み W の総和を W̄ とすると、

W̄ =

N∑
1max

Wn

N1 =N1min∫ N1max=
N1min∫ N1max

=
N1min

=
=
=
=

Wn (N1 )dN1
1
−
4π

(

)

dΩim (rim (N1 ))dN1
dN1

1 ∫ N1max−
dΩim (rim (N1 ))4π N1min1
1max− [Ωim (rim (N1 ))]N
N1min
4π
1
[Ωim (rim (N1min )) − Ωim (rim (N1max ))]
4π
]
1 [
im
Ω(im) (Rmin ) − Ωim (Rmax)
4π30(3.25)

1(4π − 0)4π
= 1=(3.26)となり、正し く規格化できていることが確認できる。
中性子数が整数であることを考慮して、中性子数が N1 の核に対する重みを求めると、
N1min + 1/2 ≤ N1 ≤ N1max − 1/2 のときは
∫ N1 + 1

W =2

N1 − 12W (N1 )dN1

[

1
1
1
=
Ωim (rim (N1 − )) − Ωim (rim (N1 + ))
4π
2
2

]

(3.27)となる。しかし 、陽子と中性子数の比を保つという条件により、中性子数 N1 の核の陽子数 Z1 は

Z1 =

N1
Z
N

(3.28)となり、N1 を整数として扱うことにしても、一般には Z1 は実数になってしまう。
整数でない陽子数・中性子数に対する平均場の解を定義することに原理的な問題はない。波動関数として BCS 変分関数の対相関ギャップがゼロの極限を用いて（ 別の言い方をすれば、偶々核に対する密度汎関数を用いて）、密度の積分値が指定された実数値になるという拘束条件のもとに変分を実行すればよい。
しかし 、∆N1 = ∆Z1 = 0.1 で核図表の計算を行うと、現実の 100 倍の個数の原子核を計算することになり計算が過大になる。そこで、実際の計算では、整数の N1 , Z1に対してのみ平均場計算を行い、その計算結果から、実数の N1 ,Z1 へ内挿・外挿を行うことにする。ここで、平均場解のエネルギーが N1 ,Z1 に対して滑らかな関数になる
ように、直前の段落で述べたとおり、奇数の核子数に対しても平均場計算では偶々核に対する波動関数（あるいは密度汎関数）を用いる。
内挿・外挿の方法は具体的には、以下に述べる通りに行う。まず、N1 が整数で Z1 が非整数の点（ 核）への内挿・外挿は２点を通る一次式で行う。つまり、E(N1 , Z1 ) =Z3 − Z1
Z1 − Z2
E(N1 , Z2 ) +E(N1 , Z3 )Z3 − Z2
Z3 − Z2(3.29)とする。ここで、Z2 ,Z3 は整数であり、(Z2 ,N1 ),(Z3 ,N1 ) に平均場解が存在するものを選ぶ。Z2 < Z1 < Z3 を満たす Z2 ,Z3 がとれる場合は、もっとも短い間隔での内挿を行う。通常は Z3 = Z2 + 1 ととることができる。とれない場合は、Z1 に最も近い Z2 ,Z3を使って外挿を行う。
Z1 が整数で N1 が非整数の点（ 核）への内挿・外挿も同様に下式で行う。

E(N1 , Z1 ) =N3 − N1
N1 − N2
E(N2 , Z1 ) +E(N3 , Z1 )N3 − N2
N3 − N2(3.30)

N1 と Z1 が共に非整数の点（核）への内挿・外挿をするには、まず、N2 < N1 < N3 、Z2 < Z1 < Z3 を満たす 4 整数の組 N2 , N3 , Z2 , Z3 で、4 点（ 核）(N2 ,Z2 ), (N2 ,Z3 ),(N3 ,Z2 ), (N3 ,Z3 ) に平均場解が存在するものを選ぶ。通常は N3 = N2 +1、Z3 = Z2 +131

を満たすようにとれる。そして、これらの４点（ 核）で平均場解のエネルギーに一致する下記の多項式を使って内挿する。

E(N1 , Z1 ) =1{(N1 − N2 )(Z1 − Z2 )E(N2 , Z2 )(N3 − N2 )(Z3 − Z2 )+(N1 − N2 )(Z3 − Z1 )E(N2 , Z3 )+(N3 − N1 )(Z1 − Z2 )E(N3 , Z2 )+(N3 − N1 )(Z3 − Z1 )E(N3 , Z3 )}(3.31)

上式は、N1 , Z1 に関する二次式であるが、N12 , Z12 に比例する項を除外してあるため、
N 軸方向、Z 軸方向には直線として振舞う。このため、片方が整数の場合は (3.29) 式
または (3.30) 式と一致する。
さて、一個の変形核を表すために重ね合わせるべき核種は、変形が大きくなれば 核
図表上で非常に遠くまで延び得る。このために、計算中に平均場解が存在しない核のエネルギーが必要となったために計算が異常終了してしまう可能性がある。これを回避するため、できれば 全く任意の核種の球形殻エネルギーを定義しておくことが望ましい。
そのために、次のような手順で、任意の (N1 , Z1 ) に対するエネルギーの値を定義す
ることにする。
まず、与えられた平均場模型のハミルト ニアンの、平均場解の存在する核の全てについてエネルギーを求める。解の存在する核の中性子数・陽子数の最小値・最大値を
mf
mf
mf
mf
Nmin
, Nmax, Zmin, Zmax
とする。
mf
mf
mf
mf
外挿の第一段階として、Nmin≤ N ≤ Nmax
かつ Zmin≤ Z ≤ Zmaxの矩形領域にある核で平均場解の存在しないものについて下記に述べる一次式による外挿を繰り返して、この矩形領域にある全点（ 核）のエネルギーを定義してしまう。
外挿の具体的な手順は下記のとおりである。外挿したい点を (N1 , Z1 ) とし 、それを1 頂点とする一辺 1 の正方形を作り、残りの 3 頂点でエネルギー値が定義済みならば 、３頂点の値で決まる N1 および Z1 の一次式を使って外挿を行い、その点のエネルギーを決める。
例えば 、図 3.6 の場合について外挿値を数式で表せば下記のとおりとなる。

E(N1 , Z1 ) = E(N1 + 1, Z1 ) + E(N1 , Z1 + 1) − E(N1 + 1, Z1 + 1)

(3.32)

このようにして外挿できる点がなくなれば 、次は、外挿で定義された点もエネルギーが定義済みであるとみなして、同じことを繰り返す。最終的には矩形領域の全点のエネルギーが定義されることになる。
外挿の第二段階として、第一段階で外挿値による充填が完了した矩形領域の外にある点については、(3.31) 式で N2 、N3 、Z2 、Z3 を下記のとおりにとって外挿を行う。
N3 = N2 + 1、Z3 = Z2 + 1 とする。
mf
mf
mf
mf
とする。
の場合は、N2 < N1 < N2 + 1, Z3 = Zmax
, Z1 > Zmax
≤ N1 ≤ Nmax
• Nminmf
mf
mf
mf• Nmin ≤ N1 ≤ Nmax , Z1 < Zmin の場合は、N2 < N1 < N2 + 1, Z2 = Zmin とする。
mf
mf
mf
mf
とする。
の場合は、Z2 < Z1 < Z2 + 1, N3 = Nmax
, N1 > Nmax
≤ Z1 ≤ Zmax
• Zmin32

Number of possible calculations of extrapolation: n=1+Z
E(N-1,Z+1)

E(N,Z+1)

E(N+1,Z+1)

E(N-1,Z)

E(N,Z)

E(N+1,Z)

E(N-1,Z-1)

E(N,Z-1)

E(N+1,Z-1)+N

Unknown Data
Known Data
Value from extrapolation
図 3.6: エネルギーデータの無い球形核種への外挿

33

mf
mf
mf
mf
• Zmin≤ Z1 ≤ Zmax
, N1 < Nminの場合は、Z2 < Z1 < Z2 + 1, N2 = Nmin
とする。
mf
mf
mf
mf
• N1 > Nmax Z1 > Zmax の場合は、N3 = Nmax Z3 = Zmax とする。
mf
mf
mf
mf
• N1 > Nmax
Z1 < Zminの場合は、N3 = Nmax
Z2 = Zmin
とする。
mf
mf
mf
mf
• N1 < Nmin Z1 > Zmax の場合は、N2 = Nmin Z3 = Zmax とする。
mf
mf
mf
mf
• N1 < Nmin
Z1 < Zminの場合は、N2 = Nmin
Z2 = Zminとする。
上述のような外挿の手順を採用する最大の理由は、結果として得られる E(N1 , Z1 ) が至る所で連続な関数になることである。もし 連続でなければ 、一個の原子核のエネル
ギーが変形度を変化させたときに不連続に変化する現象が起き、物理的に承服しがた
い模型となってしまうのである。
最後に重ね合わせの重みの再規格化について述べる。重みを N1 の最小値から最大値まで積分すれば 1 になるはずであるが、数値計算誤差のため 1 からのずれが生じる。このずれを補正するため、エネルギーを積分するために用いるのと同一の数値積分法 (例えば 本論文の場合は ∆N1 = 0.1 での台形則) で１を積分した結果で積分結果を割っておくのが望ましい。即ち下記のように置き換えるとよい。

∫ N1max∫ N1max
N1min

E(N1 )Wn (N1 )dN1

−→

N1min

E (N1 )Wn (N1 )dN1
∫ N1max
N1min(3.33)
dN1

N1 を整数に限り、有限の和で表す場合は、下記の置き換えで再規格化することになる。
N∑
1max
N∑
1max

E(N1 )Wn (N1 )−→N1 =N1min

E(N1 )Wn (N1 )N1 =N1minN∑
1max(3.34)Wn (N1 )N1 =N1min(3.34) 式で表される再規格化が不可欠な場合があり、それは、殻エネルギーを積分でなく整数についての有限和で計算するときに、
1 ∫ N1 + 2 dΩim
dN1
4π N1 − 21 dN1
1Wn (N1 ) = −とせず、(3.35)

1 dΩim(3.36)
4π dN1∑とする場合である。この場合 N1 Wn (N1 ) の 1 からのずれは 0.1 を越えることが頻繁に起きる。この大きなずれの原因は次に論じる Wn (N1 ) の振る舞いにある。
図 3.7 は N = 98, Z = 68 のエルビウム [Er] 同位体が Prolate 変形 (α2 = 0.2) および Oblate 変形 (α2 = −0.2) したことにより周辺の原子核（ 中性子数は N1 、陽子数は68N1) が混ざる重み W (N1 ) を N1 を横軸にとって描いたものである。ここで、N1 の値
98は整数に限り、各整数値に対する重みは (3.35) 式のように N1 ± 21 を端点とする区間の積分により得たものである。ただし 、積分は N1min および N1min で打ち切るため、最大および最小の N1 に対しては積分区間の幅が１より小さくなるので、重みは隣の点列からはずれた小さい値となる。
Wn (N1 ) = −34

球形 (α2 = 0) では N1 = 98 に 1 の重みが集中しているのに対し 、図 (3.7) を見ると、|α2 | = 0.2 という控えめな大きさの変形でも、プロレート 変形で 60, オブレート 変形で60 もの幅が中性子の個数にある原子核が混ざってくることがわかる。この分布の端に閉殻核を捉えるまで変形すれば殻エネルギーを大きく稼ぐことができるのであろう。
もし 、モデルに中間形状を導入しなければ (h = 1 ならば ) 混ざりの幅は中性子の個数にして数十個分にのぼり、閉殻を越えてまだ先へ延びるであろう。そのため、エネ
ルギー最小化で決定される変形度 |α2 | は、0.2 でさえ過大で 0.02 程度にしかならないであろう。これは原子核の実際の変形度 (0.3 程度) に合致しない。したがって中間形状の導入は必須であり、また、小浦らの h = 0.46 ∼ 0.5 という設定値にも大きくする
方向への変更の余地はあまりないのだろうと思われる。
重みの分布の形状を見ると、プロレート 変形の場合は N1min に、オブレート 変形の場
合は N1max に急峻なピークができている。このために、この模型では、下の閉殻核に近い核はプロレート 変形しやすく、上の閉殻核に近い核はオブレート 変形しやすいと予想される。この傾向は変形の微視的な理解と概ね一致している。
ただし 、これではオブレート 変形した核とプロレート 変形した核が同程度の個数存
在することになりそうであり、原子核のプロレート 変形優勢という事実が再現できないであろう。後述するように、これを解決するために、小浦らはプロレート 変形優勢の傾向をエネルギーの巨視的部分に付与した。我々もまたこの論文では小浦らにならう。しかし 、プロレート 変形優勢の起源はむし ろ微視的エネルギーにあると考えられるので、将来的には、微視的部分にプロレート 優勢の傾向を付与することが望ましい。
これは、変形核を球形核の重ね合わせにおきかえる処方という模型の中心部分の改変を意味する。
核の変形を与えるパラメータ α2 に制限を加える。なぜなら、核を非常に軟らかい場
合にはどこどこまでも変形しうるからである。しかし 閉殻内での変形と制限されるため α2 は、

α2 = 0.6 −A
1000(3.37)とする。この制限について、軽核では次の閉殻まで広い空間が確保できるため、α2 の範囲も広いと考えられるが、重核の場合は次の閉殻まで空間的余裕が軽核と比べて狭いと考えるために、α2 の範囲も狭くなると考えられるからである。

35

0.2

Weight W(N1)

0.18
0.16

α2= 0.2 (Prolate Shape)
α2=−0.2 (Oblate Shape)

0.14
0.12
0.1
0.08
0.06
0.04
0.02
0
40

50

60

70

80

90

100

110

120

Neutron Number N1• Prolate 変形と Oblate 変形では、傾き方が反対になる。
• 変形度を表すパラメータ α2 により、N1 の範囲の決まり方が決まる。
図 3.7: エルビウム [Er](N = 98, Z = 68) の中性子数 N1 と重み W (N1 ) の関係

36

130

3.2.2

平均変形エネルギー

変形に伴い、微視的な固有殻エネルギーだけでなく、巨視的な (液滴としての) 平
均変形エネルギーも変化することを考慮しなければならない。両者を加え合わせたものが全エネルギーだからである。
平均エネルギーへの寄与として、表面エネルギー、クーロンエネルギー、Prolate 形
優勢エネルギーの 3 種類を考慮する。表面エネルギーと sharp-cut な表面の核を使ったクーロンエネルギーを変形パラメータのべき級数で表すと
(

∆Es
∆EC

)

20
2 2
= Es
α2 + α42 + α62 + · · · ,
5
13
(
)
1 2
5 2
25 2
= EC − α2 − α4 −α + ···
527
169 6(3.38)(3.39)となる。ここで、Es と EC はそれぞれ考慮中の核種の球形状態の表面エネルギーとクー
ロンエネルギーである。これら Es 、EC には、Weizsäcker-Bethe 質量公式で使われるものを使用する。

Es = asurf A2/3 ,2EC = a C Z A−1/3(3.40)(3.41)

ここで asurf 、asI はパラメータである。
平均変形エネルギーへの 3 番目の寄与は Prolate 形優勢エネルギーの項である。このような項は他の模型には見られないが、プロレート 変形の優勢という実験的事実は、この項無しには再現できない。小浦らはこの項に下記の形を仮定した。

∆Eprl = −Cprl1 α2 A2/3 exp[−Cprl2 α22 ](3.42)

ここで、Cprl1 と Cprl2 は調節可能なパラメータである。
以下で (3.42) 式の形を検討し よう。まず第一に、この項は、α2 < 0 で正の値を、
α2 > 0 で負の値をとるべきであるから、α2 のべき展開の初項は α2 の奇数のべき乗に負の係数を乗じたものである。小浦らがこのべきを 1 乗にとった理由は次のようなものである。
1 個の核が同じ |α2 | でオブレート 変形した場合とプロレート 変形した場合の形状を比較すると、Rmin はオブレート 変形のほうがプロレート 変形より短い。
（なぜならば、オブレート 形は２つの軸方向に伸び、１つの軸方向に縮むので、体積保存により、伸び
より縮みが強くなる。プロレート 形は２つの軸方向に縮み、１つの軸方向に伸びるので、縮みより伸びが強くなる。）これがプロレート 変形優勢の理由であると仮定すると（なぜ小浦らがそう仮定するかと言えば、巨視的な理由で簡単に否定できないものが他に見当たらないからだろう）、縮みは α2 の 1 乗に比例するのでこの項も α2 の 1 乗に比例することになる。
次の因子 A2/3 は、核変形が表面効果であることを反映している。
最後因子のガウス関数は、|α2 | の値が大きい場合にこの項の増加率を緩くするために導入された。小浦らの考えでは、|α2 | が大きくなると Rmin が小さくなりすぎて、その変化はもはや全エネルギーにあまり影響を及ぼさなくなるだろうからである。
37

平均変形エネルギーは、以上の 3 つのエネルギーの和として与えられる。

Ēdef = ∆Es + ∆EC + ∆Eprl .

(3.43)

様々なパラメータの値はこのように決定される。[19]h = 0.46

as = 17.3MeV

Cprl1 = 0.28
3.2.3

asI = 5MeV,

Cprl2 = 5.

(3.44)(3.45)

変形を含めた原子核の質量の計算

混ぜ合わせの重み W (N1 ) を用いて固有殻エネルギー Ein (N, Z) を求めるには以下のようにする。球形平均場エネルギーを Emf (N, Z) と表し、球形平均場解のエネルギーの平均的 (大域的) 挙動を Eg (N, Z) と表す。それらを使って変形の固有殻エネルギーは、

Ein (N, Z) =∑W (N1 ){(Emf (N1 , Z1 ) − Eg (N1 , Z1 )) − (Emf (N, Z) − Eg (N, Z))}

N1

(3.46)と表せるとする。但し 、Z1 = ZN1 /N である。球形の場合に Ein (N, Z) = 0 となるように配慮した定義である。
Eg (N, Z) には、核図表全面分の球形平均場解のエネルギーの計算値 {Emf (N, Z)} に
単一の Weizsäcker-Bethe 公式を最小二乗法でフィットしたものを用いる。この論文で示した計算では、フィッティング対象に含めた核種は、8 ≤ N ≤ 184, 8 ≤ Z ≤ 126のうち束縛した球形解のある全ての偶々核である。
なお、小浦らの模型では、殻エネルギーは単一粒子準位の和であるため、その振る舞いは現実の原子核の束縛エネルギーとは大きく異なっており、Weizsäcker-Bethe
公式をフィット することは全く不可能である。このために、小浦らは Strutinsky 法によって平均的挙動を求めているが、そのための代償として、フェルミ準位より高い一粒子準位が多数必要になり、連続状態の扱いに苦慮していると思われる。
原子核の基底状態の殻エネルギー Esh (N, Z) は、固有殻エネルギー Ein (N , Z ) と平均
変形エネルギー Ēdef (N, Z) の和の最小値として与えられる。
[]

Ein (N, Z) + Ēdef (N, Z)
Esh (N, Z) = α min
,α ,α
2

4

6

(3.47)となる。これによって、核の全エネルギー Etot (N, Z) は球対称平均場模型および KUTY
流の近似的な変形の扱いに基づく新しい原子核質量公式によってEtot (N, Z) = Emf (N, Z) + Esh (N, Z) + Eeo (N, Z)(3.48)となる。ここで、Eeo (N, Z) は偶奇エネルギーのことでEeo (N, Z) = Ceo1 {1 + Ceo2 exp[Ceo3 (N − 2Z + Ceo4 )/A]}−1/2 A−1/2 δoddZ+ Ceo1 {1 + Ceo2 exp[Ceo3 (N − 2Z + Ceo4 )/A]}−1/2 A−1/2 δoddN− ∆Moo (A)δoddZ δoddN ,38(3.49)

∆Moo (A) = Ccc1 A−1 + Ccc1 A−4/3 + Ccc1 A−5/3 ,δoddZ(3.50)

 1 Z 奇数= 
0 Z 偶数, 1 N 奇数

(3.51)

δoddN = 
0 N 偶数(3.52)となる。
この論文では、α2 のみに値を与え、α4 と α6 は 0 とする。そして、α2 を変化させてた
ときの固有殻エネルギーがより滑らかな曲線となるように、∆N1 = 1 を等間隔に 100等分する。その一点一点での球形平均場エネルギーは線形補間によって導きだし 、その実数 N1 に対応した重み W (N1 ) も求め、∆N1 = 0.01 での数値積分として計算を行う。図 3.8 は N1 が刻み幅 1 の場合の比較と刻み幅 0.01 の場合である。刻み幅を短くすれば 、線がより滑らかになることがわかる。
得られた変形 α2 に対するエネルギー曲線の形は、oblate 側の浅い極小、球形での極大、prolate 側の深い極小 (最小) を備えている点は、もっともらしい結果だと言える。
ただし 、最小点における α2 が 0.2 程度であることは、0.3 程度という通常信じられてい値より有意に小さく、今後の検討が必要と思われる。

-1332

Energy Etot (MeV)

-1334
-1336
-1338
-1340
-1342
-1344
-1346
-1348
-0.4

-0.3

-0.2

-0.1

0

α2

0.1

0.2

0.3

図 3.8: エルビウム (N = 98, Z = 68) で α2 を変化させたとき、実数 N1 の核の全エネルギー

39

0.4

第4章

平均場模型

原子核は有限個の核子からなる量子力学的多体系である。平均場模型とは、多体の相互作用を直接解くことが非常に困難なため、1 個の核子が他の全ての核子によって作られる平均ポテンシャル (平均場) の中を独立に運動しているという猫象が近似的に成立しているとする模型である。

4.1

Skyrme-Hartree-Fock 法

4.2

Hartree-Fock 法

まずフェルミオン多体系を考える。この系はハミルト 二アン H で記述されるとする。よく知られているように、Schrödinger 方程式

H|Ψi = E|Ψi(4.1)

は、変分方程式 (variational equation)

δI[Ψ] = 0,I[Ψ] =

hΨ|H|Ψi
hΨ|Ψi

(4.2)

と同等である。すなわち、|Ψi に何らの制約をおかず、任意の変分 |Ψi → |Ψ + δΨi に対して汎関数 I[Ψ] が停留値をとるならば、そのときの |Ψi は等式 (4.1) の解であり、そのときの汎関数の値がエネルギー固有値である。
等式 (4.2) における |Ψi として試行関数 (trial function)|Φi を採用し、系の基底状態に対する近似解を求める。試行関数 |Φi として単一の Slater 行列 (Slater determinant)をとる近似法が、以下で説明する Hartree-Fock 法である。
いま 1 粒子表示を決めるため、適当な 1 粒子ポテンシャル U (q) を導入する。変数 qは空間座標、スピン座標などをまとめたものである。1 粒子 Schrödinger 方程式
}

{

h̄ 2
∇ + U (q) ϕk (q) = εk ϕk (q)
−
2M

(4.3)

の解 ϕk (q) は規格直交系を作るものとし 、ĉ†k を状態 k に粒子を生成する生成演算子(creation operator)、ĉk を消滅演算子 (annihilation operator) とする。これらはフェルミオンの反交換関係 (anti-commutation relations){ĉk , ĉ†k0 } = δk,k0 ,{ĉk , âk0 } = {ĉ†k , ĉ†k0 } = 0(4.4)

を満たすことはいうまでもない。この 1 粒子表示を使って系のハミルト 二アンは

H=

∑
k1 k2

tk1 k2 ĉ†k1 ĉk2 +1 ∑vk1 k2 k3 k4 ĉ†k1 ĉ†k2 ĉk3 ĉk4
2 k1 k2 k3 k4
40

(4.5)

と表される。右辺の第 1 項が運動エネルギーで、第 2 項が 2 粒子間の有効相互作用である。vk1 k2 k3 k4 は反対称化された有効相互作用の行列要素でΨ∗k1 (q1 )Ψ∗k2 (q2 )v(q1 , q2 )Ψk3 (q1 )Ψk4 (q2 )dq1 dq2 ,

(4.6)

1
{hk1 k2 |v|k3 k4 i − hk2 k1 |v|k3 k4 i − hk1 k2 |v|k4 k3 i + hk2 k1 |v|k4 k3 }
4

(4.7)

hk1 k2 |v|k3 k4 i =v k1 k2 k3 k4 =∫ ∫と与えられる。
Hartree-Fock 基底状態 |Φ0 i は等式 (4.3) の解のうち、エネルギー準位 εk に低い方から順番に A 個の粒子を詰めた Slater 行列式で表され、|Φ0 i =

A∏ĉ†i |0i

(4.8)

i=1

と書ける。ここで |0i は粒子がまったくない状態、すなわち真空である。Hartree-Fock
基底状態 |Φ0 i における最も高い 1 粒子エネルギーを Fermi エネルギー (Fermi energy)と呼び εF で表す。
Hartree-Fock 基底状態 |Φ0 i に対する粒子 (particle) および空孔 (hole) を定義する。εk > εF の 1 粒子状態を粒子状態 (particle state)、εk ≤ εF を空孔状態 (holestate) と呼ぶ。関数

θα =


 0,

ε k > εF(4.9)

 1, εk ≤ εFを使って、生成演算子 ĉ†k をĉ†k = (1 − θk )ĉ†k + θk ĉ†k = â†k + b̂k(4.10)と書くと、右辺の演算子 â†k は粒子状態 α に粒子を生成する演算子であり、演算子 bα は空孔状態 α の空孔を消す演算子となる。

âk |Φ0 i = 0 b̂k |Φ0 i = 0(4.11)

であるから、Hartree-Fock 基底状態 |Φ0 i は粒子および空孔に対する真空となる。 この粒子・空孔表示を用いると等式 (4.5) のハミルト ニアンはH = H0 +∑εk1 k2 (â†k1 âk2 − b̂†k1 b̂k2 ) +

k1 k2

∑

εk1 k2 (â†k1 b̂†k2 − b̂k2 âk1 )

k1 k2

1 ∑
+
vk k k k : ĉ† ĉ† ĉk ĉk :,2 k1 k2 k3 k4 1 2 3 4 k1 k2 3 4∑
1∑
H0 = hP hi0 |H|P hi0 i =tk1 k1 θk1 +vk k k k θk θk ,2 k1 k2 1 2 1 2 1 2k1
εk1 k2 = tk1 k2 + 2∑

vk1 k3 k2 k3 θk3(4.12)(4.13)(4.14)

k3と書き直すことができる。ここで等式 (4.12) における : ĉ†k1 ĉ†k2 ĉk3 ĉk4 : のように、演算子を : の記号ではさんだ表式は、粒子・空孔に関する正規積 (normal product) であ41

る。正規積とは、生成演算子を最左方に、消滅演算子を最右方に並べる演算子の積であり、Wick の定理を用いて、演算子の積を正規積の方にすると、定数項が真空期待値となる。
Hartree-Fock 法は上記の Slater 行列式 |Φ0 i を構成する 1 粒子波動関数 |φk i = |kiを変分関数として、系のハミルト ニアン H の期待値 H0 = hΦ0 |H|Φ0 i を極小にするように |φk i を決定する近似である。規格条件 hφk |φk i = hk|ki = 1 を付すための Lagrange
未定乗数として εk をとれば 、変分方程式はδ hk|t|ki +∑(hk1 k2 |v|k1 k2 i − hk1 k2 |v|k2 k1 i) − εk hk|ki

(4.15)

k2と書かれる。その結果、1 粒子状態を決める Hartree-Fock 方程式 (Hartree-Fockequation)∫h̄2 2−∇ φk (q) + U (q, q 0 )φk (q 0 )dq 0 = εk (q)(4.16)
2Mを得る。ただし 、変分方程式 (4.15) の主旨からわかるように、この方程式は空孔状態に対して導かれたものであるが、粒子状態まで拡張するものとする。すなわち、1 粒子状態 k は空孔・粒子のすべての状態を意味するものとする。また、非局所的な平均ポテ
ンシャル、すなわち Hartree-Fock ポテンシャル (Hartree-Fock potential)U (q, q 0 )は
00

U (q, q ) = δ(q − q )
−

∑
k0

∑

∫

θk 0

φ∗k0 (q 00 )v(q 00 , q)φk0 (q 00 )dq 00

k0
θk0 φ∗k0 (q 0 )v(q 0 , q)φk0 (q)(4.17)

で与えられる。Hartree-Fock 方程式 (4.16) は非線形方程式であり、平均ポテンシャ
ル U (q, q 0 ) の中に解くべき 1 粒子波動関数 φk が入っているので、自己無撞着的 (selfconsistent) に解かなければならない。
1 粒子表示をきめるために使った 1 粒子 Schrödinger 方程式 (4.3) の解 φk の代わりにHartree-Fock 方程式 (4.16) の解をとけば等式 (4.14) の式の εk1 k2 は、k1 または k2 の一方が粒子状態で他方が空孔状態のときには 0 となるから、等式 (4.12) の右辺の第 3項となる。したがって、Hartree-Fock 近似のもとで、系のハミルト ニアンはH = H0 + H1 + H2 ,H0 = hΦ0 |H|Φ0 i =
H1 =

∑
k1

H2 =

∑
k1

εk1 (â†k1 âk1 −tk1 θk1 −∑∑(4.18)
vk1 k2 k1 k2 θk1 θk2 ,(4.19)k1 k2

εk1 (b̂†k1 b̂k1 ,(4.20)k1

1 ∑vk k k k : ĉ† ĉ† ĉk ĉk :,2 k1 k2 k3 k4 1 2 3 4 k1 k2 3 4(4.21)と書かれる。H0 が Hartree-Fock 基底状態のエネルギー、H1 が 1 粒子 (空孔) ハミル
ト ニアンであり、全ハミルト ニアンのうち Hartree-Fock 近似のもので H0 や H1 に入りきれなかった残留相互作用が H2 である。[22]424.2.1

Skyrme 相互作用

原子核は多数の核子が核力と (nuclear force) いう強い短距離相互作用によって結合した多体系で、その核力は平均ポテンシャルを作り出している。しかしながら、核力の効果がすべて完全に平均ポテンシャルになるわけでわなく、平均ポテンシャルに吸収されない相互作用も少しは考えられる。これを残留相互作用 (residual interaction)と呼ばれる。基底状態や励起状態は、閉核でない不完全な殻 (オープン殻) に与えられた個数の核子が配置されるすべての可能な状態 (配位) の間に働く残留相互作用を考慮して決まる。その相互作用エネルギーが最も低くなる状態が基底状態である。
残留相互作用は通常有効相互作用 (eﬀective interaction) と呼ばれ、その現象論的モデルとして Skyrme 力と呼ばれるものがある。Skyrme 力は 2 体力と 3 体力からなり、その Skyrme 相互作用を含めたハミルト 二アン Ĥ は

Ĥ =

∑

t̂i +

∑

i

(2)

vij +

i<j

∑

(3)

vijk .

(4.22)

i<j<k

2 体力の項は、(2)

vij= t0 (1 + x0 P̂σ )δ(~r1 − ~r2 )
(
)
1
0
~
+ t1 δ(~r1 − ~r2 )k̂ 2 + k̂ 2 δ(~r1 − ~r2 )2ˆˆˆ+ t2~k · δ(~r1 − ~r2 ) + iW0 (σ~ˆ1 + σ~ˆ2 ) · k~0 × δ(~r1 −~r2 )k~0 .

(4.23)

そして、3 体力の項は、

v(3)ijk = t3 δ(~r1 − ~r2 )δ(~r2 − ~r3 )(4.24)

ˆˆパラメータ ~k と k~0 のと定義は~kˆ = 1 (∇1 − ∇2 ),21
ˆ
k~0 = − (∇1 − ∇2 ).
2(4.25)~kˆ は左側に作用する。
Hartree-Fock 方程式を得るために、Slater 行列式 |HFi のハミルト 二アンの期待値を求めなければならない。
E = hHF|H|HFi=∑
i

hi|t̂|ii +1∑
1∑
hij|v̄ (2) |iji +hijk|v̄ (3) |ijki.
2 ij
6 ijk

(4.26)

ハミルト 二アン密度の空間積分して書き直すと∫E=

d3 r Ĥ (~r ).

(4.27)

ハミルト 二アン密度は単一粒子状態の総和によって得れる一様な密度汎関数である。その汎関数をφi (~r, σ, q),43(4.28)

とおけば 、次に述べる密度を定義するのに役立つ。ここ σ = ± 21 はスピン射影を表し 、
q = ± 12 (表記 q = p, n は添字 q の固定された値として用いられる) はアイソスピンを表す、そして i は軌道関数だけを数え上げる。以下に記述するのは φi (~r, σ, q) の定義である。

1. 陽子と中性子についての密度ρq (~r) =∑|φi (~r, σ, q)|2 ,(4.29)

ı,σ2. 陽子と中性子についての運動エネルギー密度
τq (~r) =∑|∇φi (~r, σ, q)|2 ,(4.30)

ı,σ3. スピン軌道のカレント
J~q (~r) = −i∑
ı,σ,σ 0ˆ |σ 0 i],
φ∗i (~r, σ, q)[∇φi (~r, σ 0 , q) × hσ|~σ(4.31)

運動エネルギー密度が、因子 h̄/2m を含まず、通常の組合せ φ∗ ∇2 φ を必要としない。
ここで用いる、より対称的な組合せは全ハミルト 二アンと同等のものである。部分積
分型が示すようなスピンカレント 密度は極めて複雑に思える。それは挿入する単一波動関数は、本質的に ρ~k × ~σ の項、すなわち、スピンのクロス積した局所確率密度による重みとして記述することを示す。
その和は占有された全ての単一粒子状態を引き継がれるべきである。これから、占有状態の空間は時間反転に対して不変であるという仮定を作る。このことは偶偶核の静的性質において最も重要な適用と言える。形式上、仮定が意味することは、もし 状態
φi (~r, σ, q) が占有されるならば 、時間反転演算子 T̂ を使って得られる時間反転状態も、また占有されることにである。スピノール項における T̂ = −iσ は、軌道波動関数における複雑な連結を結合させる結果を使うと、

T̂ φi (~r, σ, q) = −i∑σ0hσ|σy |σ 0 iφ∗i (~r, σ 0 , q) = −2σφ∗i (~r, −σ, q).

(4.32)

その結論はスピン消滅演算子を含んだ期待値である。
∑
i,σ,σ 0ˆ |σ 0 iφi (~r, −σ 0 , q) = 0.
φ∗i (~r, σ, q)hσ|~σ(4.33)

このことを理解するために時間反転対称の式に演算子 21 (1 + T̂ ) を使う。
∑
iφ∗i (~r, σ, q)φi (~r, σ 0 , q)1∑ ∗(φi (~r, σ, q)φi (~r, σ 0 , q) + 4σσ 0 φ∗i (~r, −σ, q)φi (~r, −σ 0 , q))
2 i
)
∑(
1
|φi (~r, σ, q)|2 + |φi (~r, −σ, q)|2
= δσσ02i
1= δσσ0 ρq (~r)
2=44

(4.34)(4.35)

つまり、要望される期待値の求めることができる。
∑

ˆ |σ 0 iφi (~r, σ 0 , q)φ∗i (~r, σ, q)hσ|~σ

i

=

1∑
ˆ |σ 0 i
δσσ0 ρq (~r)hσ|~σ
2 σ,σ01ˆ} = 0= ρq (~r)Tr{~σ(4.36)
2
前式は Pauli 行列の消滅跡から起こる。この計算はその跡を明確にする。もし 2 つの占有された時間反転状態を持つ状況が満たさないとき、形式は余りにも複雑になるだろう。なぜなら、単一スカラー密度はスピンの添字を使う密度行列と取り換えらればならない。
運動エネルギーの項を含むエネルギー汎関数の展開の結果は
[(

)

(

)

]

)
h̄2
1
1
1 ( 2
Ĥ(~r) =τ + t0 1 + x0 ρ2 − x0 +
ρn + ρ2p
2m
2
2
2
1
1
+ (t1 + t2 )ρτ + (t2 − t1 )(ρn τn + ρp τp )4
8
1
1
+ (t2 − 3t1 )ρ∇2 ρ + (3t1 + t2 )(ρn ∇2 ρn + ρp ∇2 ρp )
16
32
1
1
+ (t1 − t2 )(J~n2 + J~p2 ) + t3 ρn ρn ρ
16
4
1
− W0 (ρ∇ · J~ + ρn ∇ · J~n + ρp ∇ · J~p ).
(4.37)
2
クーロン効果は、上式に含まない。それは、クーロン交換に対応する近似を使うクー
ロンエネルギーは単独に付け加えねばならない。
とても簡単な形は N = Z の核種について有効である。クーロン効果の削除、ρn =ρp = ρ と仮定できる。そしてもう一方の密度も同様にすると

h̄2
3
1
1
Ĥ(~r) =τ + t0 ρ2 + t3 ρ3 + (3t1 + 5t2 )ρτ
2m
8
16
16
1
3
1
2
+ (9t1 − 5t2 )(∇τ ) − W0 ∇τ · J~ + (t1 − t2 )J~2 ,(4.38)64
432
両方の場合、エネルギー式は単一粒子波動関数に関して変分により得られる HartreeFock 方程式に向かっての出発点であり、これは標準 Hartree-Fock の場合では積分微分型方程式を出す。

4.3

相対論的平均場模型

相対論的平均場模型の構成要素としての中間子場の意味を論じた後、模型の基礎にあたる 2 つの近似について説明する。

4.3.1

核子の媒体としての有効中間子場

相対論的平均場模型の目的は、原子核系の相対論的な記述である。原子核の内の核子間に働く直接的瞬間的な力という概念は、相対論的な定式化には適切ではない。核子
45

間相互作用は、核子とは独立した自由度である中間子場によって媒介されると考えるほうが相対論との親和性がよい。低エネルギー領域において、中間子場はそれらのスピン(内部角運動量)J 、偶奇性 (空間反転操作により生じる状態の符号変化)Π とアイソスピ
ン (スピンに類似した荷電空間内の自由度)T による指定が可能であり、小さい J と T の値をもつ中間子場が核子間相互作用を支配していることになる。そのことは、中間子の質量スペクト ルの低エネルギー見ればわかる。但し、高エネルギー領域までを視野に入れると同じ J, Π, T をもつ中間子は複数種あり、大きな J と T の値をもつ中間子場も存在する。また、核子間に働くポテンシャルエネルギーが one − boson − exchange ptentialの和として高精度に再現できることからもわかる。そこで、原子核の相対論的平均場
模型にあいて通常では、スカラー中間子 (J = 0) とベクト ル中間子 (J = 1)、アイソ
スカラー (T = 0) とアイソベクト ル (T = 1) 中間子のみを考慮する。さらに、良い偶奇性を持つ原子核状態を研究の対象とすることが多く、その場合「不自然な偶奇性」(Π = −(−1)J ) を持つカレント の期待値は消えるので、「自然な偶奇性」(Π = (−1)J )を持つ場のみを取り入れる。これは、pi 中間子 (J = 0, Π = −, T = 1) や eta 中間子(J = 0, Π = −, T = 0) の場を除外することを意味する。具体的には、

σ
ω
ρ
γ

:
:
:
:

アイソスカラー・スカラー場 Φ(xµ )アイソスカラー・ベクト ル場 Vν (xµ )
アイソベクト ル・ベクト ル場 Rν (xµ )
質量・電荷を持たないベクト ル場 Aν (xµ ), 光子を考える。 σ,ω,ρ と言った名前は、各スピン・アイソスピンを持った中間子の代表的なものの名前からとった取ったものであるが、同じスピン・アイソスピンを持つ他の種類の中間子の役割をもカバーする必要がある為に、純粋な σ,ω,ρ 中間子とは質量や核子との結合定数が異なったものとなる。ρ 場は原子核のアイソベクト ル依存性を定量的に再現する為には導入しないでは済まされず、光子ももちろん既知の電磁気相互作
用を説明する為に必要である。アイソベクト ル・スカラー場 (δ 中間子) は、模型を改善しないので選択から省かれる。

4.3.2

中間子の古典場近似

次に近似の骨格を、最も単純な例すなわち、丁度１つのスカラー場 Φ だけによって束縛された核子 (その場を ψ する) から成り立っている模型の場合について説明する。
その模型は、ラグランジアン密度(free)

L = ψ̄(γ µ ∂µ − mB )ψ + L(meson) − gσ Φψ̄ψ

(4.39)と記述される。場の理論の研究においても、このラグランジアン密度で表される模型を使うことがあるが、その場合は Φ や ψ の量子場としての振る舞いが興味の対象となる。原子核の相対論的平均場模型における意図は、それと異なっており、2 つの基礎となる近似、中間子場の平均場近似と核子のフェルミの海を無視する近似を、前提として初めて意味を持つものとしてこの模型のラグランジアンとして扱われる。

46

平均場近似は、中間子場の全ての量子的揺動を無視し、期待値だけに注目するという近似である。つまり、全ての中間子場を古典的 c− 数 (classical number) の場として扱う近似である。一般に、中間子場の扱いに関しては、平均場近似という語句と古典
場近似という語句は全く同じ意味で用いられている。この近似を形式的に表現すれば、

Φ̂ → Φ = hΦ̂i(4.40)と言うことである。その結果、核子は中間子の平均場を経由してのみ相互作用することができることになる。そのメカニズムを述べると、まず中間子場が核子を束縛する(中間子場に結合定数を乗じたものが核子の感じるポテンシャルエネルギーとなる)。核子はそれぞれ単独に中間子場に反作用を返す。その影響を核子の個数分蓄積したものが中間子場のソース項 (非斉次項) となる。このようにして核子と中間子場は互いに自己無撞着に決定されることになる。

4.3.3

負エネルギー状態の核子の無視中間子場の平均場近似は、核子の処理をも容易にする。すなわち、核子は中間子場で独立した粒子として運動することになる。核子場の演算子 ψ を

ψ=∑ϕα (xµ )âα(4.41)

αとして、単一粒子状態 α の完全系で展開して表しておく。ここで、âα は状態 α の核子に対する消滅演算子であり、ϕα (xµ ) は単粒子状態の波動関数である。そうすると、ス
カラー密度の自由核子の真空からのずれは、

h: ψ̄ψ :i =∑ϕ̄α ϕα −

α<F

∑

free
ϕ̄freeα ϕα(4.42)

α<F0と表される。ここで、ϕfree
α は Dirac 方程式の自由 (平面波) 解であり、F0 は核子数に対応するフェルミ準位で、F は与えられた原子核の核子数 A に対応するフェルミ準位である。等式 (4.42) の和は、反粒子状態の Dirac の海にある全ての核子を含んでいる。
この無限和は非常に扱いにくい。それは、全ての物理量に対して真空分極効果の評価を要求する。しかしながら、核子と中間子はそれら自身が複合粒子であり、基礎的な場の理論の構成要素ではない。したがって、それらの量子場効果の一部分だけを取り入れることの意味は明確ではないので、有効場であるという理由で量子的揺動を無視するのが、筋の通った方法論である。この目的に至る為に等式 (4.42) の総和は

h: ψ̄ψ :i = ∑α<Fϕ̄α ϕα −∑α<F0
free 
+ϕ̄freeα ϕαA∑ϕ̄α ϕα(4.43)

α=1と書き直す。ここで、F0 と F の間にある核子の状態数が A に等しくなるように F は決めなければならない。上式の第 1 項は真空分極効果を分離し 、第 2 項は原子核内の核子の分布を記述する。第 1 項の真空分極効果を無視することは、第 2 項の α = 1, · · · , A

47に対応するわずかな A 個の占有された粒子状態だけに総和を制限することを意味する。
これは no-sea 近似と呼ばれるものである。それを形式的に表せば 、

h: ψ̄ψ :i 'A∑ϕ̄α ϕαα=1となる。スカラー密度以外のタイプの密度にもこの近似が同様に適用される。

48(4.44)

第5章

新しい質量公式

前章で完成させた新しい原子核の質量公式プログラムに、球形平均場解の質量計算
値を入力として読み込ませて、変形を考慮した原子核の質量を計算して値を出力する。
そして、G.Audi と A.H.Wapstra のまとめた実験値との平均二乗誤差を出して、質量公式の評価を行う。
Audi-H.Wapstra の表にある質量のうち予測データは使用せずに精度の認められた
実験値のみを使う。また、中性子数 N と陽子数 Z は共に 8 以上の核種の質量を使う。
以上の制限により、使用可能な核種の質量は 2147 個となる。
球形平均場模型としては、さしあたり Skyrme-Hartree-Fock 法を用いる。Skyrme
相互作用には幾つかのパラメータの組が 存在する。それらの現在までに提唱された
Skyrme 力のパラメータセット から SIII を選ぶそれらの理由は、変形平均場解と実験値との質量の平均２乗誤差がすでに計算されていることである。
2147 個の核種の質量値から Weizsäcker-Bethe 質量公式の各項のパラメータを求める。その質量公式には偶奇項を含めないものとする (偶偶核データだけを使うならば含
めない、そうでなければ含めるべきである)。それは、Skyrme-Hartree-Fock 法で球形核の質量を導き出すに当たり、偶偶核、奇核、奇奇核それぞれを分けて考慮しなければならない。しかし 、本論文では奇核も奇奇核も偶偶核とみなした値だからである。
G.Audi と A.H.Wapstra の実験値と Weizsäcker-Bethe の質量公式の値の差を表したものが図 5.1 である。これから分かることは、と Weizsäcker-Bethe の質量公式は実験値と比べて、陽子数と中性子数が魔法数を中心に誤差が大きくなることがわかる。それは第 2 章にも述べたが、Weizsäcker-Bethe の質量公式が殻効果を考慮していないか
らである。
次に、Skyrme 相互作用の球形核解の質量値と Weizsäcker-Bethe の質量公式の質量の差を表したものが図 5.2 である。
図 5.1、5.2 は同一のパラメータの Weizsäcker-Bethe 公式を基準として、それぞれ実験値と球形平均場解の値を示したものである。Weizsäcker-Bethe 公式からの差とは、殻エネルギーと考えてよく、この図から現実の原子核の質量の特徴が、また理論模型の特徴が読みとれる。
図 5.1 からは魔法数を中心に安定性が大きいことが明瞭に分かる。図 5.2 でも、魔法数は同様に再現されているが、図 5.1 との違いが目立つのはむしろ、魔法数と魔法数の中間の領域である。この領域において図 5.2 の質量値は異常に大きくなっている。この領域の核では、現実の原子核では変形によってエネルギーは低下させるのに対し 、球形平均場法ではそれができないためである。
3 番目に、本論文で作成した質量公式と Weizsäcker-Bethe の質量公式の値の差を表
49

AudiWapstra-BW (MeV)
120

20
15
10
5
0
-5
-10
-15
-20

100

Z

80
60

40

20
0
0

20

40

60

80

100

N

120

140

160

180

図 5.1: G.Audi と A.H.Wapstra の実験値と Weizsäcker-Bethe の質量公式の値の差

SphericalS3-BW (MeV)
120

20
15
10
5
0
-5
-10
-15
-20

100

Z

80
60

40

20
0
0

20

40

60

80

100

120

140

160

180

N
図 5.2: skyrme S3 力の質量値と Weisäcker-Bethe の質量公式の値の差50したものが図 5.3 である。図 5.1 は図 5.2 と違い、魔法数と魔法数の中間の領域で質量が異常に大きくなってはおらず、図 5.2 よりは実験値の図 5.1 に似た分布になっている。これは、新しい質量公式が変形を考慮した成果である。
図 5.1〜図 5.3 は、Weizsäcker-Bethe 公式を基準とすることで、実験値および各種理論模型の特徴を表現したものであったが、図 5.4 と図 5.5 では、実験値を基準として、各種理論模型のそこからの差 (即ち、誤差あるいは精度) を表したものである。

S3+Koura-BW (MeV)
120

20
15
10
5
0
-5
-10
-15
-20

100

Z

80
60

40

20
0
0

20

40

60

80

100

N

120

140

160

180

図 5.3: KUTY 流質量公式の値と Weizsäcker-Bethe の質量公式の値の差

図 5.4 は、Skyrme 相互作用の球形核平均場解の質量値と G.Audi と A.H.Wapstraの実験値との核種の質量誤差を表したものである。平均二乗誤差は 5.064MeV である。魔法数付近の誤差は小さくなっているが、それ以外のところでの誤差が大きいことがわかる。すでに述べたように、これは変形核の領域である。
図 5.5 は、本論文の新公式の値と G.Audi と A.H.Wapstra の実験値を同様に比較したものである。平均二乗誤差は 1.607(MeV) となり、先程の平均二乗誤差の約 30％まで低減され、それだけ実験値に近づいたといえる。KUTY 流変形の扱い方はいく
つかのパラメータが存在するが、本研究では (3.44) 式のパラメータ aS の値の１個のみパラメータとして選び、その値を小浦らの使用した値の −34.0% にした場合に、平
均二乗誤差が一番小さくなった。図 5.4 とは異なり、魔法数近傍も、変形核領域も、同程度に誤差が小さくなっていることが見てとれる。

51

SphericalS3-AW (MeV)
120

20
15
10
5
0
-5
-10
-15
-20

100

Z

80
60

40

20
0
0

20

40

60

80

100

N

120

140

160

180

図 5.4: skyrme S3 力の質量値と G.Audi と A.H.Wapstra の実験値の差

S3+Koura-AW (MeV)
120

20
15
10
5
0
-5
-10
-15
-20

100

Z

80
60

40

20
0
0

20

40

60

80

100

N

120

140

160

180

図 5.5: KUTY 流質量公式の値と G.Audi と A.H.Wapstra の実験値の差52

第 6 章 まとめ
本論文では、球対称平均場模型に近似的な変形の扱い方を組み合わせた新しい原子
核の質量公式を提唱し 、その実用性を調査した。試行的なパラメータフィッティングを行った結果から、その可能性の高さが示された。
最近の質量公式は平均場模型に基づくもの [27] が主流であるが、変形核を扱うために計算規模が大きくなり過ぎてパラメータの最適化を十分に行うことが難しい傾向がある。一方、小浦ら [19] は、一昔前に隆盛を誇った Strutinsky 法に基づくタイプの質量公式の新規開発に際して、変形核を球形核の重ね合わせとして近似するという斬新なアイデアを発案した。しかし 、Strutinsky 法は平均場法より計算規模が小さいという利点はあるが、多くの欠点を持っている。
そこで、我々は、新しいタイプの質量公式として、平均場模型に基づくが、球形の平均場解のみを使用し 、変形の影響に限って小浦らの近似的方法を援用するというタイプのものを開発することにした。近似的な変形の扱い方の含むパラメータの個数はせいぜい数個であり、パラメータ数の激増にはつながらない。一方、計算規模は、変形解でなく球形解を求めればよいことから何桁も小さくなる。それに反比例した十分な手間をパラメータの最適化にかけることができるようになる。
本研究では、研究の第一段階として既存の平均場模型をそのまま利用し、変形の近似的取扱い方の部分も小浦らの方法を１個のパラメータ（ 表面エネルギーの係数）の値
を除いてほぼそのまま踏襲した。表面エネルギーの係数は質量の平均二乗誤差を最小化するように決定した。
実験データには G.Audi と A.H.Wapstra による実験値の表 [13] からとった 2147
個の核質量を使用した。この実験値と Skyrme SIII 相互作用 [28] を用いた HartreeFock+BCS 解の質量値との平均二乗誤差は 5.06(MeV)MeV であった。文献によるとこの相互作用で変形解を求めた結果の平均二乗誤差は 2.2 MeV であった [29]。本論文で開発した新しい質量公式の平均二乗誤差は 1.61MeV となり、表面エネルギーの係数という１パラメータの自由度を利用して、変形を正確に計算した場合よりも実験値に近い質量値を得られることが分かった。ちなみに、小浦らの KUTY 公式の場
合、2921 個の核種に対しての実験値との平均二乗誤差は 0.715 MeV であり、これを越える精度を達成することが次の目標となる。
本研究により、球形核の殻エネルギーの重ね合わせから変形核の質量を導き出すという、極めて現象論的な方法が、その提唱者らの模型だけで有効なのでなく、平均場との組み合わせにおいても十分に有効であることが示された。
今後の課題は、まず、今回は偶々核のみを扱ったが、現象論的に対エネルギーの項を導入することで、奇核・奇々核をも扱える質量公式とすることである。さらには BCS理論などにより微視的平均場模型の部分で対エネルギーを求めることも試みる価値が53

ある。
さらに、平均場模型の部分には、非相対論的模型よりは、核図表上の実験的未知領域への外挿に際し 信頼性がより高いと期待される相対論的平均場模型を使用することが望ましい。
このような改良の後に、平均場模型のパラメータを所与のものとし てでなくフィッ
ティングパラメータとして扱うことで、精度の高い質量公式を創出することへの挑戦が控えている。
近いうちに本格的な新質量公式を世に送り出せることを期待して本論文を終えたい。

54

関連図書
[1] 有馬 朗人:「原子と原子核−量子力学の世界−」 朝倉書店 (1991)pp.162-171[2] JAGDISH K.TULI,NATIONAL NUCLEAR DATE CENTER,NUCLEAR WALLET CARDS(Seventh edition),Brook Nation LaboratoryUpton,New York 119735000,USA (2005)[3] 野上茂吉郎:「原子核」 裳華房 (1973)pp.94-102[4] 八木 浩輔:「原子核物理学」 朝倉書店 (1971)[5] B.Povh,K.Rith,C.Scholz,F.Zetsche:「素粒子・原子核 物理入門」 Springer (1997)[6] 白井 行秀 福井大学 工学部 応用物理学科 卒業論文 「原子核質量の系統性 (液滴模型による解析) 」 1999 年 2 月[7] 池田 智也 清澤 聡 福井大学 工学部 応用物理学科 卒業論文 「原子核の質量公式
における対エネルギー項の改良」 2001 年 2 月[8] C.F.von,Weizsäcker,Z.Physik96,461(1935)[9] H.A.Bethe,F.Bach,Rev.Mod.Phys.8,82(1936)[10] 小川 直久,http://www.lowtem.hokudai.ac.jp/ptdice/ogawa/data/MS.doc[11] J.Rainwater.Phys.Rev.79,432(1950)[12] 金原 寿郎:「電磁気学 (1) 」 裳華房 (1972)pp.125-127[13] G.Audi, A.H.Wapstra and C.Thibault:ATOMIC MASS ADJUSTMENT (2003)
file : mass.mas03, Nucl.Phys.A729,129(2003)[14] 山田 勝美:「原子核はなぜ壊れるか」 丸善 (1987)[15] 山田 勝美、宇野 正宏、橘 孝博:「 原子質量公式−未知の原子核を理論面から探る」
日本原子力学会誌 Vol.42,No.4,pp.245-256 (2000)[16] 宇野 正宏、橘 孝博:「原子核質量の評価データと予測」 日本原子力研究所核デー
タニュース (1988)[17] Masami Yamada, Takahiro Tachibana,Partly Frozen Liquid-Drop of Atomic Nuclei,MEMOIRS OF THE SCHOOL OF SCIENCE & ENGINEERING WASEDA
UNIV.NO46(1982)55

[18] Takahiro Tachibana,Masahiro Uno,Masami Yamada,So Yamada,EMPIRICALMASS FORMULA WITH PROTON-NEUTRON INTERACTION,ATOMICDATA AND NUCLEAR DATA TABLES 39,251-258,(1988)[19] Hiroyuki Koura,Masahiro Uno,Takahiro Tachiro,Masami Yamada,Nuclear mass
formula with shell energies calculated by a new method,Nuclear Physics A674
(2000) 47-76[20] 高橋 大輔:「数値計算」 岩波書店 (1996)[21] 奥村 晴彦: 「 C 言語による最新アルゴ リズム事典」 技術評論社 (1991)p.354-356[22] 高田 健次郎 池田 清美 : 「原子核構造論」 朝倉書店 (2002)p.164-168[23] A・メシア:小出昭一郎, 田村二朗 訳 :「 メシア 量子力学 3 」 東京図書株式会社(1976)p.180-187,217-220[24] 三和 之浩 福井大学大学院 物理工学専攻 修士論文 「原子核の相対論的平均場模型」2008 年 2 月[25] B.D.Serot and J.D.Walecka, THE RELATIVISTIC NUCLEAR MANY-BODY
PROBLEM,Advances in Nuclear Physics,vol16,Plenum Press,New York(1986)[26] C.J.Horowitz and B.D.Serot, SELF-CONSISTENT HERTREE DESCRIPTION
OF FINITE NUCLEI IN A RELATIVISTIC QUANTUM FIELD THEORY,Nuclear Physics A368(1981)503-528[27] S. Goriely and J.M. Pearson, Phys. Rev. C, 77, 031301(R) (2008).
[28] M. Beiner, H. Flocard, Nguyen van Giai, and P. Quentin,Nucl. Phys. A238, 29(1975).
[29] N. Tajima, S. Takahara, and N. Onishi, Nucl. Phys. A603, 23 (1996).

56

謝辞
本論文を作成するにあたり，終始変わらぬ懇切丁寧なご指導と助言頂いた 田嶋直
樹先生、鈴木敏男先生、林明久先生には誠に感謝し ，お礼申し 上げます。特に田嶋直樹先生には最初から最後まで本研究に対し 、貴重な御時間を割いてまで根気強く丁寧に御指導頂いたこと、また本論文作成に関しても、多大な御協力を頂きましたことを、心より感謝申し 上げます。また物理工学科の諸先生方に対しても、日頃の厚い御指導に深く感謝しています。最後に２年間共に学び、励まし 合ってきた原子核理論研究室の先輩、大学院生、卒業研究配属生諸氏にも本当に感謝し 、謝辞の言葉とさせて頂きます。

2009 年 2 月
山田昌平

57

付録：

小浦らの論文からの抜粋の邦訳

以下に小浦らの論文（ 文献 [[19]] ）の第３セクション “Nuclear shell energies ofdeformed nuclei” の最初の２サブセクション “3.1 Essentials of the method oftreating deformed nuclei” と “3.2 Deformed nucleus as a superposition of
spherical nuclei” を私が日本語に翻訳したものを掲載する。この部分は、変形の近似的取り扱いの解説である。この近似は本修士論文の最も重要な要素であり、私自身の理解を深めるために、翻訳を作成しつつじっくりと時間をかけて読んだのである。この近似の私自身の観点からの説明は本修士論文の本文中で既に与えたが、私の理解に不同意な方がいた場合には、即座に原論文を併せて読んでいただければ 便利であろうと考え、ここに採録したのである。

6.0.4

変形核を取り扱う大切な方法

基底状態の質量を一般的な傾向からの寄与と、それらからのずれの和で記述する。
球形核の場合、このずれを殻エネルギーと呼ぶことはごく自然である。変形核の場合、事態はより複雑である。もし 原子核の形状が、３つの軸の比率が与えられた楕円体の場合の様に指定されたならば 、Z と N の滑らかな関数は偶奇効果を除けば 、このような核の最も低いエネルギーの一般的傾向を表す為に決められるだろう。滑らかな関数
から個別の変形核のエネルギーのずれは、変形核の構造のためである。そして、これらのずれを変形核の固有殻エネルギーと呼ぶ。
球形核に反して変形核の場合、様々な変形を考慮しなければならない。夫々の核種について、もし 変形の状態が球形の状態より低いエネルギーを与えるならば 、１つの変形だけが基底状態して実現される。この最小エネルギーで、固有殻エネルギーだけでなく滑らかな関数も考慮しなけらばならない。なぜなら滑らかな関数もまた変形に依存するからである。球形核の滑らかな関数を引いた、滑らかな関数は以下で平均変形
エネルギーと呼ぶ。すなわち最小とされるエネルギーは固有殻エネルギーと平均変形
エネルギーの和である。この和を全殻エネルギー、又は単に殻エネルギーと呼ぶ。
幾つかの仮定 (補正) の助けを用いて、変形殻の固有殻エネルギーを得る為の概略方法を与える。
基本的な仮定の１つは、ある種の量を扱う中で変形核が適当な球形核の重ね合わせと扱うことができる。その重ね合わせについて、干渉を伴わないものだけが意味される。球形一粒子ポテンシャルを取り、その中で一粒子状態を考慮しなさい。変形核状
態は原則として、一粒子状態で記述される。そして、一粒子ポテンシャルレベルの占有率は決められることができる。すなわち様々な陽子数と中性子数を持つ適当な球形
核最低エネルギーでの占有率の重ね合わせによって近似させるだろう。仮定として占58

有率を表すことの可能性と置く。これは対相関の場合よりも多少思いきった仮定である。なぜならこの時の用いられた球形核の占有率は、一般的に sharp-cut ではないか
らである。対相関の場合、占有率は階段関数によって表される。明確にすると、もし変形核での ν 番目一粒子状態の占有率 wiν (i = p, n) と表すとき、上記の仮定は、

wpν (Z, N ) =∑Wp (Z 0 ; Z, N )w0pν (Z 0 , N 00 ),(6.1)

Wn (N 0 ; Z, N )w0nν (Z 00 , N 0 ),

(6.2)

Z0

wnν (Z, N ) =∑N0と表される。ここで N 00 と Z 00 は夫々N Z 0 /Z と ZN 0 /N に最も近づく整数である。係数
wp ,wp は規格化条件を満たしている混ぜ合わせた重みである、∑

Wp (Z 0 ; Z, N ) = 1,∑

Wn (N 0 ; Z, N ) = 1,

(6.3)

N0

Z0

次の仮定は、変形核の固有殻エネルギー Ein (Z, N ) が等式 (6.2) での様に、同じ混ぜ合わせの重みで同じ球形核の陽子殻エネルギーと中性子殻エネルギーの重ね合わせとして表されることができる。つまり、

Ein (Z, N ) =∑Wp (Z 0 ; Z, N )Eps (Z 0 , N 00 ) +Z0∑

Wn (N 0 ; Z, N )Ens (Z 00 , N 0 ).

(6.4)

N0

この仮定は部分的に占有率についての最終仮定と関係付けられる。もし 球形な陽子殻エネルギーと中性子殻エネルギーが一粒子レベル εν とこれらのレベルの平均的位置
εν の間の差より引き起こされるならば 、等式 (6.4) は等式 (6.1)、(6.2) からの結果として生じる。非補正殻エネルギーはこの状態を満足させるけれども、対相関の補正は他のエネルギー寄与つまり ∆E npr と ∆E ppr を導入する。これらの寄与は核種に依存して０から数 MeV まで及び負でないのを取る。変形核では対相関から同じ寄与であるだろう。すなわち、等式 (6.4) の仮定では変形核ではこの種の核の寄与が等式 (6.4) での同じ重みを用いて、µ∆E npr ,µ∆E ppr (µ は減衰因子である) の重ね合わせとして表させれると仮定する。
重ね合わせとして表させると仮定させる第３の量は全方向で平均された陽子 (中性子)
数密度である。変形核でのこの量は ρi (r; Z, N )(i = p, n) と表し、適当な球形核 (Z 0 , N 0 )でのものを ρ0i (r; Z 0 , N 0 )(i = p, n) と表す。すなわち,

ρp (r; Z, N ) =∑Wp (Z 0 ; Z, N )ρ0p (r; Z 0 , N 00 ),(6.5)

Wn (N 0 ; Z, N )ρ0n (r; Z 00 , N 0 ),

(6.6)

Z0

ρn (r; Z, N ) =∑N0と仮定する。ここで、混ぜ合わせの重みは等式 6.1、6.2 でのものと同じである。これは方向的な平均密度が粒子レベルの占有率によって殆んど一意に決まるから自然であ59rΩoc(r)/2rimΩim(rim)/2

Rmin

図 6.1: r から Rim までの半径座標の変換の説明図と、中間形状の導入る。この主張は主殻、又は主殻に隣接した２つの一粒子状態の間の干渉が、全方向で平均される時に消滅する事実に基ずく、それはこれら２つの状態が同じ角運動量と同じパリティを持つことができないからである。
等式 (6.5)、(6.6) は混ぜ合わせの重み Wp (Z 0 ; Z, N ) と Wn (N 0 ; Z, N ) を決める為に用いられ、回転での混ぜ合わせの重みは等式 (6.4) によって固有殻エネルギーを与える。
等式 (6.1) と (6.2) で明確に表される仮定は直接計算を使わない。それらに対する物理的洞察を与えることにより、等式 (6.4)-(6.6) の仮定を支えるためだけに等式 (6.5)、(6.6) を導入した。

6.0.5

球形核の重ね合わせとしての変形核

2、3 の仮定と近似に基ずく混ぜ合わせの重み Wp (Z 0 ; Z, N ) と Wn (N 0 ; Z, N ) を決める為に、規則を与える。図 (6.1)、(6.2) はこの規則の理解を容易にするだろう。この後で Wp (Z 0 ; Z, N ) を扱うが、Wn (N 0 ; Z, N ) も同様に扱うことができる。Wp (Z 0 ; Z, N )を決める為に用いれる基礎的等式はが等式 (6.6) である。
等式 6.6 を用いる時に、全ての原子核が空間的に sharp-cut な表面を持つ近似をする。実際に核の表面はぼやけているが、表面のぼやけの側面は核種と核の形状に余り依存する様には見えないから、上の近似はあまり極端ではない。変形核に対しては表面の内部を一様な密度と仮定する。これに対して後で詳細に説明される様に、球形核についての一様な密度を取る余地がある。一様でない密度について変形核は重ね合わせとして表される。この自由度なしに計算結果として最終的に得られる変形はかなり小さいだろう。もう１つのこの自由度の支えは後で与えられる。
もし表面の内部密度が関連する全ての核について一様であれば、等式 (6.6) の使用は60

rim(N’+1/2)

rim(N’-1/2)

図 6.2: 球形核の内に含まれる中間形状核の分解説明図。斜線部は、球形核 (N 00 , Z 0 ) の一部に対応する。

比較的容易である。それ故に球形核成分の一様でない密度の場合、r から rim へ半径座標を変え、あたかも全内部密度が一様に成る様に問題を扱うことができる為に、中間形状へ核形状を対応させて変換する。この変換は中性子と陽子共に共通である。
変形核が与えられる時、占有される立方角 Ωoc と表す。つまり中心から与えられる距
離 r について、半径 r の球の外側に広がる核表面の部分は立方角 Ωoc に広がる。したがって半径座標の変化は、{[Ωoc (r )drim = 1 − h 1 −
4π

]}

dr ,

(6.7)

によって作られる。r = 0 に於いて rim = 0 である。等式 6.7 では、h(< 1) は正符号の調節できるパラメータである。もし変形核の最小半径を Rmin と表したならば 、r ≤ Rminについて rim であり、r > Rmin について rim < r である。核の中間形式はオリジナルの変形核と同じ密度を持つと仮定される。そして元々範囲 dr にある核物質は中間形状の核で対応している範囲 drim に移動される。したがって、
2Ωim (rim )rim
drim = Ωoc (r )r 2 dr ,(6.8)

を持つ。Ωim は中間形状の核での占有される立体角である。この等式は Ωim とそれに応じて中間形状を決める。
座標 rim について仮定される重要な性質は、この座標系での変形核を形作る成分として用いられる球形核が変形核のものと同じ核子数密度 ρ̄ を持つ。球形核の成分に対して核物質保存の等式は、
2drim ,
ρ0 (r)r2 dr = ρ̄rim

(6.9)

2
drim であである。ρ0 (r) は半径 r の基の球形核の密度である。r > Rmin では r2 dr > rim
るから、同じ r の範囲内で ρ0 (r) < ρ̄ を持つ。一様でない密度を持つこれらの球形核は一般的にエネルギー固有状態に対応していない事に注目すべきである。一様でない球
形核の使用が変形核の形状において都合のいい理由は、変形核の全方向で平均される61ポテンシャルは通常の球形核のポテンシャルよりも浅いからである。もし 一粒子ポテ
ンシャルがもっと浅ければ 、一粒子の波動関数はより外側に広がるだろう。このような一粒子波動関数が球形核の最低エネルギーを与えないにも関わらず、中性子・陽子
殻エネルギー Ens (Z, N ), Eps (Z, N ) は２つの原子核に対して同様であると仮定される。
rim での座標系では、等式 6.6 はρnim (rim ; Z, N ) =∑Wn (N 0 ; Z, N )ρ0nim (rim ; Z 00 , N 0 ),(6.10)

N0として書かれる。ρnim (rim ; Z, N ) は全方向で平均される中間形状核の中性子数密度である。等式 (6.10) の中のもう１つの密度関数、すなわち ρ0nim (rim ; Z 00 , N 0 ) は rim 座標での球形核の密度を表し 、次の様な階段関数になる。

ρ0nim (rim ; Z 00 , N 0 ) =


 ρ̄
 0

n

rim ≤ R(N 0 )rim > R(N 0 )R(N 0 ) = (N 0 /N )1/3 R0 .

(6.11)(6.12)

等式 (6.11) の ρ̄n は中間の座標系での全ての原子核に於ける共通な中性子密度であり、
等式 (6.12) の R0 は中性子 N を持つ球形核の半径である。一般化すると等式 (6.12) によって、中性子数 N1 と半径座標 rim を関連付ける

rim (N1 ) = (N1 /N )1/3 R0 .

(6.13)

等式 (6.6)、(6.10) の中で N 0 は実際に整数しか取ることができないけれども、N1 だ
けでなく N 0 もどのよな実数でも取ることができると仮定されることによって近似的にそれらを解明する。したがって、∫

ρnim (rim (N1 ); Z, N ) =Wn (N 0 ; Z, N )ρ0nim (rim (N1 ); Z 00 , N 0 )dN 0 ,(6.14)

として等式 (6.10) を近似することができる。そして積分は十分に広い範囲でなされるべきである。等式 (6.14) を N1 について微分すると、∫
dρnim (rim (N1 ); Z, N)= − Wn (N 0 ; Z, N )ρ̄n δ(N1 − N 0 )dN 0
dN1= −ρ̄n Wn (N1 ; Z, N ),(6.15)を得る。密度関数 ρnim (rim (N1 ); Z, N ) は占有される立体角 Ωim (im) を用いて表される。

ρ̄n Ωim (rim (N1 )).
4π(6.16)

1 dΩim (rim (N1 ))
,
4π
dN1

(6.17)

ρnim (rim (N1 ); Z, N ) =
したがって等式 (6.15)、(6.16) から、

Wn (N1 ; Z, N ) = −を得る。陽子に対しても同様な方法で Wp (Z1 ; Z, N ) を得る。したがって等式 (6.4) により固有殻エネルギーを計算することができる。等式 (6.4) の中で、Z 0 と N 0 について62

の総和はしばしば 多い数の核種まで延ばす。それはとても広範囲な核種についての球
形陽子・中性子殻エネルギーを計算したのが理由である。変形核が陽子・中性子ド リッ
プライン近傍にある時、たまにド リップライン外側の球形殻エネルギーを必要とする。
そのような場合、同じ Z 0 又は N 0 を用いてド リップラインの丁度内側にある殻エネル
ギーを使う。例えば、もし中性子ド リップラインの外側にある等式 (6.4) の Ens (Z 00 , N 0 )を必要とする時、核種 (Z1 , N 0 ) が内側にあり、中性子ド リップラインに最も近いところにある Ens (Z1 , N 0 ) を代わりに用いる。
実際に Ein (Z, N ) の計算は以下の様に行われる。半径を方向の関数 (θ, φ) に与えることにより核の形状を記述する。

R = R(θ, φ; α, s).

(6.18)

α, s は変形を詳細に述べる為のパラメータを表す。α, s の値が 一旦固定され、変形が固定されると、占有される立体角 Ωoc (r ) が決まる。し たがって等式 (6.7)、(6.8) を使って,Ωim (rim ) を得る。Ωim (rim ) については等式 (6.17) によって混ぜ合わせの重み
Wn (N1 ; Z, N ) を与える (又は対応する等式によって Wp (Z1 ; Z, N ) を与える)。最後に固有殻エネルギー Ein (Z, N ) はこれらの混ぜ合わせの重みと球形殻エネルギーを用いて、等式 (6.4) から得られる。[19]

63

付録

Program List//Nuclear mass formula of deformed nuclei#include<stdio.h>#include<math.h>#include<stdlib.h>#define Ntheta 100#define NMAX 300#define ZMAX 200#define neutron1max NMAX
#define proton1max ZMAX
#define NMIN 1#define ZMIN 1#define EPS 1.0e-8#define KMAX 4#define JMAX 1#define CHECK 1#define NON 1.0e8
typedef enum{Spherical,Prolate,Oblate,Unknow} SHAPE;//-------------------------- Function prototype ---------------------------void Read_Dats();
int skyrmeS3_dat();int NE_mass03_dat();int Alpha2_dat();int LSM();
void gauss1(double a[KMAX][KMAX+1]);
void RMS1();
double term(int TN,int N,int Z);
int Make_nucl_chart();int extrapolation();int Nucl_shape(double alpha2);
double radius0(int N,int Z);
double N_R(double R);
SHAPE SetDeformation(double alpha2b,double alpha4b,double alpha6b,double R0b);
double radius(double theta);
int rim_theta_relation();double drdtheta(double theta);
double Omega(double N1_min,double N1_max);
double Omega_N1(double N1);
double rim_N1(double N1);
double theta_N1(double Rim_N1);
double interpolate(int N,double x0,double x[],double y[]);double radius_N1(double Theta_N1);
double Omega_im_prolate(double Rim,double Theta_N1,double R);

64

double Omega_im_oblate(double Rim,double Theta_N1,double R);
double GS_shell_energy();double SUM1_Wn_E();
double Ener_g(int N,int Z);
double Ener_def(int N,int Z);
double SHIFT_1(int N1,int Z1);
double SUM2_Wn_E();
double two_points1(int x,double y,int y1,int y2);
double two_points2(int x,double y,int y1,int y2);
double four_points(double N,int N1,int N2,double Z,int Z1,int Z2);
double Weight(double n0,int i,int imax,int h);double Nucl_Ener(double Es);
double Ener_eo(int N,int Z);
void RMS0();
void RMS2();
void RMS3();
double AW_g(int N,int Z);//-------------------------- Global variable -----------------------------SHAPE shape;
int proton;//Z

int neutron;
int Opion1;//N :neutron numbers//1 or 2

int Opion2;
double alpha2;//1 or 2 or 3//alpha2,alpha4 and alpha6

double alpha4;
double alpha6;//

double R0;
double lambda;//The radius of the shpere//Ddenominator lamdba

double thetamin;
double Rmin;//The minimum value of theta//The minimum radius od deformed nucleus

double Rmax;
double Rim_min;//The maximum radius od deformed nucleus//The minimum value of Rim

double Rim_max;
double hpara=0.46;//The maximum value of Rim//Positive adjustable parameter(hpara<1)

double r0=1.2;
double Theta[Ntheta+1];//The value of r0 is 1.2fm//The polar angle

double Rim[Ntheta+1];
double PI=M_PI;//Radius in the intermediary-shape nucleus//Pi(disambiguation)

double RimN[neutron1max+1];
double N1_min;//Values of Rim in N1//The minimum value of N1

double N1_max;//The maximum value of N1

double OmegaN_min;double OmegaN_max;//Value of Omega in N1_min//Value of Omega in N1_max:proton numbers:alpha2i are parameters to specify the shpa

double NE_sky[NMAX+1][ZMAX+1];//Nuclear Energy of Skyrme S3
double NE_sky_extr[NMAX+1][ZMAX+1]; //Nuclear Energy of Skyrme S3 + extrapolation
double alpha2_dat[NMAX+1][ZMAX+1];
double alpha2_MIN[NMAX+1][ZMAX+1];//Values of alpha2 data in neutron and proton//The minimum value of alpha2

double u =931.4939;
double Mn=939.5652;//atomic mass unit//neutron mass

double Mh=938.7829;
double kel=1.433*1.0e-5;//hydrogen atom mass//electronic constant of proportion(MeV)

double C_gross[KMAX];
double AW_gross[KMAX];
double AME[NMAX+1][ZMAX+1];
double NE_m03[NMAX+1][ZMAX+1];//atomic mass excess//nuclear mass by mass.mas03

65

double Fermi[ZMAX+1][NMAX+1][JMAX+1]; //Fermi level(n)int Use_NC[NMAX+1][ZMAX+1];//Data used for nuclear chart in data of Skyrme S3// Used data is 0 or Not used data is 1

double NE[NMAX+1][ZMAX+1];
double Nmin,Zmin;
double Omega1,Omega2;
double Esh_MIN;//-------------------------- Main function--------------------------------

int main(){Opion1=2;//Opion1=1: N1 is an integeral number//Opion1=2: N1 is an real number

if(Opion1!=1 && Opion1!=2) exit(1);Opion2=2;//Opion1=1: Nuclear chart//Opion1=2: Experimental valus of mass and Calculation valus of mass//Opion1=3: Data of alpha2 and Calculation value of the minimum alpha2

if(Opion2!=1 && Opion2!=2 && Opion2!=3) exit(1);switch(Opion2){
case 1:
Nmin=NMIN;
Zmin=ZMIN;
break;
case 2:
Nmin=8;
Zmin=8;break;
case 3:
Nmin=8;
Zmin=8;break;}
Read_Dats();
LSM();
Make_nucl_chart();if(Opion2==2){
RMS2();
RMS0();}else if(Opion2==3){
RMS3();}else

exit(1);}//---------------------- Reading of Data -----------------------------------void Read_Dats(){skyrmeS3_dat();
NE_mass03_dat();Alpha2_dat();}66

//-------------------- Substitute Data of Skyrme_S3 --------------------------int skyrmeS3_dat(){//input : Z,N//output:char buff[80];int n1,z1;
int n, z;
int i=1;
double NE;//Nuclear Energy of Skyrme S3
double FLN,FLZ;
int N=0,Z=1;
int Ndrip,Zdrip;
FILE *FSK;FSK=fopen("SIII.dat","r");//Data input

for(n1=0;n1<=NMAX;n1++){for(z1=0;z1<=ZMAX;z1++){
NE_sky[n1][z1]=NON;
NE_sky_extr[n1][z1]=NON;
Fermi[n1][z1][Z]=NON;
Fermi[n1][z1][N]=NON;}}
while(fgets(buff,80,FSK)!=NULL){
if(buff[0] == ’#’) continue; // Skip comment lines
if(sscanf(buff,"%d %d %lf %lf %lf",&z,&n,&NE,&FLZ,&FLN)==5){if(n>=0 && n<=NMAX && z>=0 && z<=ZMAX){
NE_sky[n][z]=NE;
NE_sky_extr[n][z]=NE;
Fermi[n][z][Z]=FLZ;
Fermi[n][z][N]=FLN;}
else if(n<0 || z<0){
break;}
else {
fprintf(stderr,"warning : data out of array size: N,Z=%d %d\n",n,z);}}}
if(0){ // interpolation to odd nuclei : begin
for(z1=ZMIN;z1<=ZMAX;z1+=2){
for(n1=NMIN;n1<NMAX;n1+=2){
if(NE_sky[n1+2][z1]<NON && NE_sky[n1][z1]<NON ){//proton:odd ,neutron:even
NE_sky[n1+1][z1]=(NE_sky[n1+2][z1]+NE_sky[n1][z1])/2.0;}}}
for(n1=NMIN;n1<=NMAX;n1+=2){
for(z1=ZMIN;z1<ZMAX;z1+=2){
if( NE_sky[n1][z1+2]<=NON && NE_sky[n1][z1]<=NON){ //proton:even,neutron:odd
NE_sky[n1][z1+1]=(NE_sky[n1][z1+2]+NE_sky[n1][z1])/2.0;
}

67

}
}
for(z1=2;z1<=ZMAX;z1+=2){//proton:odd , neutron:odd

for(n1=2;n1<=NMAX;n1+=2){
if((NE_sky[n1+1][z1]!=0.0)&&(NE_sky[n1][z1+1]!=0.0)&&(NE_sky[n1+2][z1+1]!=0.0)&&(NE_sky[n1+1][z1+2]!=0.0)){
NE_sky[n1+1][z1+1]=(NE_sky[n1+1][z1]+NE_sky[n1][z1+1]+NE_sky[n1+2][z1+1]+NE_sky[n1+1][z1+2])/4.0;}}}} // interpolation to odd nuclei : endfor(z1=ZMIN;z1<=ZMAX;z1++){
Ndrip=0;
for(n1=NMIN;n1<=NMAX;n1++){
if(NE_sky[n1][z1]<NON){
if(Fermi[n1][z1][N]<=0.0 && Fermi[n1][z1][Z]<=0.0){
Use_NC[n1][z1]=0;
if(n1>Ndrip) Ndrip=n1;}else if(Ndrip>0 && n1-Ndrip<=4){
Use_NC[n1][z1]=0;}else{
Use_NC[n1][z1]=(int)NON;}}else

Use_NC[n1][z1]=(int)NON;}}
for(n1=NMIN;n1<=NMAX;n1++){
Zdrip=0;
for(z1=ZMIN;z1<=ZMAX;z1++){
if(NE_sky[n1][z1]<NON){
if(Fermi[n1][z1][N]<=0.0 && Fermi[n1][z1][Z]<=0.0){
if(z1>Zdrip) Zdrip=z1;}else if(Zdrip>0 && z1-Zdrip<=4){
Use_NC[n1][z1]=0;}}}}
fclose(FSK);
return 0;}//---------------

Nuclear mass from Data by Audi-Wapstra ------------------

int NE_mass03_dat(){char buff[35];
int i=1;
int Z,N;
double ME;
FILE *FAM;//atomis mass excess

FAM=fopen("atom_03.dat","r");for(Z=0;Z<=ZMAX;Z++){

68

for(N=0;N<=NMAX;N++){AME[N][Z]=NON;
NE_m03[N][Z]=NON;}}
while(fgets(buff,35,FAM)!=NULL){
if(sscanf(buff,"%d %d %lf",&N,&Z,&ME)==3){
if(N>=0 && N<=NMAX && Z>=0 && Z<=ZMAX){
AME[N][Z]=ME;
i++;
}
}
}
fclose(FAM);
i=1;
for(Z=0;Z<=ZMAX;Z++){
for(N=0;N<=NMAX;N++){if(AME[N][Z]<NON){
if((ceil(AME[N][Z])-floor(AME[N][Z]))==1.0 || AME[N][Z]==0.0){i++;
NE_m03[N][Z]=AME[N][Z]*1.0e-3+(u-Mn)*N+(u-Mh)*Z+kel*(pow(Z,2.39)-(double)Z);}}}}
return 0;}//--------------------- alpha2 Data ------------------int Alpha2_dat(){char buff[20];int i=1;
int Z,N;
double alp2;
FILE *FA;
FA=fopen("a20frdm.dat","r");for(Z=0;Z<=ZMAX;Z++){
for(N=0;N<=NMAX;N++){alpha2_dat[N][Z]=NON;}}
while(fgets(buff,20,FA)!=NULL){
if(sscanf(buff,"%d %d %lf",&Z,&N,&alp2)==3){
if(N>=0 && N<=NMAX && Z>=0 && Z<=ZMAX){alpha2_dat[N][Z]=alp2;
i++;
}
}
}
fclose(FA);return 0;69}//------------------ least squares method --------------------int LSM(){int l;
int k;//Matrix A of file//Matrix A of rank

int N,Z;double s1,s2;//Matrix A of components

double t1,t2;
int n=1;//Matrix A of 5th component

double A1[KMAX][KMAX+1];
double A2[KMAX][KMAX+1];//Matrix A//Matrix A

double a;for(l=0;l<KMAX;l++){for(k=0;k<KMAX;k++){s1=0.0;t1=0.0;s2=0.0;t2=0.0;for(N=Nmin;N<=NMAX;N++){
for(Z=Zmin;Z<=ZMAX;Z++){switch(Opion2){
case 1:a=(double)Use_NC[N][Z];break;
case 2:
a=NE_m03[N][Z];break;
case 3:a=alpha2_dat[N][Z];break;}
if(Opion2==1){
if(a<NON){
s1+=term(k,N,Z)*term(l,N,Z);
t1+=term(l,N,Z)*NE_sky[N][Z];n++;}}else if(Opion2==2){
if(a<NON && N%2==0 && Z%2==0){
s1+=term(k,N,Z)*term(l,N,Z);
t1+=term(l,N,Z)*NE_sky[N][Z];s2+=term(k,N,Z)*term(l,N,Z);
t2+=term(l,N,Z)*NE_m03[N][Z];n++;}}else if(Opion2==3){
if(a<NON && N%2==0 && Z%2==0){
s1+=term(k,N,Z)*term(l,N,Z);
t1+=term(l,N,Z)*NE_sky[N][Z];
n++;
}
}
}
}

70

n=1;
A1[l][k]=s1;
if(Opion2==2) A2[l][k]=s2;
A2[l][k]=s;}A1[l][KMAX]=t1;
if(Opion2==2) A2[l][KMAX]=t2;
A2[l][KMAX]=t2;}
gauss1(A1);gauss_j(A2);
return 0;}//-------------------- Gauss elimination part1 ---------------------void gauss1(double a[KMAX][KMAX+1]){//input : A[][]//output: X[]
int i,j,k;double t;
double x[KMAX];for(k=0;k<KMAX-1;k++){for(i=k+1;i<KMAX;i++){t=a[i][k]/a[k][k];for(j=k+1;j<=KMAX;j++){a[i][j]-=t*a[k][j];}}}for(i=KMAX-1;i>=0;i--){t=a[i][KMAX];for(j=i+1;j<KMAX;j++){t-=a[i][j]*a[j][KMAX];}a[i][KMAX]=t/a[i][i];C_gross[i]=a[i][KMAX];}/*行列が最後ど うなったか見たいときに実行*/#if CHECK==1
for(i=0;i<KMAX;i++) {for(j=0;j<KMAX+1;j++) {}}#endiffor(i=0;i<KMAX;i++){}
RMS1();}//-------------------------------------------int gauss_j(double a[KMAX][KMAX+1]){71

int k; /* pivot */int i; /* row */int j; /* column */double pivot, w;
int k1;
double x[KMAX];for(k=0; k<KMAX; k++) {pivot = a[k][k];if(fabs(pivot) < EPS) {printf("PIVOT is nearly equal zero.\n");return 1;}
k1 = k + 1;
for(j=k1; j<KMAX+1; j++) {a[k][j] = a[k][j] / pivot;}
for(i=0; i<KMAX; i++) {if(i != k) {w = a[i][k];for(j=k1; j<KMAX+1; j++) {a[i][j] = a[i][j] - w * a[k][j];}}}for(i=0; i<KMAX; i++) {for(j=0; j<KMAX+1; j++) {}//printf("\n");}}for(i=0;i<KMAX;i++) {AW_gross[i]=a[i][KMAX];}
RMS1();
return 0;}//------------------ Root Mean Square 1 -------------------------void RMS1(){//input: X[KMAX],N,Z,BE[][]//output:RMS
int A;
int l;
int n;
int N,Z;double EV;//emprical value

double CV;
double S=0.0;//calaulated value

double RMS;//Root Mean Squaren=1;
for(N=Nmin;N<=NMAX;N++){72

for(Z=Zmin;Z<=ZMAX;Z++){
EV=0.0;CV=0.0;switch(Opion2){
case 1:A=(double)Use_NC[N][Z];break;
case 2:
A=NE_m03[N][Z];break;
case 3:
A=alpha2_dat[N][Z];break;}
if(Opion2==1){
if(A<NON){
for(l=0;l<KMAX;l++){CV+=term(l,N,Z)*C_gross[l];}EV=NE_sky[N][Z];S+=pow((CV-EV),2.0);n++;}}else if(Opion2==2 || Opion2==3){
if(A<NON){
for(l=0;l<KMAX;l++){CV+=term(l,N,Z)*C_gross[l];}EV=NE_sky[N][Z];S+=pow((CV-EV),2.0);
n++;
}
}
}
}
RMS=sqrt(S/(double)n);}//----------------- term variable ---------------------double term(int TN,int N,int Z){//input : TN(term number),N,Z//output: f
int A;
double f;//mass number

A=N+Z;
if(TN==0)
f=(double)A;else if(TN==1) f=pow(A,2.0/3.0);
else if(TN==2) f=(pow((N-Z),2.0))/(double)A;else if(TN==3) f=(pow(Z,2.0))*(pow(A,-1.0/3.0));
else

exit(1);

return f;}//--------------------------------------------------------------------------

73

int Make_nucl_chart(){int i;
int I= 600;
int n,z;//Number of alpha2 that can be put in the direction of positive(the negative direction)//n is neutron number

double Es;
double P;
int A;
int imax;//alpha2 maximun * 1e3

FILE *fp;
fp=fopen("MF_SIII_D34.dat","w");extrapolation();if(0){
for(z=ZMIN;z<=ZMAX;z++){for(n=NMIN;n<=NMAX;n++){if(NE_sky_extr[n][z]<NON){
printf("%d %d %f\n",n,z,NE_sky_extr[n][z]);}}}}for(z=Zmin;z<=ZMAX;z++){proton=z;for(n=Nmin;n<=NMAX;n++){neutron=n;A=n+z;
imax=I-A;
Esh_MIN=0.0;
switch(Opion2){
case 1:P=(double)Use_NC[n][z];break;
case 2:
P=NE_m03[n][z];break;
case 3:
P=alpha2_dat[n][z];break;}
if(Opion2==1){
if(P<NON){
for(i=-imax;i<=imax;i++){if(i!=0){alpha2=i*1.0e-3;
Nucl_shape(alpha2);
rim_theta_relation();N1_min=N_R(Rim_min);N1_max=N_R(Rim_max);
Es=GS_shell_energy();}}NE[neutron][proton]=Nucl_Ener(Es);74}}else if(Opion2==2 || Opion2==3){
if(P<NON){
for(i=-imax;i<=imax;i++){if(i!=0){alpha2=i*1.0e-3;
Nucl_shape(alpha2);
rim_theta_relation();N1_min=N_R(Rim_min);N1_max=N_R(Rim_max);
Es=GS_shell_energy();}}NE[neutron][proton]=Nucl_Ener(Es);fprintf(fp,"%4d %4d %12.6f %12.6f %12.6f %12.6f %9.5f\n",z,n,NE[neutron][proton],NE_sky[neutron][proton],NE_m03[neutron][proton],AW_g(neutron,proton),alpha2_MIN[neutron][proton]);printf("%4d %4d %12.6f %12.6f %12.6f %12.6f %9.5f\n",z,n,NE[neutron][proton],NE_sky[neutron][proton],NE_m03[neutron][proton],AW_g(neutron,proton),alpha2_MIN[neutron][proton]);}}}}fclose(fp);}//------------------ Extrapolation of Data ------------------------int extrapolation(){int i,j,k,h,m,n;
int nmax;
int Z,N;
double e[8];double E1,E2,E3,E4;
nmax=1e5;for(n=0;n<=nmax;n++){m=0;for(Z=ZMIN;Z<=ZMAX;Z++){
for(N=NMIN;N<=NMAX;N++){if(NON-NE_sky_extr[N][Z]<=EPS){
k=0;
h=0;
for(i=-1;i<=1;i++){for(j=-1;j<=1;j++){if(i!=0 || j!=0){
if(N+j<=NMIN || N+j>=NMAX || Z+i<=ZMIN || Z+i>=ZMAX)
else
k++;

e[k]=0.0;e[k]=NE_sky_extr[N+j][Z+i];}}}if(e[1]<0.0 && e[3]<0.0 && e[0]<0.0){E1=e[1]*e[3]/e[0];
h++;

75

}else{E1=0.0;}if(e[1]<0.0 && e[4]<0.0 && e[2]<0.0){E2=e[1]*e[4]/e[2];h++;}else{
E2=0.0;}if(e[6]<0.0 && e[4]<0.0 && e[7]<0.0){E3=e[6]*e[4]/e[7];k++;}else{
E3=0.0;}if(e[6]<0.0 && e[3]<0.0 && e[5]<0.0){E4=e[6]*e[3]/e[5];h++;}else{
E4=0.0;}if(h>0){NE_sky_extr[N][Z]=(E1+E2+E3+E4)/((double)h);
m++;
}
}
}
}
if(m==0) break;}}//--------------------------- Set alpha2 ------------------------------int Nucl_shape(double alpha2){//input : alpha2b,alpha4b,alpha6b,R0b,neutron,Rim_min,Rim_max//output: Sum_W_N
int n1;//neutron1

int ret;int i;//Return value from function SetDeformation

double alpha2b;
double alpha4b;
double alpha6b;
double R0b;alpha2b=alpha2;alpha4b=0.0;
alpha6b=0.0;
R0b=radius0(neutron,proton);
ret=SetDeformation(alpha2b,alpha4b,alpha6b,R0b);if(ret==Spherical);//printf("Shape is Spherical.\n");else if(ret==Oblate) ;//printf("Shape is Oblate.\n");else if(ret==Prolate) ;//printf("Shape is Prolate.\n");else if(ret==Unknow){
fprintf(stderr,"error:Shape is Unknow.\n");76

exit(1);}return 0;}//----------------------

Radius of the sphere -------------------------

double radius0(int N,int Z){//input : N,Z//output: R0
return r0*pow((double)(N+Z),1.0/3.0);}//-------------- Neutron1_min or Neutron1_max --------------double N_R(double R){//input : R,R0//output: (double)NZ*pow(R/R0,3.0)
return (double)neutron*pow(R/R0,3.0);}//---------------------- Deformation --------------------------SHAPE SetDeformation(double alpha2b,double alpha4b,double alpha6b,double R0b){//input : alpha2b, alpha4b, alpha6b, R0b, Rmin, Rmax, shape//output: (global variables)lambda, thetamin, Rmin//(return value) shape

int i;
double r1,r2,rmin;
double l;double theta;
int nincr=0;
int ndecr=0;
alpha2=alpha2b;alpha4=alpha4b;alpha6=alpha6b;
R0=R0b;l=1.0+(3.0/ 5.0)*alpha2*alpha2+(3.0/ 9.0)*alpha4*alpha4+(3.0/13.0)*alpha6*alpha6;
if(l<0){fprintf(stderr,"error:l=%f\n",l);exit(1);}else{
lambda=pow(l,1.0/3.0);}
for(i=0;i<=Ntheta;i++){theta=M_PI/(Ntheta*2.0)*i;
r1=radius(theta);
if(i==0){rmin=r1;
thetamin=theta;}else{
if(r1<rmin){rmin=r1;
thetamin=theta;77}if(r1<r2)else if(r1>r2)

ndecr++;
nincr++;}
r2=r1;}
Rmin=rmin;
Rim_min=Rmin;if(ndecr==0){
if(nincr==0)
else

shape=Spherical;
shape=Oblate;}else{
if(nincr==0)

shape=Prolate;else

shape=Unknow;}
return shape;}//----------------- Radii as a function of theta --------------------double radius(double theta){//input : theta,lambda, R0//output: radius
double r;
double P2,P4,P6;//Legendre polynomials

double c,c2,c4,c6;
c=cos(theta);c2=c*c;c4=c2*c2;
c6=c4*c2;
P2=
P4=(3.0*c2-1.0)/ 2.0;
(35.0*c4 -30.0*c2+3.0)/ 8.0;

P6=(231.0*c6-315.0*c4+105.0*c2-5.0)/16.0;
r=R0*(1.0+alpha2*P2+alpha4*P4+alpha6*P6)/lambda;
static int ncall=0;
ncall++;
if(ncall%100==1){//fprintf(stderr,"radius:%f %f %f %f %f %f\n",alpha2,alpha4,alpha6,lambda,theta,r);}
return r;}//------------------------- Rim of theta----------------------------------

int rim_theta_relation(){//input : Ntheta, hpara, thetamin, Rmin//output: Rim[]
int i,j;double hj;
double sum;double s;
double R;
double theta;
double Ntheta2=100;78

double R1;
double rim_max;
double COS;
for(i=0;i<=Ntheta;i++){
Theta[i]=(M_PI*i)/(2*Ntheta);}
Rim[0]=0.0;
for(i=0;i<Ntheta;i++){sum=0.0;hj=(Theta[i+1]-Theta[i])/Ntheta2;
for(j=0;j<=Ntheta2;j++){theta=hj*j+Theta[i];COS=cos(theta);
if(Prolate==shape){//prolate shapeif(j==0||j==Ntheta2) sum+=(1.0-hpara*COS)*drdtheta(theta)/2;
else
sum+=(1.0-hpara*COS)*drdtheta(theta);}else if(Oblate==shape){//oblate shapeif(j==0||j==Ntheta2) sum+=(1.0-hpara*(1.0-COS))*drdtheta(theta)/2;
else}else exit(1);sum+=(1.0-hpara*(1.0-COS))*drdtheta(theta);}
sum*=hj;
Rim[i+1]=sum+Rim[i];}if(fabs(thetamin-0.0)<EPS)
else if(fabs(thetamin-M_PI*0.5)<EPS)

s=Rmin;s=Rmin-Rim[Ntheta];

else{
fprintf(stderr,"error:thetamin=%f\n",thetamin);exit(1);}for(i=0;i<=Ntheta;i++){Rim[i]+=s;}for(i=0;i<=Ntheta;i++){
R1=Rim[i];if(i==0){
rim_max=R1;}else{
if(R1>rim_max)

rim_max=R1;}}
Rim_max=rim_max;
return 0;}//---------------------- dr divied by dtheta ----------------------------------double drdtheta(double theta){//input : theta,alpha2,alpha4,alpha6,lambda,R0//output: drdthe
double drdthe;
double s,c,c2,c4;
s=sin(theta);

79c=cos(theta);c2=c*c;c4=c2*c2;
drdthe=(R0*s*c*(-24.0*alpha2+60.0*alpha4-105.0*alpha6+(-140.0*alpha4+630.0*alpha6)*c2-693.0*alpha6*c4))/(8.0*lambda);
static int ncall=0;
ncall++;
if(ncall%100==1){//fprintf(stderr,"drdtheta:%f %f %f %f %f %f\n",alpha2,alpha4,alpha6,lambda,theta,drdthe);}return drdthe;}//------------------------

Omega_im ---------------------------

double Omega(double N1_min,double N1_max){//input : N1_min,N1_max//output: Omega
int i;
int min;int max;
double Omega;min=ceil(N1_min);max=floor(N1_max);
Omega=Omega_N1(N1_min);for(i=min;i<=max;i++){
Omega=Omega_N1((double)i);}
Omega=Omega_N1(N1_max);
return 0;}//---------------------------- Omega(N1) -------------------------------double Omega_N1(double N1){//input : N1//output: Omega
double Rim_N1;
double Theta_N1;
double R_N1;
double Omega;
Rim_N1=rim_N1(N1);
Theta_N1=theta_N1(Rim_N1);
R_N1=radius_N1(Theta_N1);
if(shape==Prolate)else if(shape==Oblate)

Omega=Omega_im_prolate(Rim_N1,Theta_N1,R_N1);
Omega=Omega_im_oblate (Rim_N1,Theta_N1,R_N1);

else

exit(1);

return Omega;}80

//-------------- Rim of neutron1 or proton1(thr integral number) ----------double rim_N1(double N1){//input:N1,R0
return R0*pow(N1/(double)neutron,1.0/3.0);}//-------------------------

Theta_N1 -----------------------------

double theta_N1(double Rim_N1){//input: Rim_N1//(used global variables):Ntheta,Rim,Rim[],Theta[]
return interpolate(Ntheta,Rim_N1,Rim,Theta);}//--------------- Lagrange’s polynomial interpolation formula ---------------double interpolate(int N,double x0,double x[],double y[]){//input : N,x0,x[],y[]//output: interpolate
int i,j,k,n,m;
int i1,i2;
double p,q;//i1=i-3, i2=i+2

i=0;
i1=i-3;i2=i+2;if(x[0]<x[N]){if(x0-x[0]<1.0e-18&&x[N]-x0<1.0e-18){
return 0;}else{
for(n=0;n<=N;n++){if(x[i]<=x0){
i=n;}else break;}if(i<3)m=3-i;

else if(i>N-2)
else

m=N-2-i;m=0;

i1=i-3+m;i2=i+2+m;p=0;
for(k=i1;k<=i2;k++){q=1;
for(j=i1;j<=i2;j++){if(j!=k) q=q*(x0-x[j])/(x[k]-x[j]);}p=p+y[k]*q;}}}else if(x[N]<x[0]){if(x0-x[N]<1.0e-18&&x[0]-x0<1.0e-18){
return 0;}else{
for(n=0;n<=N;n++){if(x[N-i]<=x0){

81

i=n;
}else

break;}
if(i<2)else if(i>N-3)

m=-3+i;m=N+2-i;else
i1=i-3-m;m=0;i2=i+2-m;p=0;
for(k=i1;k<=i2;k++){q=1;
for(j=i1;j<=i2;j++){if(j!=k) q=q*(x0-x[N-j])/(x[N-k]-x[N-j]);}p=p+y[N-k]*q;}}}else{
printf("error:3\n");exit(3);}return p;}//-------------------------- R of theta ---------------------------------double radius_N1(double Theta_N1){//input : Theta_N1//output: radius(Theta_N1)
return radius(Theta_N1);}//------------- The occupied solid angle of prolate shape ------------double Omega_im_prolate(double Rim,double theta,double R){//input : Rim,Theta_N1,R,hpara//output: Omega_P
int i;
double c;double hpc;
double Omega_P;
c=cos(theta);
Omega_P=((1.0-c)*R*R)/((1.0-hpara*c)*Rim*Rim);
return Omega_P;}//------------- The occupied solid angle of oblate shape -------------double Omega_im_oblate(double Rim,double theta,double R){//input : Rim,Theta_N1,R,hpara//output: Omega_O
int i;
double c;double hpc;
double Omega_O;82

c=cos(theta);
Omega_O=(c*R*R)/((1.0-hpara*(1.0-c))*Rim*Rim);
return Omega_O;}//------------------ The ground-state shell_energy --------------------double GS_shell_energy(){//input : neutron,proton,NE,skyrme[][],WeightN[]//output: Emf,E_weight
int N1MIN, N1MAX;
double Edef;
double Esh;double dEs;
double Emf;//Energy of Spherical Relativistc Mean Field model(SRLF)

Emf=NE_sky[neutron][proton];switch(Opion1){
case 1:
dEs=SUM1_Wn_E();
Edef=dEs+Ener_def(neutron,proton);
break;
case 2:
dEs=SUM2_Wn_E();
Edef=dEs+Ener_def(neutron,proton);
break;}
Esh=Emf+Edef;
if(Esh_MIN>Esh){alpha2_MIN[neutron][proton]=alpha2;
Esh_MIN=Esh;}return Esh_MIN;}//----------------------------------------------------------------------double SUM1_Wn_E(){
int i,imax;
int N1MIN,N1MAX;
int Z1f, Z1c;
int N1;
double Z1;
double Emf_n1z1;//Energy of SRLFin N1 and Z1

double Eg;
double E=0.0;//Energy that removes shell effect from SRLF//Sum total of product of WeightN and Emf_n1z1 - Eg

double Wn;double sum_Wn=0.0;
N1MIN=ceil(N1_min);
N1MAX=floor(N1_max);if(N1MIN-N1_min>0.5)

N1MIN=N1MIN-1;

if(N1_max-N1MAX>0.5)

N1MAX=N1MAX+1;

83imax=N1MAX-N1MIN;
for(N1=N1MIN;N1<=N1MAX;N1++){Z1=proton*(N1/(double)neutron);
Z1f=floor(Z1);Z1c=ceil(Z1);if(Z1c-Z1f==0){
if(NE_sky_extr[N1][Z1f]<NON && NE_sky_extr[N1][Z1c]< NON){
Emf_n1z1=NE_sky_extr[N1][Z1c];}else{//When both of the two values exist//When both of the two values do not exist

Emf_n1z1=NON;}
Eg=Ener_g(N1,Z1c);}else if(Z1c-Z1f==1){//Z1 is a real number

if(NE_sky_extr[N1][Z1f]<NON && NE_sky_extr[N1][Z1c]<NON){
Emf_n1z1=two_points1(N1,Z1,Z1f,Z1c);}else{
Emf_n1z1=NON;//When both of the two values exist//When both of the two values do not exist}
Eg=(Z1c-Z1)*Ener_g(N1,Z1f)+(Z1-Z1f)*Ener_g(N1,Z1c);}else{
exit(1);}i=N1-N1MIN;
Wn=Weight((double)N1MIN,i,imax,1);
sum_Wn+=Wn;if(Emf_n1z1>EPS)
else

E+=NON;//When the value of Emf_n1z1 doesn’t exist
E+=Wn*((Emf_n1z1-NE_sky_extr[neutron][proton])-(Eg-Ener_g(neutron,proton)));}return E;}//----------------- Nuclear mass formula of Weizsacker ---------------------double Ener_g(int N,int Z){//input : N,Z
int A;//output: E//mass number

double f0,f1,f2,f3;
double E;
A=N+Z;
f0=A;
f1=pow(A,2.0/3.0);f2=(pow((N-Z),2.0))/A;
f3=(pow(Z,2.0))*(pow(A,-1.0/3.0));
E=C_gross[0]*f0+C_gross[1]*f1+C_gross[2]*f2+C_gross[3]*f3;
return E;}//---------------- The average defotmation energy --------------double Ener_def(int N,int Z){//input : N, Z//output: E_def

84int A;//Mass number

double E_def;double E_s;//The average defotmation energy//The surface energy

double E_C;
double dE_prl;//Coloumb energy//The energy tofavor prolate shapes

double dE_s;
double dE_C;
double A2,A4,A6;
double a_s = 17.3*(1.0-0.34); //(MeV)
double a_sI= 5.0;//(MeV)
double a_C = C_gross[3];double C_prl1=0.28;//(MeV)

double C_prl2=5.0;
A=N+Z;
A2=alpha2*alpha2;
A4=alpha4*alpha4;
A6=alpha6*alpha6;
E_s=a_s*pow(A,2.0/3.0)-a_sI*pow((N-Z),2.0)*pow(A,-4.0/3.0);
E_C=a_C*Z*Z*pow(A,-1.0/3.0);
dE_s= E_s*(0.4*A2+

A4+(20.0/ 13.0)*A6);

dE_C=-E_C*(0.2*A2+(5.0/27.0)*A4+(25.0/169.0)*A6);
dE_prl=-C_prl1*alpha2*pow(A,2.0/3.0)*exp(-C_prl2*A2);
E_def=dE_s+dE_C+dE_prl;
return E_def;}//----------------------------------------double SUM2_Wn_E(){
int N1MIN,N1MAX;
int i,j;int imax;
int h=100;
double z,n;double z0;
int zf,zc,nf,nc;
double dEmf;
double dEg;//Energy obtained from Substitute Data of Skyrme_S3//Energy obtained from Nuclear mass formula of Weizsacker

double e1,e2;
double Wn;double dEs;
double sum_dEs=0.0;//Value of product of Wn and dEmf - dEg//Sum total of product of Wn and dEmf - dEg

double sum_Wn=0.0;
N1MIN=ceil(N1_min*h);
N1MAX=floor(N1_max*h);if((double)N1MIN-N1_min*h>0.5)

N1MIN=N1MIN-1;

if(N1_max*h-(double)N1MAX>0.5)

N1MAX=N1MAX+1;

85

imax=N1MAX-N1MIN;
for(i=0;i<=imax;i++){n=(N1MIN+i)/((double)h);z=proton*(n/(double)neutron);
zf=floor(z);zc=ceil(z);nf=floor(n);nc=ceil(n);dEg=Ener_g(n,z);
Wn=Weight(N1MIN/(double)h,i,imax,h);sum_Wn+=Wn;if(nc-nf==0 && zc-zf==0 && NE_sky_extr[nf][zf]<NON){
dEmf=NE_sky_extr[nf][zf];dEg=Ener_g(nf,zf);}else if(nc-nf==0 && zc-zf==1 && NE_sky_extr[nf][zf]<NON && NE_sky_extr[nf][zc]<NON){
dEmf=two_points1((int)n,z,zf,zc);
dEg =(zc-z)*Ener_g(n,zf)+(z-zf)*Ener_g(n,zc);}else if(nc-nf==1 && zc-zf==0 && NE_sky_extr[nf][zf]<NON && NE_sky_extr[nc][zf]<NON){
dEmf=two_points2((int)z,n,nf,nc);dEg =(nc-n)*Ener_g(nf,z)+(n-nf)*Ener_g(nc,z);}else if(nc-nf==1 && zc-zf==1 && NE_sky_extr[nf][zf]<NON && NE_sky_extr[nf][zc]<NON&& NE_sky_extr[nc][zf]<NON && NE_sky_extr[nc][zc]<NON){
dEmf=four_points(n,nf,nc,z,zf,zc);
e1 =((double)nc-n)*Ener_g(nf,zf)+(n-(double)nf)*Ener_g(nc,zf);
e2 =((double)nc-n)*Ener_g(nf,zc)+(n-(double)nf)*Ener_g(nc,zc);
dEg =((double)zc-z)*e1+(z-(double)zf)*e2;}else{
sum_dEs+=NON;}
dEs=Wn*((dEmf-NE_sky_extr[neutron][proton])-(dEg-Ener_g(neutron,proton)));sum_dEs+=dEs;}
return sum_dEs;}//------------------------------------------------------------------double two_points1(int x,double y,int y1,int y2){//input : N, N1, N2, Z//output: E
return ((double)y2-y)*NE_sky_extr[x][y1]+(y-(double)y1)*NE_sky_extr[x][y2];}//------------------------------------------------------------------double two_points2(int x,double y,int y1,int y2){//input : N, N1, N2, Z//output: E
return ((double)y2-y)*NE_sky_extr[y1][x]+(y-(double)y1)*NE_sky_extr[y2][x];}//-------------------------------------------------------------------double four_points(double N,int N1,int N2,double Z,int Z1,int Z2){//input : N, N1, N2, Z, Z1, Z2//output: E
double EZ1,EZ2;86

double E;
double n1,n2;
n1=N-(double)N1;n2=(double)N2-N;
EZ1=n2*NE_sky_extr[N1][Z1]+n1*NE_sky_extr[N2][Z1];EZ2=n2*NE_sky_extr[N1][Z2]+n1*NE_sky_extr[N2][Z2];E =((double)Z2-Z)*EZ1+(Z-(double)Z1)*EZ2;
return E;}//---------------------- Weight(N1) --------------------------double Weight(double n0,int i,int imax,int h){double n;double Weight;
if(i==0){n=n0+((double)i+0.5)/(double)h;Omega1=Omega_N1(N1_min);Omega2=Omega_N1(n);}else if(i==imax){
Omega1=Omega2;
Omega2=Omega_N1(N1_max);}else{
n=n0+((double)i+0.5)/(double)h;Omega1=Omega2;
Omega2=Omega_N1(n);}Weight=(Omega1-Omega2);
return Weight;}//----------------- Nunulear Mass --------------------------double Nucl_Ener(double Es){
returnEs+Ener_eo(neutron,proton)-kel*pow(proton,2.39);}//------------------ The even-odd term ------------------double Ener_eo(int N,int Z){//input : N, Z//output: E_eo
int A;
double E_eo;//Mass number//The even-odd term

double M_oo;double del_oddZ;
double del_oddN;//Parametersin the even-odd term
double C_eo1=
double C_eo2=

13.90;
0.94;//(MeV)//(MeV)

double C_eo3=
double C_eo4=

6.50;
2.00;

//(MeV)87

double C_oo1=

55.0;

double C_oo2=-120.0;
double C_oo3= 65.0;
A=N+Z;
if((Z%2)==1) del_oddZ=1;
else
del_oddZ=0;//for odd Z//for even Z

if((N%2)==1) del_oddN=1;//for odd

else//for even N

del_oddN=0;

N

M_oo=C_oo1*pow(A,-1.0)+C_oo2*pow(A,-4.0/3.0)+C_oo3*pow(A,-5.0/3.0);
E_eo=C_eo1*pow((1.0+C_eo2*exp(C_eo3*(N-2*Z+C_eo4)/(double)A)),-0.5)*pow(A,-0.5)*(del_oddZ+del_oddN)-M_oo*del_oddZ*del_oddN;
return E_eo;}//------------------ Root Mean Square 2 -------------------------void RMS0(){//input: NE_m03[N][Z],NE_sky[N][Z]//output:RMS
int i;
int N,Z;double S=0.0;
double RMS;//Root Mean Square

i=0;
for(N=Nmin;N<=NMAX;N++){
for(Z=Zmin;Z<=ZMAX;Z++){
if(NE_m03[N][Z]<NON && NE_sky[N][Z]<NON &&N%2==0 && Z%2==0){
S+=pow((NE_m03[N][Z]-NE_sky[N][Z]),2.0);
i++;
}
}
}
RMS=sqrt(S/(double)i);
printf("Skyrme and Audi-Wapsta : RMS=%10.8f(MeV) %5d\n",RMS,i);}//------------------ Root Mean Square 3 -------------------------void RMS2(){//input: NE[N][Z],NE_m03[N][Z]//output:RMS
int i;
int N,Z;double S=0.0;
double RMS;//Root Mean Square

i=0;
for(N=Nmin;N<=NMAX;N++){
for(Z=Zmin;Z<=ZMAX;Z++){
if(NE[N][Z]<NON && NE_m03[N][Z]<NON && N%2==0 && Z%2==0){
S+=pow((NE[N][Z]-NE_m03[N][Z]),2.0);//printf("%d %d %f %f %f\n",N,Z,CV,EV,fabs(CV-EV));
i++;

88

}
}
}
RMS=sqrt(S/(double)i);
printf("This work and Audi-Wapsta : RMS=%10.8f(MeV) %5d\n",RMS,i);}//------------------ Root Mean Square 4 -------------------------void RMS3(){//input: alpha2_MIN[N][Z]<NON,alpha2_dat[N][Z]//output:RMS
int i;
int N,Z;double S=0.0;
double RMS;//Root Mean Square

i=0;
for(N=Nmin;N<=NMAX;N++){
for(Z=Zmin;Z<=ZMAX;Z++){if(alpha2_MIN[N][Z]<NON && alpha2_dat[N][Z]<NON && N%2==0 && Z%2==0){
S+=pow((alpha2_MIN[N][Z]-alpha2_dat[N][Z]),2.0);
i++;
}
}
}
RMS=sqrt(S/(double)i);
printf("alpha2 : RMS=%10.8f(MeV) %5d\n",RMS,i);}//----------------- Nuclear mass formula of Weizsacker ---------------------double AW_g(int N,int Z){//input : N,Zdouble A;//output: E//mass number

double f0,f1,f2,f3;
double E;
A=(double)N+(double)Z;
f0=A;
f1=pow(A,2.0/3.0);f2=(pow((N-Z),2.0))/A;
f3=(pow(Z,2.0))*(pow(A,-1.0/3.0));
E=C_gross[0]*f0+C_gross[1]*f1+C_gross[2]*f2+C_gross[3]*f3;
return E;}

89