原子核の相対論的平均場模型プログラムの開発と検証
2008 年 2 月
18 年度入学 10 番 三和 之浩
原子核の相対論的平均場模型とは、核子 (陽子・中性子の総称) の自由度に加えて、核子間に働く相互作用の媒介として複数種の中間子場を陽に取り入れた模型であり、相対論的な現象論的手法として注目を集めている模型である。陽子間のクーロン斥力が陽子の密度分布を源とする Poisson 方程式の解として定義される電位ポテンシャルで表されるのと同様に、２核子の間に働く各スピン・アイソスピンチャンネルの相互作用
が 、それぞれ適当な種類の核子密度を源とする遮蔽効果のある Poisson 方程式の解として定義される σ 、ω 、ρ などの中間子場の振幅に比例したポテンシャルで表されるのである。
本論文では、原子核の相対論的平均場模型の解を、密度分布に球対称性を課して求める数値計算プログラムを開発し 、模型に関する古典的な文献のひとつに掲載された種々の計算結果との比較をした。そして、両者の違いが十分に小さいことから、プログラムの検証に成功したと結論づけた。
原子核の相対論的平均場模型における中間子場は、古典場扱いすることを前提にしてその性質 (中間子１個の質量や核子との結合定数などのパラメータの値) を調節された有効場である。量子場の理論からこの平均場模型を導くためには２つの近似が必要である。第一の近似は核子間の相互作用の媒介として導入された各種の中間子の場のもつ量子力学的揺動を無視してその平均値だけに注目するという、中間子場の古典場近似である。第二の近似は、核子がスピン 21 の粒子の従う相対論的な波動方程式である Dirac 方程式に従うとする場合に現れる負のエネルギーの解を無視するという所謂no-sea 近似である。本論文ではこれらの近似について簡潔な解説を試みた。
本論文ではさらに、原子核の相対論的平均場模型の簡潔にしてかつ完全な定式化を与えた。密度分布に球対称性を課す場合のプログラムを開発するという目的、また、開発されたプログラムを理解するという目的に必要な情報は完全に含まれている。
また、数値解法上の重要ポイントも解説した。特に、r = 0 および r → ∞ での解析
解への接続と、両端から延長した２つの解を連続に接続する為に、エネルギー固有値を探索する効率的な手順を論じた。特に 、２つの解を接続する地点 (matching point)の選択において、非相対論的に言えば 、エネルギーがポテンシャルエネルギーより低くなる地点を選ぶと、計算機上の実数表現の有限性のため、可能な最大精度までエネ
ルギーを調節しても matching point に他の領域を圧倒する高さをもったピークを伴った不正確な解しか得ることができないことに注意を喚起した。また、自己無撞着解を求める際の出発点とする核子のポテンシャルの初期値の良い設定法も提案した。
さらに付録として、本研究分野の基礎知識である原子核や相対論的波動方程式等について学習したことを簡潔にまとめたものを収録した。開発したプログラムのソース
コード リストも付した。

1

原子核の相対論的平均場模型
プログラムの開発と検証

2008 年 2 月
福井大学 大学院 工学研究科 物理工学専攻

18 年度入学 10 番 三和 之浩

目次
序章

3

第1章
1.1
1.2
1.3

模型の物理的位置付け
核力の媒体としての有効中間子場 . . . . . . . . . . . . . . . . . . . . .
中間子場の古典場近似 . . . . . . . . . . . . . . . . . . . . . . . . . . .
負エネルギー状態の核子の無視 . . . . . . . . . . . . . . . . . . . . . .

4
4
5
5

第2章
2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8

模型の定式化
アイソスピン . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
核子の性質 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
考慮するボソンの種類 . . . . . . . . . . . . . . . . . . . . . . . . . . .
本論文で使用する特別な行列のまとめ . . . . . . . . . . . . . . . . . .
自然単位系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
有効場のラグランジアン密度 . . . . . . . . . . . . . . . . . . . . . . .
核子場の従う方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
中間子場の従う方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
8
9
9
10
11
12
12

第3章
3.1
3.2
3.3
3.4
3.5

球対称な系の取り扱い
核子の波動関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
核子の状態を指定するための量子数 . . . . . . . . . . . . . . . . . . . .
軌道の表 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
動径波動関数の従う方程式 . . . . . . . . . . . . . . . . . . . . . . . . .
各種の核子密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13
13
13
14
14
15

第 4 章 核子の波動関数の求め方
4.1 手順の概略 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 r → +0 での漸近形 . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 r → +∞ での漸近形 . . . . . . . . . . . . . . . . . . . . . . . .
4.3 エネルギー固有値の求め方 . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Matching condition . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 エネルギー固有値探索のアルゴ リズム . . . . . . . . . . . . . .
4.3.3 Matching point の選択 . . . . . . . . . . . . . . . . . . . . . . .
4.4 動径グリッド のとり方 . . . . . . . . . . . . . . . . . . . . . . . . . . .

16
16
16
16
17
17
17
18
18
18

第 5 章 中間子場の求め方

20
1

第6章
6.1
6.2
6.3

自己無撞着解の求め方
自己無撞着解を得るためのアルゴ リズム . . . . . . . . . . . . . . . . .
原子核の全エネルギーの表式 . . . . . . . . . . . . . . . . . . . . . . .
Vs , VV の初期値の設定法 . . . . . . . . . . . . . . . . . . . . . . . . . .

22
22
22
22

第 7 章 結果と検証
7.1 208 Pb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.2 スペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 40 Ca . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 密度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.2 スペクトル . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
24
24
25
29
29
30

第 8 章 まとめ

32

参考文献

32

謝辞

35

付 録 A 原子核
A.1 質量と半径 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 核力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 結合エネルギー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.4 魔法数と殻構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36
36
38
40
42

付 録 B 相対論的波動方程式
B.1 相対論的古典力学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 クライン -ゴルドン方程式 . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3 デ ィラック方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.4 中心力ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47
47
49
51
53

付 録 C 全角運動量の固有波動関数

58

付 録 D 4 次の Runge-Kutta 法

62

付 録 E 波動関数の漸近形の導出

63

付 録 F プログラムリスト

67

2

序章
原子核の相対論的模型 [1] は近年、多くの注目を浴びている。それは原子核を 「中間子場を経由して相互作用するデ ィラック方程式に従う核子の系」 として記述する。
原子核の構造や動的性質の非常に顕著な特徴の多くは、すでに単純な模型、特に本論文の主題である相対論的平均場模型の中で記述できるように見える。それは，原子核の強いスピン -軌道力を説明し 、原子核の飽和性を簡単に再現する (文献 [2, 3, 4] 参照)。
模型の相対論的な特徴は、低エネルギーの核構造データという確度の高い情報を用いて高温高密度の原子核物質の性質の推定をすることができるという展望を抱かせる (文献 [4] 参照)。相対論的視点に対する支持は、高エネルギー陽子-原子核散乱のスピン偏極量を説明することにおける、相対論的インパルス近似の成功から来た [14, 15, 16]。
Brueckner-Hartree-Fock 法のレベルでの相対論的な多体論は、核子-核子散乱によって独立に決定できる相互作用を用いて定量的に結合エネルギーや原子核密度の平衡密度を説明することにおいて、最終的に大きな成功を収めた [17, 18, 19]。相対論的な原子核物理についての詳しい解説は文献 [6] で議論されている。相対論的平均場模型で中間子場が古典的 c-数 (classical number) 場として扱われたり、ディラックの海の効果が除外されることについては次章で詳しく説明する．
相対論的平均場模型は、Hartree-Fock 法または時間依存 Hartree-Fock 法の理論的枠組みにおいて使用されることで原子核の構造と低エネルギーの動力学を記述することに大変な成功を収めてきた Skyrme 力 [20] または、Gogny 力 [21] の様な非相対論的な対応物と競合する。そのような力は、原子核媒体中での核子-核子散乱振幅をパラメータ化
する為に仮定されたり、唯一平均場近似においてのみ用いられるべき有効力 (eﬀective
interaction) として考慮されるべきである。８つ前後の少数のパラメータを調節することで、原子核の基底状態の特性、表面振動、巨大共鳴、核分裂、核融合や重イオン反応を定量的に記述することができる。相対論的平均場模型は、核子間に直接力が働くとするのではなく (有効理論としての) 中間子的な自由度を介して相互作用が起きるとするタイプの模型の相対論的な一般化として捉えることもできる。相対論的平均場模型は、いくつかの相対論的効果を付加的に取り込んだ非相対論的模型に較べて、より
自然にかつ同程度に柔軟で強力に現象を再現する力があるということが示されてきた。
本論文の目的は、原子核の相対論的平均場模型の数値解を、密度分布に球対称性を課して求めるための C 言語プログラムを作成し 、その動作の検証を文献 [22] に掲載された種々の計算結果との比較により行ったことの報告である。論文の構成は、最初に模型の物理的意味を説明した後、模型の数学的定式化、数値計算の手順および特に工夫を要した点の説明を行い、動作の検証結果について述べる。付録として開発したプ
ログラムのソースコード の他、原子核の一般的性質の概観、相対論的波動方程式の説
明等についての解説を補足した。

3

第 1 章 模型の物理的位置付け
本章では相対論的平均場模型の構成要素としての中間子場の意味を論じた後、模型の基礎にある２つの近似 —中間子場の古典場近似と負エネルギー状態にある核子の無視— について説明する。

1.1

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

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

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

1.2

中間子場の古典場近似

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

L = ψ̄ (iγ µ ∂µ − mB ) ψ + L(free)meson − gσ Φψ̄ψによって記述される。場の理論の研究においても、このラグランジアン密度で表される模型を使うことがあるが 、その場合は Φ や ψ の量子場としての振る舞いが興味の対象となる。原子核の相対論的平均場模型における意図は、それとは異なっており、２つの基礎となる近似、中間子場の平均場近似と核子のフェルミの海を無視する近似を前提として初めて意味を持つものとしてこの模型のラグランジアンをとらえる。即ち、それらの近似のもとで用いられる有効ラグランジアンとして扱われるのである。
平均場近似は、中間子場の全ての量子的揺動を無視し 、期待値だけに注目するという近似法である。換言すると、全ての中間子場を古典的な c-数 (classical number) の場として扱う近似である。一般に、中間子場（ボソン場）の扱いに関しては、平均場近似という語句と古典場近似という語句は全く同じ意味で用いられており、ニュアンスに一切の差はない。この近似を形式的に表現すれば 、

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

1.3

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

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

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

α<F

∑

free
ϕ̄free
α ϕα

(1.1)

α<F0

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

h: ψ̄ψ :i = ∑α<F0ϕ̄α ϕα −∑free 
+ϕ̄freeα ϕαA∑ϕ̄α ϕαα=1

α<F0として書き直してみる。ここで F0 と F の間にある核子の状態数が A に等しくなるよ
うに F は決定されなければならない (図 1.1 参照)。上式の第１項は真空分極効果を分
離し 、第２項は原子核内の核子の分布を記述する。第１項の量子場効果を無視することは、第２項の α = 1, · · · , A に対応する僅か A 個の占有された粒子状態だけに総和を制限することを意味する。これは no sea 近似 (フェルミの海が空っぽであるときと同じ扱いをする近似という意味) と呼ばれるもので、形式的に表せば 、

h: ϕ̄ϕ :i ∼A∑ϕ̄α ϕαα=1となる。スカラー密度以外のタイプの密度にもこの近似は同様に適用される。
図 1.1 にディラック方程式のスペクトルを模式的に示した。この図において、F0 は、
総核子数 0 に対応するフェルミ準位であり、F は原子核が存在する場合の核子数 A に対応するフェルミ準位である。この図で、もし陽子と中性子の両方が F まで満たされていれば 、この図の表している原子核は 16 O である。なお、相対論的平均場模型ではスカラーポテンシャルとベクトルポテンシャルが 、正エネルギーでは相殺する相対符号を与えられ 、負エネルギーでは強め合う相対符号を与えられる為、負のエネルギーでの核中のポテンシャルの上昇は正エネルギーでの下降より大幅なものとなっている。

6

図 1.1: 動径座標 r に依存した外場中の核子のディラック方程式のスペクトル (模式図)。

7

第 2 章 模型の定式化
2.1

アイソスピン

中性子と陽子という２種類の粒子が 、実は核子という１種類の粒子の２つの異なる状態であると見なすとき、荷電空間 (アイソ空間) という仮想的な空間内の大きさ 21 のスピンに類似した構造をもつ自由度を想定する。このスピンをアイソスピン（ 荷電ス
ピン ）と呼び 、その第３軸方向の成分（ t3 ）の固有状態を中性子および陽子の純粋な状態であると定義する。原子核分野と素粒子分野ではこのラベルと中性子・陽子の対応が正反対であるが、本論文では原子核分野のラベル法に従い下記の通りとする。表 2.1
表 2.1: 粒子とラベル
中性子

t3
τ3

1
2

1

陽子
− 12
−1

で、τ はパウリ行列であり、アイソスピン t と t = 21 τ の関係にある。本論文の定式化
においては、t よりも τ を用いて数式を構成することが多い。

2.2

核子の性質

核子は、内部自由度として{スピン
S = 12
アイソスピン T = 12を持ち、静止質量の実験的に知られている精密な値は、

2

陽子 : mp = 938.27200 MeV/c ,2

中性子 : mn = 939.56533 MeV/cである [25]。但し 、プログラムの検証のための計算を行うときは、比較対象の文献 [22]に合わせて、mp = mn = 939 MeV/c2 とする。

8

2.3

考慮するボソンの種類

模型に取り入れる (核子間の相互作用の媒介となる) ボソン場を、表 2.2 に列挙した。
中間子以外に、電磁相互作用を媒介する光子 (γ) も含まれている。表 2.2 において、Sはスピン、T はアイソスピン、m は静止質量 (MeV/c2 ) 、g 2 は後述の核子との結合定数の２乗である。光子の結合定数は gγ = e であり、２乗して h̄c で割り無次元化すると微細構造定数 α = 1/137.036 の 4π 倍に等しくなる。
表 2.2: 中間子の種類

σ
ω
ρ
γ

2.4

S
0
1
1
1

g2
(無次元数)h̄c

M (MeV/c2 )
520
783
770
0

T
0
0
1
0

109.626
190.431
16.30654π
137.036

本論文で使用する特別な行列のまとめ

本論文の数式で使用する特別な行列を以下にまとめておく。

σ

=

(
(

σ1

=

α

=

(
(

αi

=
[

β

=
(

γ0

=

)

σ1 , σ2 , σ3
0 1
1 0

)

(

0 −i
i 0

, σ2 =))(, σ3 =

1 0
0 −1

)

α1 , α2 , α3
0 σi
σi 0
I 0
0 −I
I 0
0 −I

)

(i = 1, 2, 3)
]

(

, I =
)

1 0
0 1

(

, γi =

)

0 σi
−σi 0

)

なお、γ 行列の上記の表現は標準表現 (Standard representation) 、あるいは Dirac 表現、あるいは Pauli-Dirac 表現と呼ばれているが、他にカイラル表現 (Chiral representation,
Wyle 表現とも言う) や Majorana 表現という上記とは異なった表現が使われることもあるので、数式を他の文献と比較する際には、個々の文献で採用しているガンマ行列の表現方法を確かめることが必須である。

9

2.5

自然単位系

以下に記述する数式は、自然単位系で表されている。自然単位系ではh̄ = c = 1(無次元数としての１）となる。
自然単位系では、唯一つの単位で全ての物理量が測られる。原子核物理なら、その単位としてエネルギーの単位 MeV が便利である。この場合、例えば 、質量 m に数値を代入する必要があれば 、m を kg 単位で表したときの数値 (陽子なら mp = 1.672 × 1027 kg) でなく、静止質量エネルギー mc2 を MeV 単位で表した数値 (陽子なら mp = 938.271MeV)を代入するのである。結合定数の２乗 g 2 は自然単位系では無次元なので、表 2.2 の値 (ρ中間子なら gρ2 = 109.626) を使えばよい。また、密度 ρ の数値は m−3 単位で表したとき3の値 (原子核の飽和密度なら 1.6×10−46 m−3 ) を使うのではなく、ρ (h̄c) を MeV3 単位で3表したときの数値 (原子核の飽和密度なら 0.16fm−3 · (197MeV · fm) = 1.2 × 106 MeV3 )を使うことになる。
しかし 、密度を MeV3 単位で表すことは一般的ではなく、わかりにくい。やはり密度は fm−3 単位で表したいものである。そのためには、g 2 の値として MeV·fm 単位で表した数値を用いることにすればよい。自然単位系では、

h̄c = 1かつh̄c = 197MeV · fm
なので197MeV · fm = 1
である。したがって、g 2 の数値として σ 中間子の場合は、

g 2 = 109.626 = 109.626 × 197.327MeV · fm
= 21632.2MeV · fm
を用いればよい。このように、原子核物理学では数式は自然単位系で表しても、数値としての扱いでは MeV と fm の２つの単位を混在させるのが便利である。
次節で導入するラグランジアン密度を通常の単位系で表すには、LN については∂µ → h̄c∂µ ,m → mc2と置き換え、Lσ , Lω , Lρ については

m→

mc
h̄
10

と置き換えればよい。他の項は変える必要はない。
ちなみに物理量 P の単位の次元を [P] で表すことにし 、E, L をそれぞれエネルギーと長さの単位とすれば
[

ψ̄ψ

]

=

[φ] =
[g] =

[
[
[

]

L−3 ,
E 2 L− 2
1

1
2

E L

1

1
2

]

]

,

である。ψ と g の次元は、クーロン力の場合の関係式
[

にならって

]

e2 1
= [eA0 ] = [E]4πε0 r[g21rと定義すれば上記のように定まる。

2.6]= [gφ] = [E]

有効場のラグランジアン密度

原子核の相対論的平均場模型で用いる有効場のラグランジアン密度は下記の通りである。

L = LN + Lσ + Lω + Lρ + Lγ + Lc.
核子部分は、

LN = ψ̄ (iγµ ∂ µ − m) ψ
で与えられる。σ 中間子部分は、}
1{(∇φσ )2 + m2σ φ2σ2で与えられる。ω 中間子部分は φω を ω 中間子場の時間成分 (第０成分) として
}
1{
Lω =
(∇φω )2 + m2ω φ2ω
2で与えられる。ρ 中間子部分は φρ を ρ 中間子場の時間成分で、アイソスピンの３軸成
分が０の成分として}
1{Lρ =(∇φρ )2 + m2ρ φ2ρ
2で与えられる。光子場部分は A0 を光子場の時間成分 (即ち、電位/e) として
1
Lγ = (∇A0 )22で与えられる。核子場と中間子場・光子場の結合 (Coupling) は
1 − τ3
Lc = −gσ ψ̄φσ ψ − gω ψ̄γ0 φω ψ − gρ ψ̄γ0 τ3 φρ ψ − eψ̄
γ0 A0 ψ2
で与えられる。

Lσ = −

11

2.7

核子場の従う方程式

核子は下記の Dirac 方程式に従う。

{−iα · ∇ + β (m + Vs) + VV} ψi = Ei ψi .
ここで、Vs はスカラーポテンシャルであり、下式で与えられる。

Vs = −gσ φσ .
VV はベクトルポテンシャルの時間成分であり、下式で与えられる。
VV = gω φω + gρ τ3 φρ + e

1 − τ3 0A.
2

Dirac 方程式については付録 B を参照せよ。

2.8

中間子場の従う方程式

古典的中間子場は遮蔽された（ screened ） Poisson 方程式に従う。また、質量がゼロである光子の場は狭義の（ proper ）Poisson 方程式に従う。各々の場の源 (source term)は後述する各種の核子密度に結合定数を掛けたもので与えられる。
(()∆ − m2σ φσ = −gσ ρs,)∆ − m2ω φω = −gω ρV,()∆ − m2ρ φρ = −gρ ρ3 ,∆A0 = −eρp.

12

第 3 章 球対称な系の取り扱い
3.1

核子の波動関数球対称系では核子の波動関数は下記のような形を持つ。


G

(r)



3
i $jntΦljm (θ, ϕ, ms) rψ$jmnt3 (r, θ, ϕ, ms, mt) =δt3 mt ,F$jnt3 (r)−Φl0 jm (θ, ϕ, ms)
r

但し 、

1
l = j + $,21
l0 = j − $,2∑1h l, m1 , , m2 | jm iYl,m1 (θ, ϕ) δmsm2
2
m1 ,m2
1= h l, m − ms, , ms | jm iYl,m−ms (θ, ϕ)21である。ここで、< l, m1 , 2 , m2 |j, m > はクレプシュ・ゴーダン係数、Yl,m (θ, ϕ) は球面調和関数である。詳しくは付録 C を参照せよ。
Φljm (θ, ψ, ms) =3.2

核子の状態を指定するための量子数球対称系では、核子の状態は、j, $ でラベルされる。これらのとりうる値は下記のとおりである。j は全角運動量であり、$ は Dirac 方程式に特有の量子数である。
{

j = 12 , 32 , 52 , · · ·$ = ±1.

(半奇数).

この他に、エネルギー固有値を指定するには動径波動関数の主要成分 G のノード 数 n(0 以上の整数値をとる）などが必要である。$ = −1 のとき、小さい成分 F のノード数は n に等しい。$ = 1 のときは、F のノード 数は n + 1 である。
さらに、方程式を簡明に表すための副次的な量として下記のものを使用する。
()

1
κ=$ j+= ±1, ±2, ±3, · · · ,2
l = j + 21 $ (G の軌道角運動量)l0 = j − 21 $ (F の軌道角運動量)13}= 0, 1, 2, · · · .

3.3

軌道の表

球対称系の Dirac 方程式の状態を指定する量子数と、よく知られた Schrödinger 方程式の軌道を表す記号との対応表を以下に示しておく。
表 3.1: 軌道の表

$
−1
1

l l0
0 1
1 0

κ 軌道を表す記号
−1
s1
2
1
p1

3
2

−1
1

1
2

2
1

−2
2

p3
2
d3

5
2

−1
1

2
3

3
2

−3
3

d5
2
f5

7
2

−1
1

3
4

4
3

−4
4

f7
2
g7

9
2

−1
1

4
5

5
4

−5
5

g9
2
h9

11
2

−1
1

5
6

6
5

−6
6

h 11
2
i 11

13
2

−1
1

6
7

7
6

−7
7

i 13
2
j 13

j
1
2

3.4

2

2

2

2

2

2

2

動径波動関数の従う方程式

核子の動径波動関数は下記の方程式に従う。

d
dr

(

G
F

)

(

=

− κr , µ + ε
κµ − ε,
r

)(

G
F

)

但し 、

µ = m + V s, ε = E − V Vとした。
なお、large component に G 、small component に F という文字を使用するのは文献[5] 、[6] などであり、本論文もこれに従う。一方、large component に F 、small componentに G という文字を使用する文献もあるので、注意が必要である。例えば 、文献 [9] 、[10] 、[11] がそうである。

14

上記の方程式のもつ対称性のひとつとして、κ および ² の符号を逆転させると、F とG が入れ替わることが挙げられる。なお、κ の符号は $ の符号をそのまま受け継いだものである。

3.5

各種の核子密度i を軌道の番号、ni を i 番目の軌道に入っている核子の個数 (0 ≤ ni ≤ 2ji+1 ) として、
下記のような各種の核子密度を定義する。
}
1 ∑ {
2
2
n
|G
(r)
|
−
|F
(r)
|
i
i
i
4πr2 i
}
1 ∑ {
2
2
n
|G
(r)
|
+
|F
(r)
|
ρV (r) =iii
4πr2 i{}1 ∑2
2
ρ3 (r) =
n
(τ
)
|G
(r)
|
+
|F
(r)
|
i
3 i
i
i
4πr2 i()}
1 ∑1 − τ3 {22
ρp (r) =
|G
(r)
|
−
|F
(r)
|
n
i
i
i
4πr2 i2iρs (r) =なお、i 番目の軌道が中性子の軌道なら((τ3 )i = 1,

1 − τ3
2

)

(

)

=0
i

陽子の軌道なら(τ3 )i = −1,

1 − τ32=1iである。ρs はスカラー密度、ρV はバリオン密度 (核子流の時間成分、即ち、通常の密度) 、ρ3 はアイソベクトル密度のアイソ空間の第３軸成分、ρp は陽子密度である。

15

第 4 章 核子の波動関数の求め方
4.1

手順の概略

境界条件を与えて Dirac 方程式の動径座標部分を４次の Runge-Kutta 法 (付録 D 参照)
で解けば 、動径波動関数を求めることができる。エネルギー E の値を仮定して、r → +0での正則解の漸近形から出発して r の増加する方向へ解き進めた解と、r → +∞ での正則解の漸近形から出発して r の減少する方向へ解き進めた解とを適当な点で接続し 、波動関数をリスケールすることで G 、F の２成分がともに連続に接続できるとき、仮定した E はエネルギー固有値となっているのである。連続に接続できなかった場合に、次にどのような値に E を変更すればよいかを指示するアルゴ リズムは後で説明する。

4.2

境界条件

4.2.1

r → +0 での漸近形

r = 0 でのテーラー展開は下記のとおりである。
()
ε2 − µ2 l+3
r + O rl+5
2 (2l + 3)()
l+1+κ l l+3+κ
F =r −(ε − µ) rl+2 + O rl+4
ε+µ2 (2l + 3)

G = rl+1 −但し 、{

ε = E − VV (r = 0)µ = m + Vs (r = 0)とした。導出方法は付録 E に詳細に記述した。
注意を払うべき点は、F の第１項の係数が $ = −1 のとき 0 になることである。
(

1
1
...
l+1+κ = j+ $+1+$ j+22= j + 1 + $ (j + 1)
= (j + 1) ($ + 1)

16)

しかし 、第２項の係数は、$ = −1 のとき、l + 3 + κ = (l + 1 + κ) + 2 = 2 なので 0 ではない。よって、F の展開は少なくとも先頭の２項を使用すべきである。
また、G の展開も２項を使用するべきである。そうしなければ 、s 波 (l = 0) の密度が接続点で不連続な１階微分を持つことになる。
なお、補足すると、VV, Vs を r2 の巾級数として展開した場合の r2 の項を考慮すると、前記の G,F の展開の第２項が変更されることになる。現在は、r2 の項の係数が小さいと期待して定数項しか考慮していないのであるが 、理想的には全ての項について展開の次数を揃えることが望ましい。将来この研究を発展させる人への提言として、VV,Vsの r2 の項の係数を含めて G と F のテーラー展開の係数の決定をやり直すという課題が残っていることを、ここに付記しておく。

r → +∞ での漸近形

4.2.2

r → +∞ での漸近形は下記の通りである。
√

G = e− m −E r ,
√

F = −

2

2

m − E −√m2 −E 2 re.
m+E

導出方法はやはり付録 E に記述してある。

4.3

エネルギー固有値の求め方

4.3.1

Matching condition

G1 (r) , F1 (r) を r → +0 で正しい漸近形を持つ解、G2 (r) , F2 (r) を r → +∞ で正しい漸近形を持つ解とする。
２つの解が一致する条件は、任意に取り得る matching point r = rm にて G，F 共に連続に接続できることである。G，F は連立１階常微分方程式の解なので，関数値が一致すれば (連続に接続できれば ) 任意の階数の微分も一致するからである。従って接続条件は
G2 (rm)G1 (rm)=F1 (rm)F2 (rm)である。

D = G1 (rm) F2 (rm) − F1 (rm) G2 (rm)
とすれば

D=0
が接続条件である。

17

4.3.2エネルギー固有値探索のアルゴリズム

D 6= 0 のとき、固有値 E0 は E より大きい側にあるのか小さい側にあるのかは D の正負からは判断できない。しかし
D0 = G1 (rm) G2 (rm) {G1 (rm) F2 (rm) − F1 (rm) G2 (rm)}と定義される D0 を用いると以下に述べる通り判断できる (証明は省略する)。
エネルギー E が与えられたときの G のノード 数 n(0 < r < rm に G1 (r) の持つノード 数と rm < r に G2 (r) の持つノード 数の和) と求めたい固有状態のノード 数 n0 とを比較して、求めるべき固有値 E0 の存在範囲について以下のことが言える。

n < n0なら n > n0なら

E < E0E > E00 D < 0 なら E < E0n = n0なら  D0 = 0 なら E = E0





 D 0 > 0 なら E > E
0

4.3.3

Matching point の選択

Schrödinger 方程式に還元して言えば 、matching point で E <potential energy(遠心力ポテンシャルを含む) であると、数値解法としては失敗することがある。これは、ポテンシャルバリア内においては激しい指数関数的発散が生じる為、倍精度実数の精度一杯に E の値を調節しても、図 4.1 の上段に模式的に示したように、G1 ,G2 共に r = rm
において最大の大きさを持つような状況しか得られないことがあり得るからである。一方、図 4.1 の下段に示したように、E >potential energy(遠心力ポテンシャルを含む)である場所に matching point を設定すると、このような問題は起きない。
なお、過去においては実数の桁あふれへの対処が面倒な問題としてよく挙げられていたが 、現在の標準的な規格では倍精度実数は 10300 以上の大きさの数でも扱えるので、本研究では、桁あふれは特に問題とはならなかった。

4.4

動径グリッド のとり方

Dirac 方程式を 4 次の Runge-Kutta 法で解いたり、Poisson 方程式を高次の積分公式
で求めたりする際には、波動関数や平均場を動径座標上のグリッド 点 (標本点) 上の値で表現する。本研究では、以下のように、等間隔グリッド を採用した。
ri = h · i,

0 ≤ i ≤ 3000,

h = 0.01 fm .
非相対論的計算では原子核の全エネルギーを 100 keV 程度の精度で計算するには 、
h = 0.2 fm 程度で十分と言われている。しかし 、本研究では Dirac 方程式の精度がまだ余り議論されていないことから、安全の為、１桁小さい値を使用した．

18

E

G1

G2
potential

0

rm

r

E
G1

G2

0

r
rm

図 4.1: Matching point の選択 (上：失敗例 下：成功例)

なお、将来的には等間隔のグリッドでなく、r = 0 付近で等比数列に近付くようなグリッドを採用することが望ましい。その理由は、遠心力ポテンシャルが r → +0 で発散することに 、数値計算的に正し く対処することが望ましいからである。等間隔グリッド では r → +0 の漸近形にどこで接続するかが微妙な調整を要する課題として残っているのだが 、等比数列グリッド の採用でその調整は全く不要となる。等比数列グリッドについての詳細は文献 [7] を参照せよ。

19

第 5 章 中間子場の求め方
Screened Poisson 方程式()∆ − m2 φ = −gρ

に、球対称性を課すと次のような常微分方程式が得られる
(

)

2 d
d2
+
− m2 φ (r) = −gρ (r) .
2
dr
r dr

この方程式の解は以下のとおりである

φ (r) = ge−mr ∫ rsinh (mr) ∫ ∞ −mr0
sinh (mr0 ) ρ (r0 ) r0 dr0 + geρ (r0 ) r0 dr0 .
mr 0
mr
r

クーロンポテンシャルの場合は m → 0 とすると得られる
∫ ∞
g∫r
0
02
0
φ (r) =
ρ (r ) r dr + g
ρ (r0 ) r0 dr0 .
r 0rρ (r) は ρ (−r) = ρ (r) として r < 0 側へ滑らかに解析接続される。この性質を利用して r = 0 付近を含めて高精度の数値積分を行うことができる。本論文では、隣接する６個のデータ点に r の５次多項式を当てはめ、３個目と４個目のデータ点の間の積分を、この５次多項式の積分値で近似した。詳しくは文献 [8] を参照せよ。
また、ρ (r) は r → ∞ で速やかに 0 に収束するので、積分範囲は有限の r まででよい。本論文では r = 30 fm まで積分している。
プログラムチェックのために行った数多くの数値計算から一例を図を示す。入力として与えた電荷密度 ρ は
ρ= (

)

1

(

1 + exp r−R
× 1 + exp −r−Raa][
) fm−3とした。但し 、パラメータ R, a はそれぞれR = 7.0fm, a = 0.7fmである。
図 5.1 はこの電荷密度であり、図 5.2 は６点公式で求めたクーロンポテンシャルである。参考までに台形公式を用いて求めたものも図 5.3 に示した。

20

図 5.1: 電荷密度

図 5.2: クーロンポテンシャル (６点公式)
500

Coulomb potential

450
400

V (MeV)

350
300
250
200
150
100
50
0
0

10

20

30

40

r (fm)

図 5.3: クーロンポテンシャル (台形公式)

2150

第 6 章 自己無撞着解の求め方
自己無撞着解を得るためのアルゴリズム

6.1

自己無撞着解を得るための手順は、まず、Vs 、VV に初期値を仮定した上で以下の１〜４の手順を繰り返すことである。

1. Vs 、VV に対する核子のエネルギー固有値 Ei と波動関数を求める。
2. Ei の低い状態から占拠核子数 ni を割り当て、各種の核子密度を得る。
3. 核子密度を源とする中間子場の古典解を求め、得られた中間子場を組み合わせて
Vs 、VV を作る。
4. ３. で得られた Vs 、VV が１. で使用した Vs 、VV と異なるなら、１. に戻る。逆に差が非常に小さければ 、自己無撞着解が得られたことになる。

原子核の全エネルギーの表式

6.2

自己無撞着解が得られたら、原子核の全エネルギー E は以下の表式で計算すること
ができる

E =

∑
i

ni Ei −

1∫ ∞{−gσ φσ (r) ρs (r) + gω φω (r) ρv (r) + gρ φρ (r) ρ3 (r)2 0}+eA0 (r) ρp (r) 4πr2 dr.

6.3

V , V の初期値の設定法
sV

空間的に一様な系では ∆φ = 0 なので()∆ − m2 φ = −gρ

より、

φ=

g
ρ
m2

22

を得る。

ρn = 0.16

N ( −3 )Z ( −3 )fm, ρp = 0.16
fm
AAと仮定すればρV = ρn + ρp = 0.16,. ρ × 0.9,
ρs =
. V
N −Z
ρ3 = ρn − ρp = 0.16
Aが導かれ 、ポテンシャルの深さはVV(0) =

gρ2
gω2
ρ
+
τ
3 2 ρ3 ,m2ω VmρVs(0) = −

gσ2
ρ,
m2σ s

となる。そこでVV (r) =Vs (r) =VV(0)(1 + exp − r−R
a
Vs(0)
(

1 + exp − r−Ra) +Vcoul (r))とすればよい。
但し 、

a = 0.7 fm,
1

R = r0 A 3 ,r0 = 1.2 fm,とし 、クーロンポテンシャルは陽子に対してのみ働き、
{( )2 }
 Z e2 3 − 1 r

Vcoul (r) =  4πR2
 Ze

2

2

R

(r ≤ R) ,(r > R) ,4πrと近似するとよい。

23

第 7 章 結果と検証
7.1

208

7.1.1

密度

Pb

作成したプログラムを用いて 208 Pb の中性子および陽子のバリオン密度および スカ
ラー密度を計算し 、図 7.1 に示した。模型のパラメータは文献 [22] と同じである。

図 7.1: 208 Pb の密度

原子核内部のバリオン密度は、中性子は赤の実線、陽子は緑の実線で示してある。ス
カラー密度は、中性子は青の実線、陽子は紫の実線で示してある。中性子のバリオン
密度は、最も荒い見積り値

0.16 ×

N
A

126
208
= 0.097fm−3= 0.16 ×24

の 90 ％程度となっている。また、陽子のバリオン密度は0.16 ×Z
82= 0.16 ×A
208= 0.063fm−3とほぼ一致している。また、中性子・陽子共にスカラー密度はバリオン密度の約 90 ％である。
図 7.1 を見ると、neutron は中心から 1.5 fm にかけて密度が濃くなり、proton は薄くなっていることがわかる。これはフェルミ準位間近の軌道の特徴が現れているのである。その後は両者とも、5 fm もしくは 6 fm まで安定しており、その後 9 fm 辺りでは
1
中心部分の 100以下にまで減少している。

7.1.2

スペクト ル

表 7.1 は 208 Pb の中性子の一粒子準位について、本研究での計算値と文献 [22](表中でSerot と表記) の値を比較したものである。表 7.2 は陽子の一粒子準位についての比較である。
また、図 7.2 は表 7.1 と表 7.2 をプロットしたもので、× は今回作成したプログラム
による計算結果で、実線は文献 [22] の結果を繋いだもので、下が中性子で上が陽子を表している。また、横軸はエネルギー準位の番号 (低いものから 1, 2, 3, · · · と番号を割り振った) である。
Serot のスペクトルは文献 [22] の Fig.13 から読み取ったもので、まず Fig13 をスキャンでコンピュータに読み込み、始めに y 軸の最小と最大のプロット点を、それからそれぞれの軌道のプロット点を読み取りプログラムで変換した。準位を表す水平な線分の幅は 0.2MeV であり、読み取りの精度はその 21 倍程度として 0.1MeV である。また、
Fig.13 は手書きではなくプロッタによる作図なので、作図の精度は読み取りの精度よりはるかに高いと考えられる。
図 7.2 を見ると、neuton,proton 共にうまく再現できている。最も大きくずれているのは neutron の下から３番目の軌道 (1p 1 ) で、ずれは 340keV である。この大きさのず
2
れはグラフの読み取り誤差より明らかに大きく、ど ちらかの計算の不備であると考え
られる。

25

表 7.1: スペクトル表 (Neutron)

orbital Serot(MeV) This work(MeV) diﬀerence(MeV)
1s 1-64.983-64.808-0.1752
1p 3-58.139-57.979
-0.16
2
1p 1
-57.543-57.203
-0.34
2
1d 5
-49.679-49.691
0.012
2
1d 3
-47.766-47.781
0.015
2
2s 1
-44.025-44.064
0.039
2
1f 7
-40.411-40.409-0.0022
1f 5-36.883-36.854-0.0292
2p 3-32.419-32.244-0.1752
2p 1-30.846-30.822-0.0242
1g 9-30.421-30.535
0.114
2
1g 7
-24.895-24.979
0.085
2
2d 5
-20.941-20.890-0.0512
1h 11-20.431-20.397-0.0342
2d 3-18.731-18.672-0.0592
3s 1-17.838-17.883
0.045
2
1h 9
-12.779-12.794
0.015
2
2f 7
-10.356-10.361
0.005
2
1i 13
-10.101-10.262
0.161
2
2f 5
-7.635-7.756
0.121
2
3p 3
-7.21
-7.2960.086
2
3p 1-6.402-6.375-0.02722g 9-1.2182
1i 11-1.1032
1j 15-0.3692
3s 1-0.08723d 5-0.0072

26

表 7.2: スペクトル表 (Proton)

orbital Serot(MeV) This work(MeV) diﬀerence(MeV)
1s 1-52.110-52.001-0.1092
1p 3-45.988-45.873-0.1152
1p 1-45.011-44.961-0.0502
1d 5-38.209-38.117-0.0922
1d 3-36.083-35.991-0.0922
2s 1-31.45-31.460.01
2
1f 7-29.239-29.277
0.038
2
1f 5
-25.456-25.457
0.001
2
2p 3
-20.184-20.156-0.0282
1g 9-19.801-19.762-0.0392
2p 1-18.824-18.736-0.0882
1g 7-13.935-13.941
0.006
2
1h 11
-9.939-9.893-0.04622d 5-9.174-9.1970.023
22d 3-6.963-7.0190.056
2
3s 3-5.900-5.857-0.0432
1h 9-2.063227

図 7.2: 208 Pb の Serot との比較 (下：neutron 上：proton)

28

7.2

40

7.2.1

密度

Ca

40

Ca についても 208 Pb のときと同様に自作のプログラムで計算し 、図 7.3 に示した。
また、文献 [22] の Fig.8 から読み取った ρp(HARTREE) との同じ図中に描き込んだ。読
み取った方法は 208 Pb のスペクトルのときと同様で、今度は x 軸についてもプロット点を読み取った。
中性子は赤線 (バリオン密度) と青線 (スカラー密度) 、陽子は緑線 (バリオン密度) と紫線 (スカラー密度) で Serot の陽子密度は橙線で示してある。これらは荒い見積り値の0.16 ×N
Z
20= 0.16 ×= 0.16 ×A
A40= 0.08fm−3とほぼ一致している。また、ρp (HARTREE) と陽子のバリオン密度は完全に重なっており、本論文と文献 [22] で完全な一致を見た。

図 7.3: 40 Ca の密度

Ca の密度については、中心付近では proton の |G|2 + |F |2 と neutron の |G|2 − |F |2の濃さは同程度であり、どの密度も 1.5 fm 辺りまでは減少し 2 fm 前後では安定してい1以下にまで減少していることがわかる。
るがその後、6 fm 辺りでは中心部分の 100
40

29

7.2.2スペクト ル

40

Ca のスペクトルも 208 Pb と同様に、Serot らの文献 [22] の Fig.12 から読み取ったス
ペクトルと比較してみた。表 7.3 は陽子の各軌道の準位を、文献 [22] と本研究とで比較したものである。また、図 7.4 は表 7.3 をプロットしたもので、× や実線、横軸についても 208 Pb の図 7.2 と同じである。
表 7.3: スペクトル表 (proton)

orbital Serot(MeV) This work(MeV) diﬀerence
1s 1-47.11-46.674-0.4362
1p 3-30.868-30.703-0.1652
1p 1-25.355-25.230-0.1252
1d 5-15.322-15.088-0.2342
2s 1-7.623-7.321-0.3022
1d 3-7.027-6.708-0.31922s 1-0.6432

図 7.4 を見ると 40 Ca では、一番下と下から５、６番目の軌道 (1s 1 , 2s 1 1d 3 ) でずれが222大きい。このずれもまた、ど ちらかの計算の不備を示唆していると思われる。

30

図 7.4: 40 Ca の Serot との比較 (proton)31

第 8 章 まとめ
本論文では、原子核の相対論的平均場模型の解を、密度分布に球対称性を課して求める数値計算プログラムを開発し 、文献 [22] に掲載された種々の計算結果と比較してその正しさを検証した。
第一章では、原子核の相対論的平均場模型における有効場の概念を説明し 、量子場の理論から平均場模型を導くために必要な２つの近似について論じた。第一の近似は核子間の相互作用の媒介として導入された各種の中間子の場のもつ量子力学的揺動を無視してその平均値だけに注目するという、中間子場の古典場近似である。第二の近似は 、核子がスピン 12 の粒子の従う相対論的な波動方程式である Dirac 方程式に従うとした場合に現れる負のエネルギーの解を無視するという所謂 no-sea 近似である。
第２章では、原子核の相対論的平均場模型の簡潔にしてかつ完全な定式化を与えた。
続く第３章では、密度分布に球対称性を課す場合の取扱いについて、状態を指定するための量子数のセットから動径波動関数の従う微分方程式に至るまで、簡潔なまとめを示した。
第４章および第５章では、数値解法上の重要ポイントとして、r = 0 および r → ∞での解析解の形や、両側から解析接続して得られた２つの解が連続に接続できるという条件を満たすようにエネルギー固有値を探索する効率的な手順などについて説明した。６章では、核子の状態と中間子場との間に無撞着性が成立する自己無撞着解を求めるアルゴ リズムを説明した。
最後の第７章においては開発したプログラムの結果を文献に記載された種々の計算結果と比較して違いが十分に小さいことから、プログラムの検証に成功したと結論づけた。
さらに付録として、本研究分野の基礎知識として研究を始める前に学習に時間を費した、原子核や相対論的波動方程式等についてのまとめを収録した。開発したプログ
ラムのソースコード リストも付した。

32

参考文献
[1] P.-G.Reinhard, “The relativistic mean-ﬁeld description of nuclei and the nuclear
dynamics” Report on Progress in Physics 52, 439 (1989)[2] H.P.Duerr, “Relativistic Eﬀects in Nuclear Forces”, Phys. Rev. 103, 469 (1956)[3] L.D.Miller and A.E.S.Green, “Relativistic Self-Consistent Meson Field Theory ofSpherical Nuclei”, Phys. Rev. C 5, 241 (1972)[4] J.D.Walecka, “-A theory of highly condensed matter-” Annals of Physics, Volume
83, Issue 2, 491 (1974)[5] W.Greiner, “Relativistic Quantum Mechanics - Wave Equations -”, Springer Verlag, Berlin (1990)[6] B.D.Serot and J.D.Walecka, “The Relativistic Nuclear Many-Body Problem”, Advances in Nuclear Physics, vol.16, Plenum Press, New York (1986)[7] 秋山大樹, 上川暢介, 「 Dirac-Fock 法による原子の基底状態の電子配置」 福井大学
工学部 物理工学科 卒業論文, 2008 年 2 月.
[8] 山田昌平 福井大学工学部 物理工学科 卒業論文 「原子核の平均場模型」, 2007 年 2
月.
[9] A. メシア, 量子力学 vol.3 小出昭一郎 他訳, 東京図書 (1972) p.180 - 187, 217 - 220.
[10] 西島和彦, 「相対論的量子力学」 (1973)[11] C.Itzykson and J.-B.Zuber, “Quantum Field Theory”, McGraw-Hill (1985)[12] 有馬朗人, 原子と原子核 (朝倉書店, 1982) p.105, 106, 147 - 151, 157 - 159, 162,
163.
[13] A. Bohr and B. R. Mottelson,Nuclear Structure vol.1 (W. A. Benjamin , Inc. ,1972) p.224, 325.
[14] B.C.Clark et al, Phys. Rev. Lett. 50, 1644 (1983).
[15] B.C.Clark et al, Phys. Rev. C 28, 1421 (1983).
[16] J.R.Shepard et al, Phys. Rev. Lett. 50, 1443 (1983).
33

[17] R.Brockmann and R.Machleidt, Phys. Lett. 149B, 283 (1984).
[18] B.TerHaar and R.Malﬂiet, Phys. Lett. 172B, 10 (1986).
[19] B.TerHaar and R,Malﬂiet, Phys. Rev. Lett. 56, 1237 (1986).
[20] P.Quentin and H.Flocard, Ann. Rev. Nucl. Part. Sci. 28, 523 (1978).
[21] J.Decharge and D.Gogny, Phys. Rev. C 21, 1568 (1980).
[22] C.J. Horowitz and B.D. Serot, Nuclear Physics, A368, 503-528 (1981).
[23] K.Holinde, Phys. Rep. 68, 122 (1981).
[24] R.Machleidt et al. Phys. Rep. 149, 1 (1987).
[25] Particle Physics Booklet, Particle Data Group, LBNL and CERN (2002).
[26] H. L. Anderson et al. Phys. Rev. 187, 1565 (1964).

34

謝辞
本論文を書くにあたり、終始変わらぬ懇切丁寧な指導と助言いただいた鈴木敏男
先生、林明久先生、田嶋直樹先生には感謝の念に耐えません。鈴木敏男先生には学部４年生の卒業研究配属以来、修士課程２年生の６月頃までの指導教員として、大学院研究科長としての激務の間をぬって御指導いただきました。その後、修士論文提出までは田嶋先生が指導教員として、論文完成まで担当していただきました。稚拙な意見にも耳を傾けていただき、的確な助言で研究をより良い方向へと導きいただいたことを、心より感謝申し上げます。また物理工学科の諸先生方に対しましても、日頃の厚い御指導に深く感謝しています。最後に１年間共に学び 、励まし合ってきた原子核理論教室の大学院生、卒業研究配属生諸氏にも本当に感謝し 、謝辞の言葉とさせていただきます。

2008 年 2 月
三和 之浩

35

付 録A

原子核

本付録の執筆に際しては参考文献 [12] の第 11 章を参考にした。

A.1

質量と半径

同じ原子番号 Z をもつ元素の集まりを同位体 (アイソトープ ) という。同位体は化学的にはほとんど 同じ性質であり、物理的手段で各元素を同位体に分けて、その原子量を測るとほとんど 整数である。この整数 A を原子核の質量数という。12 C の原子量の1/12 を単位として原子量を測ることがあり、この 1 単位は

1 amu = 1.66057 × 10−24 g
1 amu c2 = 931.5016 MeV
で、この単位を原子質量単位という。
水素の陽子や、そのほぼ 4 倍の質量をもつヘリウムの原子核を他の原子核に衝突させることによって、原子核の大きさを知ることができる。このヘリウムの原子核を α粒子とよび 、また標的にした原子核のことを標的核という。

図 A.1: ラザフォード 散乱

一般に、入射粒子と標的核の電荷をそれぞれ ze, Ze とすると、この間にクーロン力が働く。クーロン力による散乱は、ラザフォード 散乱 [12] とよばれ 、強度 I0 、速度 v 、
電荷 ze 、質量 m なる粒子の散乱後の流れは、散乱角 θ, ϕ 方向への単位立方角あたり

I (θ, ϕ) = I0z 2 Z 2 e4
14
2
44m v sin (θ/2)
36

(A.1)

の強度をもつことが予想される (図 A.1)。入射粒子がもっとも原子核に近づく距離 rmin
は
{

1
zZe2
1
+
rmin =
mv 2
sin (θ/2)

}

(A.2)

で、210 Po の α 崩壊によって生じる α 粒子のエネルギーは 5.4 MeV である。この α 粒
子を Au の原子に衝突させた結果、(A.1) 式がよく成り立っていた。rmin は θ = π のときもっとも小さくなる。

mv 2= 5.4 MeV
2
Z (Au) = 79e21=
(微細構造定数)h̄c137h̄c = 197 × 10−13 MeV cm
を用いると4 × 79e2
mv 2
4 × 79 × 197 × 10−13=cm
137 × 2 × 5.4= 43 fmrmin =となる。但し 、1 fm= 10−13 cm である。即ち Po から自然放射能で放出される α 線を用いた実験では、Au の原子核の半径は 43 fm より小さいという結論が得られる [12]。
ところが、加速して十分エネルギーを上げた α 粒子を原子核に衝突させると、ある角度以上で散乱強度がラザフォード の (A.1) 式からずれてくる。図 A.2 は 208 Pb を標的にし 、104 MeV のエネルギーをもった α 粒子を衝突させたときの散乱強度と、(A.1) 式の値の比を、散乱角 θ の関数として描いてある。散乱強度は、θ ≤ 10 °でほぼラザフォー
ド 散乱と一致するが、θ > 10 °ではそれよりずっと小さくなる。これは、θ ≤ 10 °の角度で散乱されるとき、α 粒子は原子核に接触しないが 、θ > 10 °では接触するため、α核子と標的核間の相互作用がもはや純粋なクーロン力ではなくなるからであると解釈することができる。そこで、θ = 10 °のときの rmin の値から原子核の半径と α 粒子の半径の和が求まる。
Eα = mv 2 /2 = 104 MeV, θ = 10 °を (A.2) 式に代入すると

rmin =

2 × 82 × 197
× 12.5 = 14.1 fm
137 × 2 × 104である。かくして原子核 208 Pb の半径はほぼ 1 fm の 10 倍ぐらいであることがわかった。
α 粒子や陽子を種々の原子核に衝突させ、上述のような分析や、さらに量子力学による定量的な研究を行うことによって、質量数 A なる原子核の半径 R は

R = r0 A1/3

(A.3)

r0 = 1.4 ∼ 1.5 fm

(A.4)

37

図 A.2: 208 Pb に α 粒子を衝突させたときの散乱強度の比で与えられることがわかった。このことは、原子核の体積 V が

V =

4π 3
r0 A
3

となり、質量数 A に比例していることを示している。この性質を原子核の体積の飽和性という [12]。

A.2

核力
3

原子核の半径は 1.2A1/3 fm 程度であり、体積は V = (4π/3) (1.2) A fm3 ぐらいで、
3
中性子または陽子 1 個はほぼ (4π/3) (1.2) fm3 の体積を占めている。つまり、粒子は半径 1.2 fm の球の中に 1 個入っていることになる。核内での 2 核子間の距離は、多少の重なりを考慮に入れれば 、2 fm 程度である。陽子や中性子をこのくらいの距離に保ちつつ束縛させるもの ( 核力 ) の起源は何であろうか。1935 年、湯川秀樹は 2 粒子間で光子を交換するとクーロン力が生じることから類推して、2 核子間である未知の粒子が交換され核力が生じると考えた。いまこの粒子の質量を mµ とする。1 つの核子がこの粒子を放出すると、系のエネルギーは mµ c2 だけゆらぐ。このようなエネルギーのゆらぎ 、不確定さが持続する時間 τ は、不確定性関係に従えば τ ≈ h̄/mµ c2 の程度である。
この時間にこの粒子が走れる距離は

rµ =

h̄
mµ c
38

(A.5)

の程度である。(A.5) 式は質量 mµ なる粒子のコンプトン波長の 1/2π にほかならない。
rµ を 2 fm に等しいとすると、mµ は ( 核子のコンプトン波長の 1/2π が 0.21 fm である
ことを用いれば ) 核子の質量の約 1/10 倍である。事実、1947 年に Powell 、Occhialini 、
Lattes らが宇宙線にさらした乾板中に電子の 300 倍 ( 核子の 0.16 倍 ) の質量の粒子が200 倍くらいの粒子にこわれている飛跡を見出した。重い方は原子核と強く相互作用
し 、軽い方はほとんど 原子核と相互作用しなかった。重い方はパイ中間子、軽い方はミューオンと名づけられ、パイ中間子の質量は陽子の 1/6.7 倍であった。またパイ中間子には電荷が ±e のもの (π ± ) と、中性なもの (π 0 ) があり、それぞれの質量はmπ± c2 = 139.57 MeV, mπ0 c2 = 134.96 MeVである。従って π のコンプトン波長の 1/2π( それを λ̄) は

λ̄ = 1.4 fmとなる。
原点に e なる点電荷があるときのクーロンポテンシャル V (r) は∇2 V (r) = −4πeδ (r)(A.6)

を満す。量子力学では −h̄2 ∇2 = p2 であるから、この方程式の左辺の − (ch̄) 倍は c2 p2であり、それは光子のエネルギーの二乗である。パイ中間子は質量 mπ をもつから、そのエネルギーの二乗は c2 p2 + (mc2 ) である。そこで (A.6) 式にならってパイ中間子が生み出すポテンシャルは
2

{

(

mπ c
∇ −
h̄
2

)2 }

V (r) = 4πf δ (r)(A.7)と書けそうである。ここで f は核子とパイ中間子の結合の強さを表す定数である。この方程式 (A.7) の解は

V (r) = −fexp (−r/rπ )rになる。
2 陽子間のクーロンポテンシャルエネルギーが eV (r) であるように、2 核子間のパイ中間子交換力は

φ (r) = f V (r)exp (−r/rπ )= −f 2rと書ける。e−x /x の形を湯川型という。陽子の陽子による散乱の実験からこの f を決め
ると

f2= 0.081 ± 0.02h̄c39である。ここで h̄c で割ったのは次元のない数にするためである。一方、電荷 e に対す
る同じような式は、微細構造定数でe21== 0.0073h̄c137である。この 2 つの数を比較すればパイ中間子と核子の結合の方が電磁場より 10 倍も強いことがわかる。そこで前者を強い相互作用という [12]。

A.3

結合エネルギー

原子核は Z 個の陽子と N = A − Z 個の中性子からできているとすると、その質量はほぼ

Zmp + N mn

(A.8)

である。但し 、mp 、mn は陽子と中性子の質量である。ところが現実の原子核の質量
M (Z, N ) は (A.8) 式より小さい。それは核力による結合エネルギーのためである。これを BE(Z, N ) と書くとBE (Z, N ) = {Zmp + N mn − M (Z, N )} c2である。例えば重陽子は mp + mn より 2.22 MeV/c2 だけ軽くなっている。図 A.3 にはBE(Z, N ) を A で割った 1 核子あたりの結合エネルギーが示してある。これを見ると、
1620軽い原子核は一般に結合エネルギーが小さいが 、42 He 、126 C 、 O 、 Ne 、· · · などは例
外で大きい。1 核子あたりの結合エネルギーは鉄 Fe やニッケル Ni で最大になり、8.7
MeV に達する。この BE/A は A が 250 になると 7.5 MeV ぐらいまで減るが 、平均は 8
MeV ぐらいとみられる。そこで大ざ っぱに

BE (Z, N ) ' 8A MeVとしてよい。これは核内で 1 核子がもつ結合エネルギーはほぼ 8 MeV で質量数にほとんど よらないことを意味していて、この現象を原子核の結合エネルギーの飽和現象という。
実験によれば 、2 陽子、2 中性子状態には束縛状態がなく不安定であるが 、陽子と中性子からなる系は重陽子 21 H(≡ d) とよばれ 、安定である。その束縛エネルギーは()

BE 21 H = 2.22 MeVであり、中性子 − 陽子間の力の方が同種核子間の力より強い。また同じ種類の核子が多くなると、パウリの禁止則によって、後から核に取り入れられる核子の状態は他のものと違う必要があるので、運動エネルギーの高い状態に入らなければならない。そこで、なるべく中性子の数と陽子の数の差を小さくした方がエネルギーは低く、これらの効果によるエネルギーを対称エネルギーという。これをasym (N − Z)2 /A
40

図 A.3: BE(Z, N ) を A で割った 1 核子あたりの結合エネルギーと書く。asym < 0 ならこの項は明らかに N = Z のとき最大である。
原子核内の核子がもと軌道角運動量の総和を L 、スピン角運動量の総和を S とする
とき、I = L + S なる全角運動量を原子核のスピンという。実験によれば 、Z も N の
偶数の原子核のスピンは例外なく 0 である。この事実を説明するには、原子核の基底状態では 2 陽子は必ず対になって合成角運動量が 0 になり、2 中性子についても同じことがいえると考えればよく、2 つの同種核子間には角運動量を 0 にするような核力が働くと考える。この力を対相互作用といい、この結果ほぼ

N = 偶数,N = 偶数,N = 奇数,Z = 偶数Z = 奇数}Z = 偶数2δ

N = 奇数,Z = 奇数0δ

だけ余分なエネルギーがつけ加わる。通常、A = 奇数の場合を基準にして、N = 偶数、Z = 偶数のとき δ 、N = 奇数、Z = 奇数のとき −δ がつけ加 わるものとする。また δ は

δ=−

130
MeV
A

程度である。

41

図 A.4: 結合エネルギーの実験値と液滴模型による計算値との差

以上、すべての項を加えると半経験的質量公式 (Weizsacker-Bethe の式 )


 δ


偶偶核
(N − Z)ac Z
BE (A) = aV A + aS A + asym
+ 0+
A = 奇核1
2A
A3
 −δ 奇奇核
2

2
3

2

(A.9)

となる。ここで aV A は体積エネルギー、aS A は表面エネルギー、ac Z 2 /A1/3 は陽子間のクーロン力によるエネルギーを表している。aV 、aS 、asym 、ac をパラメータとして BE の実験値に一番合うように決めると、例えば

aV = 15.56 MeV,aS = −17.28 MeV,

asym = −46.57 MeV,

ac = −0.697 MeVというパラメータの組が得られる [12]。

A.4

魔法数と殻構造

液滴模型によれば 、核内では核子は互いに衝突しあい結合エネルギーを交換しあっていて、核子が 1 つの原子核に捕獲されるや否や、結合エネルギーは原子核全体に分配されて熱平衡状態に達すると考えられている。従って、原子核内で 1 つの核子が一定の軌道に沿って他の核子と衝突することなく運動することは、全く予想されなかった。
図 1.4 には、Z = 36 をもつ Kr 原子核の結合エネルギーの実験値と、半経験的な (A.9)
式との差が示してある。この図 A.4 を見ると、N = 50 のときこの差が他のときと大きく異なっていて、実験値が理論値より大きくなることがわかる。即ち、これは液滴状態より安定になることを意味している。また図 A.5 には同じ領域の偶偶核の第 1 励起状態が示してある。この図 1.5 から、N = 50 、N = 82 および Z = 50(Sn) の原子核の42

励起エネルギーが高いことがわかる。このようなことからZ} = 2, 8, 20, 28, 50, 82, 126
N(A.10)をもつ原子核は特に安定であることが結論される。この数を原子核の魔法数という [12]。

図 A.5: 第一励起エネルギー

A.2 節で述べたように核力は近距離力である。従って核内のある核子に働く他核子からの核力は、その核子の周辺の密度 ρ (r) に比例する。そこで 1 核子に働くポテンシャ
ル V (r) は
V (r) = V0 ρ (r)と書ける。ρ (r) は内部でほぼ一定、表面で急激に減少する。そこで、V (r) は図 A.6 に示したような形になり、これをさらに近似して、数学的に取り扱いやすい調和振動子
ポテンシャルに等しいとする。即ち

V (r) = −V0 +mω 2 r2
2(A.11)とする。ここで −V0 は一定の深さ、m は核子の質量、ω は調和振動子の角速度である。

1 つの核子が、(A.11) 式という平均ポテンシャル内を運動しているとすると、この核子に対するシュレデ ィンガー方程式は)(mω 2 2h̄2 2∇ − V0 +r φ (r) = Eφ (r)−2m2(A.12)

である。r2 = x2 + y 2 + z 2 と書くと、この方程式は x 、y 、z 方向の別々の部分に分けられる。そこで変数分離法を用いればφ (r) = φnx (x) φny (y) φnz (z)43図 A.6: 密度と平均ポテンシャル
(

)

3
h̄ω
2
N0 = nx + ny + nz
ε = N0 +であることがわかる。但し 、² = E + V0 である。ここで調和振動子に対するシュレディ
ンガー方程式の固有値(

E = n+

)

1
h̄ω,
2

n = 0, 1, 2, · · ·を用いた。
一番エネルギーの低い状態は、量子数が(nx , ny , nz ) = (0, 0, 0)で

3
ε = h̄ω
2をもつ。この状態にスピン上向き、下向き 2 個の陽子または中性子が入り、この軌道を満員にするのが Z = 2 、N = 2 が魔法数である理由である。このように、1 つの軌道が満員になることを原子の場合と同じく閉殻になったといい、次に N0 = 1 の殻では(1, 0, 0) ,(0, 1, 0) ,44

(0, 0, 1)

という 3 つの状態が縮退している。スピンの上下を考慮に入れると、この軌道には陽子または中性子が 6 個入れる。そこで、Z = 2 + 6 = 8 、N = 8 が魔法数になる。全く同様に、Z (N ) = 8 + 12 = 20 、20 + 20 = 40 、さらに 70 、112 が魔法数になることが予想される。

図 A.7: スピン軌道力による分離

図 A.8: 魔法数の説明

ところが (A.10) 式によれば 2 、8 、20 までは同じであるが 、40 の代わりに 50 、70 の代わりに 82 でなければならないし 、28 も魔法数であった。この矛盾を解決したのがMayer と Jensen であった。
1949 年、Mayer と Jensen は独立に 、重い原子の構造で重要な役割を果たしたスピ
ン軌道力が 、原子核内の平均ポテンシャルでも強いと考えれば 、魔法数を説明できることを発見した。スピン軌道力とはスピン磁気モーメントと磁場 H との相互作用エネ
ルギー

W =−

eh̄2 1 d
U (r) (S · L)m2 c2 r drのことで、スピン − 軌道相互作用ともいい、これによりエネルギー準位の分離が起こる。ここで説明を簡単にするために、スピン軌道力をξ (L · S)と書き、ξ は定数であるとする。この L は軌道角運動量、S はスピン角運動量である。
また調和振動子型から少しずれていることを考慮に入れて、平均ポテンシャルに aL2

45

という項がつけ加わり a < 0 とする。この 2 つの相互作用がつけ加わると、(A.12) 式は{}h̄2 2 mω 2 2−∇ +r + aL2 + ξ (L · S) ϕ = εϕ2m2となる。この式の固有状態は

1n, L, J = L ± , m
2という 4 つの固有値で指定される。そして固有値は()

3εnLJ = 2n + L +h̄ω + aL (L + 1)
2{

J = L + 12
+ξ L2
−ξ L+1J = L − 122である。原子の場合は ξ > 0 で J = L − 1/2 が低かったが 、原子核では ξ < 0 であるとすると、N = 3 に属する軌道は図 A.7 のように 0f7/2 、1p3/2 、1p1/2 、0f5/2 の 4 本の軌道に分離する。このうち 1 番低いのは 0f7/2 である。J = 7/2 の状態は m = ±7/2 、±5/2 、±3/2 、±1/2 をもつから、8 個の陽子または中性子で閉殻になる。20 にこの 8 を加えると、まさに 28 の魔法数が得られたのである。
それでは 50 や 82 の魔法数はど うであろうか。図 A.8 のように、N = 4 殻に属する
0g 状態がスピン軌道力で 0g9/2 と 0g7/2 に分離し 、0g9/2 は N = 3 の 0f5/2 や 1p1/2 と競争するくらい低くなる。N = 3 殻が閉じると陽子または中性子数が 40 になる。一方、
0g9/2 は 10 個の陽子または中性子が入る。そこで合わせて 40 + 10 = 50 の魔法数が得られる。全く同様に、0h11/2 が N = 5 殻から離れて N = 4 殻と一緒になることによって 82 が 、0i13/2 が N = 6 殻から離れて N = 5 殻に入ることによって 126 の魔法数が説明されるのである [12]。
このように、強いスピン軌道力が存在すれば魔法数が説明できることがわかる。ス
ピン軌道力の起源は、まだ明らかではないが 、本論文で取り上げる相対論的効果がその一つの候補である。

46

付 録B

相対論的波動方程式

本付録の執筆にあたっては文献 [9] を参考にした。

B.1

相対論的古典力学

電磁場 (ϕ, A) のなかにあって、静止質量 m と電荷 e をもつ古典粒子の力学的性質を復習する。
この粒子の速度を v で表すと

dr
dt

v≡

のように書ける。また、相対質量 M と狭義の運動量 π を

M≡√

m
,
1 − v2π ≡ Mvのように定義する。(M, π) は 1 つの 4 元ベクトル π µ であって、そのノルムは m2 、即ちM 2 − π 2 = m2

(B.1)

である。(M > 0)。場がないとき、この粒子は一様な直線運動を行う。つまり、v = 定数である。電磁場があるとき、この粒子がたど る軌道は方程式dπ= e [ε + v × H] ≡ F
dt

(B.2)

を満たし 、これが質点の相対論力学の基本方程式である。この ε は電場、H は磁場を表していて、ベクトル F はローレンツ力と呼ばれる。方程式 (B.2) から出る結果のうち

dM= (v · F ) = e (v · ε)dtd (r × π)
=r×F
dt

(B.3)
(B.4)

をあげておく。これはそれぞれ質量と狭義の運動量のモーメントの変化の法則を与える。この粒子の固有時間 τ を導入すれば 、方程式 (B.3) 、(B.4) の関係を共変形式に書き表すことができる。そうすると τ は
1

dτ = (dxµ dxµ ) 2 =47√

1 − v 2 dt

によって定義される。これにより 4 元速度

dxµ
u ≡
≡
dτ
µ

(

dt dt
,v
dτ dτ

)

(uµ uµ = 1)を定めると、これに m を掛けた積が狭義の運動量の 4 元ベクトルになる。即ち

π µ ≡ muµ ≡ (M, π)となる。方程式 (B.2) と (B.3) は、形式上共変な方程式

dπ µ= eF µν uν ,
dτ
あるいは

duµ
e
= F µν uνdτmに同値である。ここで F µν は電磁テンソルである。
これらの運動法則をラグランジュまたはハミルトンの形式論から導き出すこともできる。運動量 p はエネルギー E とともに 1 つの 4 元ベクトル pµ をつくり、これは

pµ = π µ + eAµ(B.5)という関係によって π µ に結びつけられる。即ちE = M + eϕ, p = π + eAである。また、関係 (B.5) と (B.1) 式に合わせて、ハミルトン関数はH ≡ eϕ +√()

p − eA2 + m2

(B.6)

であるとする。この定義から、ハミルトンの正準方程式

dr
π
=
,
dt
M

dp
= −e grad (ϕ − v · A)dtが導かれる。第 1 式は速度の定義で、第 2 式は方程式 (B.2) に同値である。このことは、
場 ε と H の定義式

ε = −∇ϕ −

∂A
,
∂x0

H = rot A

を使い、

dA
=
dt

(

)

∂
+ v · grad A∂tという事実を利用すれば 、容易に証明される [9]。

48

B.2

クライン -ゴルド ン方程式

電子に対する相対論的波動方程式の研究はスピンの存在によって複雑になるので、まずスピン 0 の 1 粒子、例えば π 中間子に対する相対論的波動方程式をつくってみる。このような粒子は内部自由度をもたず、その波動関数 Ψ は変数 r と t だけを含む。粒子の質量を m 、電荷を e で表し 、これが電磁ポテンシャル Aµ ≡ (ϕ, A) のなかにあると考える。
波動方程式を求めるには、経験的な方法をとり対応原理を使えば 、古典近似の成立
条件が満されるとき、古典運動法則にもど ることが保証される。
シュレデ ィンガー方程式の対応規則

E −→

i∂
,
∂t

p −→ −i∇(B.7)

を使う。pµ ≡ (E, p) とおけば 、この規則はもっと簡単に

pµ −→ i∂ µと書き表される。
ハミルトニアンを表す (B.6) 式から√

E = eϕ +(p − eA)2 + m2

(B.8)

が得られ 、これに規則 (B.7) を当てはめると、波動方程式()∂
i − eϕ Ψ =
∂t

[(

)

]1

2
2
1
∇ − eA + m2 Ψiが導き出される。この方程式には 2 つの重大な欠点があり、第 1 に空間座標と時間座標の間に対称性がないため、相対論的不変性およびそれから生じる諸性質が明瞭に現れていないこと、第 2 に右辺の演算子が平方根であるため、場 A が 0 であるときを除け
ば 、実際上扱いようがないことである。
これら 2 つの困難を避けるために 、対応の手続きの出発点として関係 (B.1) 式をとる。これからは(E − eϕ)2 − (p − eA)2 = m2(B.9)

が得られる。この関係は (B.8) 式に同値ではなく、もっと一般的な関係

E = eϕ ±√(p − eA)2 + m2に同値である。符号 + だけが実の古典解に対応するが 、符号 − は負の質量をもつ解を表し 、物理的意味は全くない。それゆえ関係 (B.9) 式を出発点にとれば 、負の質量に対する解を導入することになる。
これらのことを保留した上で、(B.9) 式に対応の手続きを施せば 、クライン − ゴルド
ン方程式

(
)2
(
)2
∂
1
 i
− eϕ −∇ − eA  Ψ = m2 Ψ

∂t

i

49

を得る。これはまた
(

)

[

]

Dµ Dµ + m2 Ψ ≡ (∂µ + ieAµ ) (∂ µ + ieAµ ) + m2 Ψ = 0という形にも書かれ 、この相対論的不変性は明瞭である。
この方程式の解釈については議論を簡単にするため、場が 0 の場合に限る。このとき方程式は単に
(

□+m

2

)

(

1 ∂2
□ ≡ 2 2 − ∆ ≡ ∂µ ∂ µ
c ∂t

Ψ=0

)

(B.10)と書き表される。ここで□はダランベール演算子である。
これは時間に関して二階の微分方程式である。それゆえ、ある時刻より後の各瞬間において Ψ が完全に決められるためには、初めの時刻における Ψ と ∂Ψ/∂t とが同時に知られていなければならない。この困難は、与えられた瞬間における系の力学的状態
は関数 Ψ だけによって表されるのではなく、2 つの関数 Ψ と ∂Ψ/∂t の組、あるいは 2
つの線形結合

Φ=Ψ+

i ∂Ψ
,
m ∂t

χ=Ψ−

i ∂Ψ
m ∂t

によって表されると考えると容易に乗り越えられる。このことは、系の状態が 2 つの成分 Φ 、χ をもつ 1 つの関数によって表されることと同じである。この波動関数は時間に関して一階の微分方程式にしたがう。この方程式をクライン -ゴルドン方程式から導き出すのは容易である。非相対論的極限において、この粒子のエネルギーは明らかに静止質量 m に等しいので

i

∂Ψ
' mΨ
∂t

が得られ 、従って χ ¿ Φ である。即ち 2 つの成分のうち一方が他方に比べて無視できるようになり、シュレデ ィンガーの非相対論的な理論がもう一度得られる。この理論によれば 、スピン 0 の 1 粒子の力学的状態は 1 つの成分から成る波動関数によって表されるのであった。
ところで、波動関数と解釈し うるためには、存在の確率密度 P と流れの密度 j を定
義し 、これらが連続の条件∂P+∇·j =0∂t(B.11)

を満たすようにする必要がある。記号 j µ ≡ (P, j) を使えば 、(B.11) はまた∂µ j µ = 0となる。
Ψ と Ψ∗ は両方とも方程式 (B.10) を満たすから

Ψ∗ (□Ψ) − (□Ψ∗ ) Ψ = 0
50(B.12)

が成り立ち、ダランベール演算子の定義を使えば 、これから∂µ [Ψ∗ (∂ µ Ψ) − (∂ µ Ψ∗ ) Ψ] = 0を得る。左辺の括弧のなかに比例するように j µ をとれば 、連続の条件は満たされる。
ここで非相対論的極限において普通の定義が得られるように比例定数を定めるためにjµ =i[Ψ∗ (∂ µ Ψ) − (∂ µ Ψ∗ ) Ψ] ,
2m

即ち
[

]

i
∂Ψ ∂Ψ∗P(r ,t) =
Ψ∗
−
Ψ ,
2m
∂t
∂t

j (r ,t) =1[Ψ∗ (∇Ψ) − (∇Ψ∗ ) Ψ]
2im

(B.13)と定義する。
(B.13) 式を見ると、密度 P (r, t) が正の定符号ではないことがわかる。ここにクライ
ン -ゴルドン方程式がもつ最大の難点の 1 つがある。
もう 1 つの難点は、上と無関係ではないが ¿ 負のエネルギーをもつ解 À の可能性
からもたらされる。例えば 、場のない方程式の解として平面波

Ψ = exp [−i (Et − p · r)]を求めるとき、この式を方程式 (B.10) 式に代入すれば 、条件√

E = ± p2 + m2√を得る。それゆえ、負のエネルギー − p2 + m2 をもつ解が存在する。このような解が生じる原因はすでに指摘したように、負の質量を理論に導入したことによる。
これらの困難を全て処理するには、Pauli と Weisskopf に従って、4 元ベクトル j µ の解釈と平均値の定義を修正することである。それによって理論を解釈しなおすと、ej µ
は流れの密度の 4 元ベクトルであり、そのうち eP (r, t) は電荷密度であるので、(B.12)
は電荷保存の方程式となる。これに反して、粒子の個数は保存されないのは、逆の電荷をもつ粒子の対が消滅したり生成したりする可能性と解釈される。このように解釈すれば 、この理論は電荷に関する理論となり、粒子に関する理論ではなくなる。これに反してデ ィラックの理論においては、正の定符号をもつ密度 P を定義することができる [9]。

B.3

ディラック方程式

ここでデ ィラックに従い、非相対論的量子力学から類推することによって、電子に対する相対論的波動方程式をつくることにする。
非相対論的理論における電子は二成分をもつスピノールによって表され 、回転によるスピノールの変換法則は角運動量 1/2 の変換法則である。これと同じように、相対論的理論における電子はいくつかの成分をもつ波動関数によって表され 、これらの成51

分はローレンツ変換から成る一層広い群に関して、はっきり決まった変換法則を満たさなければならない。波 Ψ の s 番目の成分を ψs (r, t) と表せば 、Ψ は 1 列から成る行列




ψ1


 ψ2 
Ψ= .  .. 
ψNの形に書き表される。
非相対論的な場合と同じように、与えられた時刻において波 Ψ は軌道変数 r と内部変数即ちスピン変数 s (s = 1, 2, · · · , N ) の関数とみなすことができる。このような波は
1 つの状態ベクトル | φ (t) i を表し 、この状態の [ 定められる ] 空間 ε は軌道変数の空間
ε(0) とスピン変数の空間 ε(s) のテンソル積ε = ε(0) ⊗ ε(s)である。波 Ψ は適当な表現においてこのベクトルを表すものとする。即ち

Ψ(r ,s,t) ≡ ψs(r ,t) = h rs | Ψ(t) iを表す。同じ類推を続けて、点 r に存在する確率密度をP(r ,t) =N∑|ψs |2s=1

によって定義する。これだけ仮定すれば 、波動方程式は必ずi∂Ψ= HD Ψ
∂t(B.14)という形をとり、HD は状態空間のエルミート演算子になる。実際、Ψ は各瞬間における電子の力学的状態を完全に定めるのであるから、波動方程式は時間に関して一階でなければならない。さらに、HD がエルミートであることによって、P (r, t) に関する定義の一貫が保証される。
ここで求めているのは相対論的方程式であるから、さらにこの方程式が空間座標と時間座標に関して、ある形式的な対称性を示すこと、即ち空間の変数についても一階であることだ。
ここでまず、場がないときの 1 電子の場合を考える。このときハミルトニアンは並進によって不変でなければならないから、r を含まない。それゆえ、前に述べた仮定を全て考慮すれば 、これをHD = α · p + βmという形に書き表すことができる。ここで p は対応法則 (B.7) によって定められた意味をもつ。即ち p = −i∇ である。そして α ≡ (αx , αy , αz ) と β はスピン変数だけに作用する 4 個のエルミート演算子を表す。記号 E ≡ i∂/∂t を採用すれば 、この波動方程式は[E − α · p − βm] Ψ = 0
52(B.15)

と書かれる。α と β を定めるには、対応原理にもとづく。即ち、この方程式の解はクライン -ゴルドン方程式[]

E 2 − p2 − m2 Ψ = 0(B.16)に従うことを求める。(B.15) 式に演算子 [E + α · p − βm] を左から掛けると、二階の方程式

E 2 −

∑(

αk

)2 (

pk

)2

− β 2 m2 −

k

∑(

)

α k α l + α l α k p k pl −

∑(

k<l

)



αk β + βαk mpk  Ψ = 0kを得る。
この方程式が (B.16) 式と同じになるためには
(

αk

)2

= 1,

αk αl + αl αk = 0

β 2 = 1,αk β + βαk = 0(k 6= l) ,(B.17)(B.18)という関係を満たせばよい。即ち、4 個の演算子 β 、α が互いに反交換し 、それぞれの二乗が 1 に等しければよい。エルミート行列 β 、α が関係 (B.17) 、(B.18) 式を満たすように選ばれたとき、方程式 (B.15) はデ ィラック方程式になる。
自由電子に関するこの方程式から、電磁場 (ϕ, A) のなかの 1 電子に関するデ ィラッ
ク方程式に移るにはE −→ E − eϕ,p −→ p − eAという置き換えを行う。[ ここで e は電子の電荷 (e < 0)]。その結果[(E − eϕ) − α · (p − eA) − βm] Ψ = 0を得る。即ち[()]∂
i − eϕ − α · (−i∇ − eA) − βm Ψ = 0∂tを得る。この方程式を (2.14) 式と比較すれば 、場が存在するときのデ ィラックのハミ
ルトニアンとしてHD = eϕ + α · (p − eA) + βmを得る [9]。

B.4

中心力ポテンシャル

中心力による静ポテンシャル V (r) のなかにあるディラック粒子のハミルトニアンはµ− 原子のとき

HD ≡ α · p + βm − V (r)53であり、原子核のときは引力を表す Lorentz スカラーポテンシャル V s(r) と、斥力を表す Lorentz ベクトルポテンシャル V 0(r) を用いてHD ≡ α · p + β (m − g sV s (r)) + g vV 0 (r)と表すことができ、この g v 、g sはそれぞれのポテンシャルの強さを表すパラメータである。これらのハミルトニアンは回転と反転によって不変であり[HD , J ] = 0, [HD , P ] = 0と表される。そのため決まった角運動量と偶奇性をもつ固有解を求める。
解 Ψ を次のように
(

Φ
χ

Ψ=

)

書く。この Φ, χ は r および z 軸方向のスピン成分 µ の関数で
(

Φ≡

ψ1
ψ2

)

(

χ≡

,

ψ3
ψ4

)

と表される。解 Ψ は J 2 , Jz , P に共通な固有関数であるとすると、各運動量を固定する
量子数を (JM ) で表す。L による相対論的な偶奇性を表す量子数 $ を
{

$=

1

+1 (偶奇性 (−)J+ 2 の状態に対して)
1−1 (偶奇性 (−)J− 2 の状態に対して)と定める。このとき
(

J

2

Φ
χ

)

(

= J (J + 1)

Φ
χ

)

(

,

JZ

Φ
χ

)

(

=M

Φ
χ

)

(

,

P

(0)

Φ
χ

)

(
J+ $
2

= (−)

Φ
−χ

)

(B.19)
を得る。スピン 1/2 と L 位の球面調和関数との結合によってつくられ た全角運動量M(JM ) の関数を yLJ(θ, φ, µ) とすると、これは偶奇性 (−)L をもつ関数である。ここで、L は軌道角運動量、J は 1 粒子の全角運動量、M はその z 成分である。また角運動量の結合法則により、L は 2 つの値

1L ≡ J + $,2

1
L0 ≡ J − $2(B.20)Mと yLM0 J は反対の偶奇性をもち、前者は偶奇性 (−)J+$/2 、しか取りえない。2 つの関数 yLJ
後者は偶奇性 (−)J−$/2 である。また (B.19) 式より、Φ は (r, θ, ϕ, µ) の関数で角運動量Mの積に等しい。同じように χ は r(JM ) と偶奇性 (−)J+$/2 をもつので、r の関数と yLJ
Mの関数と yL0 J の積に等しい。
J+$/2をもつ状態を表すとき、ΨM
ΨM$J は$J によって角運動量 (JM ) と偶奇性 (−)

1
ΨM
$J =
r

(

M
F (r) yLJiG (r) yLM0 J

54)(B.21)

という固有関数で表すことができる。この式の L と L0 は (B.20) 式によって与えられ 、F と G は r の動径波動関数である。
ここで固有値問題
M
HD ΨM$J = EΨ$J(B.22)を解く。それには、演算子 HD のなかで ¿ 角 À 変数と動径変数を分離するとよい。
動径運動量 pr と ¿ 動径速度 À αr を次のように導入する

pr ≡ −i

1 ∂
r
r ∂r

αr ≡ (α · r̂) =(B.23)
ρ1 (σ · r).
r

(B.24)

また、Pauli 行列の等式から(α · r) (α · p) = (σ · r) (σ · p) = (r · p) + i (σ · L) = rpr + i (1 + σ · L)が導き出される。この式に左から αr /r を掛け、αr2 = 1 を使えば

(α · p) ≡ αr)(i
pr + (1 + σ · L)r

(B.25)を得る。演算子 1 + (σ · L) は

1 + (σ · L) = J 2 +

1
− L2
4

2
と示される。ΨM$J に対する L の作用は演算子(

1
J + $β
2

)(

)

1
1 1
J + $β + 1 ≡ J (J + 1) + + $β (2J + 1)2
4 2の作用と同じである。これより

1M(1 + σ · L) ΨM$J = − $ (2J + 1) βΨ$J
2(B.26)となり、(B.25) 式と (B.26) 式を代入すれば 、固有値方程式 (B.22 は、µ− 原子の場合
αr pr −(i$ J + 21
r

) 



M
β  + mβ − V (r) ΨM$J = EΨ$J(B.27)という形になる。ここで固有関数をその表示式 (B.21) に、pr と αr をそれぞれの定義式(B.23) と (B.24) 式に置き換えてM= −yLM0 J(σ · r̂) yLJ
M(σ · r̂) yLM0 J = −yLJ

55

を使えば 、(B.27) 式の各項は(

1 ∂= (α · r̂) −irΨMr ∂r $J

α r p r ΨM
$J

(

=

0
σ · r̂
σ · r̂
0

i ∂
= −
r ∂r
i ∂
=
r ∂r
(

αr

i$ J + 12
r

)

(

(

i$ J + 12=(r

i$ J + 12

=

(

r2

i$ J + 12=

M−iG (r) yLJ−F (r) yLM0 J

)(
)(
)(

r2
(

mβΨM
$J = m
m
=
r

(

(

i ∂
−
r ∂r

M
iG (r) yLJ
F (r) yLM0 J

(

βΨM
$J

)[

)

)

,

)(

M
iG (r) yLJ−F (r) yLM0 J

1
r

(

M
F (r) yLJ
−iG (r) yLM0 J

I 0
0 −I

)(

0
−σ · r̂
σ · r̂
0

)

)]

)

0
σ · r̂
σ · r̂
0

I 0
0 −I

M
F (r) yLJiG (r) yLM0 J

)

)

1
r

(

M
F (r) yLJiG (r) yLM0 J

M
F (r) ylJiG (r) ylM0 J

)

)

,
M
F (r) yLJiG (r) yLM0 J))となる。従って、ΨM$J の上の二成分と下の二成分を比較することにより
(

−

1 d
G (r) +r dr


i$ J + 12

)

m
V (r)
E
F (r) −
F (r) = F (r)
r
r
) r

G (r) +r2(

1

i$ J + 2d− G (r) = (E − m + V (r)) F (r)
+
,
dr
r
(

i d
F (r) +r dr

i$ J + 12
r2)

F (r) −(

imiV (r)iE
G (r) −
G (r) =
G (r)
r) 
r
r

1

i$ J + 2d F (r) = (E + m + V (r)) G (r)+drrが得られる。まとめると µ− 原子の連立一階微分動径方程式 [9] は
(



1

)

$ J+2
d
 G (r) = (E − m + V (r)) F (r) ,
−
+
dr
r
(




d
+
dr

$ J + 12r

(B.28)) F (r) = (E + m + V (r)) G (r)56(B.29)

と書け、このポテンシャル V (r) は∫V (r) =

ρp(r0 )|r − r 0|

dr 0である。これらと同様に原子核について計算すると、原子核の連立一階微分動径方程
式は
(



1

)

$ J+2
d
 G (r) = (E − V − m − V ) F (r)
−
+
V
s
,
dr
r
(




d
+
dr

$ J + 12r)

 F (r) = (E − V + m + V ) G (r)
V
s

と書ける。

57

付 録C

全角運動量の固有波動関数

スピン 1/2 と L 位の球面調和関数との結合によってつくられた全角運動量 (JM ) の
関数は
∑

M
=
yLJ

(LML 1/2 MS |JM ) YLML ξMS

ML MSと表すことができる。ML は軌道角運動量の z 成分、MS はスピンの z 成分であり、(LML 1/2 MS |JM ) は Clebsch-Gordan 係数といわれ、YLML は球面調和関数、ξMS はス
ピンの波動関数である。
Clebsch-Gordan 係数とは、L と 1/2 の合成によって J の状態をつくるとき、その波
動関数の係数のことである。

L2 YLML = L (L + 1) h̄2 YLMLLz YLML = ML h̄YLML
に対して
(
2

2

S ξMS = S (S + 1) h̄ ξMS ,

1
S=
2

)

Sz ξMS = MS h̄ξMSであればよい。これらに対して

L = r × p,
(

σx =

0 1
1 0

)

(

,

σy =

1
S = h̄σ,
2

0 −i
i 0

)

(

,

σz =

1 0
0 −1

)

とおけば

[Lx , Ly ] = ih̄Lz ,[Ly , Lz ] = ih̄Lx ,[Lz , Lx ] = ih̄Ly[Sx , Sy ] = ih̄Sz ,[Sy , Sz ] = ih̄Sx ,[Sz , Sx ] = ih̄Sy

と同様に

58

を満足するので 、上の行列表式を採用する。σ は Pauli のスピン行列という。そうす
ると
(

ξ1 =
2

1
0

(

)

,

ξ− 1 =
20
1)と書けばよく、これをスピナー (スピノール ) という。このようにおくと

1
Sz ξ 1 = h̄ξ 12
2 2となり、ξ1/2 は Sz の固有値 h̄/2 をもつ関数である。ここで

σx 2 = σy 2 = σz 2 = 1

より、S は

1 2 2
h̄ σ
4
= S (S + 1) h̄2S2 =である。次に J 2 , Jz の固有関数をつくる。そのために

YLML = | LML i,
ξMS = | 1/2 MS i(

1
MS = ±
2)と書き、まず | LL i| 1/2 1/2 i を調べる。ここでL± = Lx ± iLy , S± = Sx ± iSy
とおくとJ 2 = (L + S)2= L2 + S 2 + 2LS= L2 + S 2 + 2Lz Sz + (L+ S− + L− S+ )と書ける。軌道角運動量について

Lz YLML = ML h̄YLML ,L± | LML i =√(L ∓ ML ) (L ± ML + 1)h̄| LML ± 1 iであったのと同様に、スピンについても
√(

)(

)

√(

)(

)

1 1
1 1
S+ | 1/2 − 1/2 i =+− + 1 h̄| 1/2 1/2 i
2 2
2 2= 1h̄| 1/2 1/2 i,
1 1
1 1
+
− + 1 h̄| 1/2 − 1/2 i
2 2
2 2= 1h̄| 1/2 − 1/2 i

S− | 1/2 1/2 i =59

となるべきである。これらを使えばJ 2 | LL i| 1/2 1/2 i =()

L2 + S 2 + 2Lz Sz + L+ S− + L− S+ | LL i| 1/2 1/2 i
[
]
13 2
1 2
2
= L (L + 1) h̄ +h̄ + 2L h̄ | LL i| 1/2 1/2 i
22
2
(
)
3
= L2 + L + + L h̄2 | LL i| 1/2 1/2 i
4
(
)(
)
1
3 2
= L+
L+
h̄ | LL i| 1/2 1/2 i22となる。従って | LL i| 1/2 1/2 i は、J = L + 1/2 の固有値をもつ J 2 の固有関数であり、さらにJz | LL i1/2 1/2 = (Lz + Sz ) | LL i| 1/2 1/2 i
)
(
1
= L+
h̄| LL i| 1/2 1/2 i2より、固有値 (L + 1/2) h̄ をもつ Jz の固有関数でもある。このときは
M =L+ 1yL J=L+21 = | LL i| 1/2 1/2 i
2= YLL ξ 1
2である。次に J− = L− + S− をこれに演算してみると

J− | LL i| 1/2 1/2 i = (L− + S− ) | LL i| 1/2 1/2 i√
2Lh̄| LL − 1 i| 1/2 1/2 i + h̄| LL i| 1/2 − 1/2 i== | A i.
| A i に対して、同様にすると
J 2 | A i = J (J + 1) h̄2 | A i, J = L +Jz | A i =

(

)

1
2

1
L−
h̄| A i
2となり、即ち | A i は J 2 と Jz の固有関数である。さらに J− | A i = | B i を求めると | B i
は J 2 と Jz の固有関数で固有値はJ 2 | B i = J (J + 1) h̄2 | B i, J = L +Jz | B i =

(

)

1
2

3
L−
h̄| B i2となるので、この操作を繰り返せば

1
2
(
)
1
3
1
1
M = L + ,L − ,L − ,···,− L +2222
J = L+

60

M
の yLJ
が求まる。次に規格化をする。例として | A i について

B 0 = C| A iとおき、規格化をすると

h B0 | B0 i = 1
でなければならない。ここで

h 1/2 1/2 | 1/2 1/2 i = 1,
h LL − 1 | LL i = 0,h LL | LL i = 1
であるから

C = ±√

1
2L + 1となる。ここでは + をとることにする。よって
M =L− 1

| B 0 i = yL J=L+21
2

(√
)
1
= √
2L| LL − 1 i| 1/2 1/2 i + | LL i| 1/2 − 1/2 i
2L + 1√
2L
1=
YLL ξ− 1
YLL−1 ξ 1 + √2
2
2L + 1
2L + 1となる。ここで√2L= h LL − 1 1/2 1/2 | L + 1/2 L − 1/2 i
2L + 1
1√= h LL 1/2 − 1/2 | L + 1/2 L − 1/2 i
2L + 1とおき、これらの係数を C.G.(Clebsch-Gordan) 係数という。この C.G. 係数は

φ=∑Cα | α iのように、全波動関数 φ を | α i で表したときの係数だから∑

Cα 2 = 1α又は∑(LML 1/2 MS |JM )2 = 1

ML MS

を満足する。

61

(C.1)

付 録D

4 次の Runge-Kutta 法

4 次の Runge-Kutta 法と呼ばれる数値積分法を説明する。まず次のような一階の微分方程式を解くことを考える
dy= f (x, y).
dx

(D.1)

4 次の Runge-Kutta 法は以下の操作を繰り返して、微分方程式を数値的に求める方法
である。図 6.1 は Runge-Kutta 法を図示したものである。手順をおって説明する。

図 D.1: Runge-Kutta 法

1. 点 (x0 , y0 ) で傾き f (x0 , y0 ) で ∆x だけ進むと y は k1 だけ変化することになる。
2. この線上の点 (x0 + ∆x/2, y0 + k1 /2) で傾きを計算し 、(x0 , y0 ) から ∆x だけ進むと y は k2 だけ変化することになる。
3. 同じように (x0 , y0 ) から傾き (x0 + ∆x/2, y0 + k2 /2) で ∆x だけ進むと k3 の変化となる。
4. (x0 , y0 ) から傾き (x0 + ∆x, y0 + k3 ) で ∆x だけ進むと k4 の変化となる。
5. 以上の計算結果の加重平均をとって点 (x1 , y1 ) を得る。
これはちょうど 、x = x0 における y(x0 + ∆x) の 4 次までのテーラー展開に一致して
いる。

62

付 録E

波動関数の漸近形の導出

本付録では、核子波動関数の r = 0 でのテーラー展開および r → ∞ での漸近形の導出法を詳述する。ε = E − VV (r) 、µ = m + Vs (r) 、κ = $ (J + 1/2) とおくと、(B.27)
式は
[

]

d
h̄c
κ
F (r)
G (r) =+
ε + µ dr r(E.1)

と書ける。(E.1) 式を (B.28) 式に代入すると[]

d2
κ (κ + 1)(ε2 − µ2 )
− 2+
F
(r)
=
F (r),
dr
r2
h̄c

即ち
[

]

d2
κ (κ + 1) (ε2 − µ2 )−+F (r) = 0
dr2
r2
h̄c

(E.2)

が求まる。ここで $ = 1 のとき κ = (J + 1/2) > 0 で、L = J + 1/2 、即ち、κ (κ + 1) =L (L + 1) 、よって κ = L である。一方、$ = −1 のとき κ = − (J + 1/2) < 0 で 、
L = J − 1/2 、即ち κ (κ + 1) = L (L + 1) であるから κ = −L − 1 である。従って (E.2)
式は
[

]

d2
L (L + 1) (ε2 − µ2 )−+F (r) = 0
dr2
r2
h̄c

となる。
(i) r → 0 の場合(E.3) 式は ε2 − µ2 < 0 のとき原点で正則な解はないので
ε2 − µ2 > 0
でなければならない。k を任意の定数としk2 =(ε2 − µ2 )
h̄c

とおくと
[

]

L (L + 1)
d2−+ k 2 F (r) = 0
2
2
dr
r
63

(E.3)

の原点で正則な一般解はF (r) = c0 rjL (kr)(E.4)

である。c は定数、jL は球 Bessel 関数である。よって、r → 0 では

F (r) ∝ rL+1

(E.5)

となる。この式を (E.1) 式に代入すると

G (r) =h̄c[(L + 1) + k] rLε+µが求まる。$ = 1 のときは κ = L であるから

G=h̄c (2L + 1) Lrε+µ(L = J + 1/2)となる。一方、$ = −1 のときは κ = − (L + 1) であるから、G = 0 となってしまう。
そこで、(E.4) 式を (E.1) 式に代入して計算すると
[

]

d
h̄c
L+1 0
G (r) =−
c rjL (kr)ε + µ dr
r
[
]
h̄cc0
djL (ρ)=ρ− LjL (ρ), (ρ = kr)
ε+µ
dρ

(E.6)

のようになる。ここで djL (ρ) /dρ は、公式を用いるとdjL (ρ)h̄c=[LjL−1 (ρ) − (L + 1) jL+1 (ρ)]
dρ
2L + 1と表され 、(E.6) 式は
[

]

L
L+1
h̄cc0
G (r) =
ρjL−1 (ρ) −ρjL+1 (ρ) − LjL (ρ)ε + µ 2L + 1
2L + 1と書ける。まず、L 6= 0 のとき

LL
ρjL−1 (ρ) = LjL (ρ) −ρjL+1 (ρ)
2L + 1
2L + 1であるから、(E.7) 式は

h̄cc0(−) ρjL+1 (ρ)
ε+µ
h̄cc0
1
= −
kr
ρL+1 + · · ·ε + µ (2L + 3)!!
1
k L+2 L+2
= −h̄cc0r+ ···(2L + 3)!! ε + µ

G (r) =

64

(E.7)

のように書ける。次に、L = 0 のときは、(E.7) 式から

h̄cc0(−) ρj1 (ρ)
ε+µ
k2 r2
= −h̄cc0
+ ···
ε+µ 3

G (r) =となる。
以上から、$ = −1, L = J − 1/2 のときは F (r) = rL+1 とすれば 、

h̄ck 2 L+2r+ ···2L + 3 ε + µ
1 ε − µ L+2= −r
2L + 3 h̄c

G (r) ' −

(E.8)

と表せる。
(ii) r → ∞ のとき
このとき (B.28) 式、(B.29) 式は−dG (r)
E−m=
F (r) ,drh̄c
dF (r)
E+m
=
G (r)
dr
h̄c

(E.9)
(E.10)と書ける。(E.10) 式を変形してG (r) =

h̄c dF (r)E + m dr(E.11)と書き、(E.9) 式に代入すると

d2 F (r)m2 − E 2
=
F (r)
dr2
h̄c
が求まる。従って√

F = e− m −E r/h̄c22(E.12)と書くことができる。これを (E.11) 式に代入すると√
h̄cm2 − E 2 −√m2 −E 2 r/h̄c
G = −
e
m+E
h̄c
√
m − E −√m2 −E 2 r/h̄c
e
= −
m+E
となる。
以上をまとめると、(i) r → 0 のとき

F (r) = rL+1
G (r) = h̄c(2L+1) r L −ε+µ
1 ε−µ L+2r
2L+3 h̄c65($ = 1)($ = 1)(E.13)

(ii) r → ∞ のとき√

F (r) = e− m −E r/h̄c√

G (r) = −22m − E −√m2 −E 2 r/h̄cem+E以上の境界条件を用いる。

66

付 録F

プログラムリスト/*
Reletivistic mean-field model for spherical symmetric systems
2007/1/15: rmf3a.c (completed and renamed)
2007/11/1,2,6,7 : rmf2c.c (self consistent field calculation)
2007/10/10-23,31: rmf2b.c (given Woods-Saxon form potential)
2007/9/3: renamed from rmf1-z.c to rmf2a.c
2007/5/21: interpolate density
2007/5/7: both for protons and neutrons
2007/4/17: proton density after sorting levels
2007/2/7: sort E
2007/1/10: proton density in main
2006/12/27: prepareDp
2006/12/22: Ebig,global Dp
2006/12/8,10-16 : created*/#include <stdio.h>#include <stdlib.h>#include <math.h>// Arrays to store the energies and wavefunctions of the eigenstates#define NEUTRON 0 // idxt = 0 for neutrons#define PROTON 1 // idxt = 1 for protons#define SIGMA 0// idxb = 0 for sigma meson#define OMEGA 1// idxb = 1 for omega meson#define RHO 2// idxb = 2 for rho meson#define COULOMB 3 // idxb = 3 for photon#define NBOSON 4// 0 <= idxb < NBOSON#define IDXT 2// idxt (index t) = 0..IDXT//
idxt = NEUTRON (=0) : neutron//
idxt = PROTON (=1) : proton#define IDXP 2// idxp (index p) = 0..IDXP-1//
idxp = 0 : varpi=-1, L=J-1/2//
idxp = 1 : varpi= 1, L=J+1/2#define IDXJ 8// idxj(index j)=0..IDXJ-1//
J = idxj+1/2#define NODE 5// node=0..NODE-1//
node=0 for the gournd state#define IDXR 3001 // idxr(index r)=0..IDXR-1

67

// "s" stands for "stored in memory"
double Es[IDXT][IDXP][IDXJ][NODE];// energy (MeV)
double Fs[IDXT][IDXP][IDXJ][NODE][IDXR];// large component of the radial wavefunction
double Gs[IDXT][IDXP][IDXJ][NODE][IDXR];// small component of the radial wavefunction// r = idxr*dr : idxr = 0..IDXR-1
double Vs[IDXT][2*IDXR-1];// scalar potential (MeV) Vs[i] = Vs(r=dr*i/2)
double Vv[IDXT][2*IDXR-1];// vector potential (MeV) Vv[i] = Vv(r=dr*i/2)
double Dens[2*IDXT][2*IDXR-1];// Dens[0][idxrd] : neutron |F|^2+|G|^2// Dens[1][idxrd] : proton |F|^2+|G|^2// Dens[2][idxrd] : neutron |F|^2-|G|^2// Dens[3][idxrd] : proton |F|^2-|G|^2// r = idxrd*dr/2 : idxrd = 0..2*IDXR-2
double phi[NBOSON][2*IDXR-1];// classical boson fields
double source[NBOSON][2*IDXR-1];// source term of boson fields// physical constants
const double HbC = 197.326960;// h-bar c [MeV fm^2]
const double fineStructureConstant = 1.0/137.036000; // fine structure constantconst double Ebig = 9.99e+30;// global variables
int idxtGlobal;
int NZ[IDXT]={126,82}; // the numbers of neutrons and protons (N and Z)
double Mass[IDXT] = {939.0, 939.0/*939.56533, 938.27200*/};// neutron and proton mass [MeV/c^2]
double BosonMass[NBOSON] = {520.0, 783.0, 770.0, 0.0}; // boson mass [MeV/c^2]
double squaredCoupling[NBOSON] = {109.626, 190.431, 16.3065/*65.226*/, 0.0}; // [hbar*c]// N.B.
16.3065=65.226/4
double Coupling[NBOSON]; // Coupling constant, to be defined in setParameters
double Rin;// (fm) The radius of connection to the Taylor series around r=0.
double Rout; // (fm) The radius of the connection to an asymptotic solution
for large r.
double Rmatch; // (fm) The radius where the forward and backward solutions meet.
double RadialGridSpacing; // (fm) Copy to a local variable "dr" and use "dr".

// paramters to construct the initial potential
double VsStr[IDXT] = {-373.0,-373.0}; // 398.5 ;// (MeV) Strength of the scalar potential
double VvStr[IDXT] = {274.0,274.0}; // 340.0 ;// (MeV) Strength of the vector potential// Positive (negative) strengths mean repulsive (attractive) potentials.
double R_ws=7.0; // (fm) radius of the Woods-Saxon potential
double r0_ws=1.2; // (fm) R=r0_ws*A^(1/3.0)
double a_ws=0.7; // (fm) surface thickness of the Woods-Saxon potential
struct RadialGridIndex {int i ; // inside point, where outward solution starts
int m ; // matching point, where inward and outward solutions meetint o ; // outside point, where the inward solution starts

68} Idxr;// NB) IDXR is macro, Idxr is structure, idxr is int

typedef struct{
double E ;int idxp ;int idxj ;int idxn ;} orbital;
orbital e[IDXT][IDXP*IDXJ*NODE];// prototype declarations of functions
int setParameters();
int writeParameters(FILE *FO);int prog1();
int prog2(int idxt, orbital *e, double *etot, int idxp, int idxj, int idxn);
int writeWavefunction(FILE *FO,int idxt, int idxp, int idxj, int idxn);// int swap : defined before function sort, only by which it is called.
int sort(orbital *e);
int initialPotential(int sw);
int idx2qn(int sw,int idxp,int idxj, double *J,int *w,int *L,int *kappa);
int setRadialGrid(int sw);
double solveDirac(int idxp,int idxj,int idxn);// retuned value = Energy eigenvalue (MeV)// int rungeBCi : defined before function runge, only by which it is called.
// int rungeBCo : defined before function runge, only by which it is called.
int runge(int sw, double E, int idxp, int idxj, int idxn, int *node, double*maco);
int Poisson(int imax,double h,double *source,double *phi);
int screenedPoisson(int imax,double h,double m,double source[],double phi[]);
int interpolate(double *y, int n);//------------------------------------------------------------------------------------int main(){setParameters();
writeParameters(stdout);
prog1();}//------------------------------------------------------------------------------------int setParameters(){
int idxb;
for(idxb=0;idxb<NBOSON;idxb++){Coupling[idxb]=sqrt(squaredCoupling[idxb]*HbC);}
Coupling[COULOMB]=sqrt(4*M_PI*HbC*fineStructureConstant);}//------------------------------------------------------------------------------------int writeParameters(FILE *FO)
{ // writes the physical parameters of the model to a stream

69

fprintf(FO,"neutron
mass=%10.4f (MeV/c^2)\n",Mass[NEUTRON]);
fprintf(FO,"proton
mass=%10.4f (MeV/c^2)\n",Mass[PROTON]);
fprintf(FO,"sigma meson mass=%10.4f (MeV/c^2) g^2=%10.4f (hbar*c)\n",BosonMass[SIGMA],squaredCoupling[SIGMA]);
fprintf(FO,"omega meson mass=%10.4f (MeV/c^2) g^2=%10.4f (hbar*c)\n",BosonMass[OMEGA],squaredCoupling[OMEGA]);fprintf(FO,"rhomeson mass=%10.4f (MeV/c^2) g^2=%10.4f (hbar*c)\n",BosonMass[RHO ],squaredCoupling[RHO ]);fprintf(FO,"hbar*c=%12.6f (MeV*fm^2)\n",HbC);
fprintf(FO,"fine structure constant=1/%12.6f\n",1.0/fineStructureConstant);
fprintf(FO,"Coupling Constant : sigma=%.6e omega=%.6e rho=%.6e coulomb=%.6e\n",Coupling[SIGMA],Coupling[OMEGA],Coupling[RHO],Coupling[COULOMB]);// fprintf(FO,"vector pot=%10.3f scalar pot=%10.3f (MeV)\n",VvStr,VsStr);// fprintf(FO,"Woods-Saxon parameters: R=%10.5f a=%10.5f\n",R_ws,a_ws);}//------------------------------------------------------------------------------------int prog1(){/*
Calculation of levels and densities fulfilling selfconsistency condition*/
int idxt,i,idxp,idxj,idxn,idxr,idxrd,degeneracy,num,iter;
int imax=IDXP*IDXJ*NODE;
double E,sumspe[IDXT],eboson,etot;// Energy [MeV], sum of single-particle energy
double r,dr,drb; // radial coordinate (fm) and its grid spacing (drb=drb/2)
int iprint;
double diffVs[IDXT],diffVv[IDXT]; // maximum (for r) of the discrepancy in potentials
double diff,diff2;
double *Fsptr, *Gsptr;
FILE *FFD, *FBD, *FPO, *FWF;
FFD=fopen("fdens.dat","w"); // Fermion (=nucleon) densities
FBD=fopen("bdens.dat","w"); // Boson densities
FPO=fopen("pot.dat","w");// Potential energies
FWF=fopen("wf.dat","w");// Nucleon wavefunctions
printf("# N=%d Z=%d\n",NZ[NEUTRON],NZ[PROTON]);
setRadialGrid(1);initialPotential(1);dr=RadialGridSpacing;
drb=dr/2;
fprintf(FPO,"# iter %d (initial potential)\n",0);iprint=0.1/drb; if(iprint<1) iprint=1;
for(i=0;i<2*IDXR-1;i++){if(i<=30||i % iprint == 0){
fprintf(FPO,"%9.5f %.7e %.7e %.7e %.7e %.7e %.7e %.7e %.7e\n",i*drb,Vs[NEUTRON][i],Vv[NEUTRON][i],Vs[PROTON][i],Vv[PROTON][i],0.0, 0.0, 0.0, 0.0);

70

}
}
fprintf(FPO,"\n");for(iter=1;iter<=1000;iter++){ // loop for selfconsistency : begin
for(idxt=0;idxt<IDXT;idxt++) prog2(idxt,e[idxt],&sumspe[idxt],idxp,idxj,idxn);
printf("iter=%d prog2 finished. E=%e %e\n",iter,sumspe[NEUTRON],sumspe[PROTON]);
printf("# J=idxj+0.5, L=idxj+idxp, varpi=2*idxp-1, kappa=varpi*(J+0.5)\n");for(idxt=0;idxt<IDXT;idxt++){num=0;
for(i=0;i<imax;i++){if(e[idxt][i].E < Ebig){degeneracy=2*e[idxt][i].idxj+2; num+=degeneracy;
printf(" #1#%3d idx(t,p,j,n)=%d %d %d %d E=%9.4f m= %+10.5f dgn=%3d sum=%3d\n",i,idxt,e[idxt][i].idxp,e[idxt][i].idxj,e[idxt][i].idxn,e[idxt][i].E,e[idxt][i].E-Mass[idxt],degeneracy,num);if(0){
fprintf(FWF,"# %d %d %d %d %d\n",idxt,i,e[idxt][i].idxp,e[idxt][i].idxj,e[idxt][i].idxn);
Fsptr=Fs[idxt][e[idxt][i].idxp][e[idxt][i].idxj][e[idxt][i].idxn];
Gsptr=Gs[idxt][e[idxt][i].idxp][e[idxt][i].idxj][e[idxt][i].idxn];
for(idxr=0;idxr<IDXR;idxr++){fprintf(FWF,"%9.5f %.5e %.5e\n",idxr*dr,Fsptr[idxr],Gsptr[idxr]);}
fprintf(FWF,"\n");}}//>if(e)}//>for(i)}//>for(idxt)
if(0){
for(idxrd=0;idxrd<2*IDXR-1;idxrd++){r=idxrd*drb;
fprintf(FFD,"%9.5f %.7e %.7e %.7e %.7e\n",r,Dens[0][idxrd], Dens[1][idxrd], Dens[2][idxrd], Dens[3][idxrd]);}fprintf(FFD,"\n");}// boson field
for(i=0;i<2*IDXR-1;i++){source[SIGMA][i] = (Dens[2][i]+Dens[3][i])*Coupling[SIGMA];
source[OMEGA][i] = (Dens[0][i]+Dens[1][i])*Coupling[OMEGA];
source[RHO][i] = (Dens[0][i]-Dens[1][i])*Coupling[RHO];
source[COULOMB][i] = Dens[1][i]*Coupling[COULOMB];

71

}screenedPoisson(2*IDXR-1,drb,BosonMass[SIGMA]/HbC,source[SIGMA],phi[SIGMA]);
screenedPoisson(2*IDXR-1,drb,BosonMass[OMEGA]/HbC,source[OMEGA],phi[OMEGA]);
screenedPoisson(2*IDXR-1,drb,BosonMass[RHO]/HbC,source[RHO],phi[RHO]);
Poisson(2*IDXR-1,drb,source[COULOMB],phi[COULOMB]);if(0){
for(i=0;i<2*IDXR-1;i++){r=i*drb;
fprintf(FBD,"%9.5f %.7e %.7e %.7e %.7e\n",r,phi[SIGMA][i],phi[OMEGA][i],phi[RHO][i],phi[COULOMB][i]);}fprintf(FBD,"\n");}// total energy
eboson=0;
for(i=1;i<2*IDXR-1;i++){eboson+=i*i*(-phi[SIGMA ][i]*source[SIGMA ][i]+phi[OMEGA ][i]*source[OMEGA ][i]+phi[RHO][i]*source[RHO][i]+phi[COULOMB][i]*source[COULOMB][i]);}eboson*=4*M_PI*drb*drb*drb;
etot=sumspe[NEUTRON]+sumspe[PROTON]-eboson/2;
printf("# total energy(MeV)=%.6f (%.6f %.6f %.6f)\n",etot,sumspe[NEUTRON],sumspe[PROTON],eboson);// new potential for the next iteration
for(idxt=0;idxt<IDXT;idxt++){diffVs[idxt]=0; diffVv[idxt]=0;}
fprintf(FPO,"# iter %d\n",iter);for(i=0;i<2*IDXR-1;i++){double Vsigma, Vomega, Vrho, Vcoul, Vsnew[IDXT], Vvnew[IDXT];
const double p=0.5; // damping factor
Vsigma=-Coupling[SIGMA] *phi[SIGMA][i] ;Vomega= Coupling[OMEGA] *phi[OMEGA][i] ;Vrho = Coupling[RHO]*phi[RHO][i];Vcoul = Coupling[COULOMB]*phi[COULOMB][i];Vsnew[NEUTRON]=Vsigma;
Vsnew[PROTON ]=Vsigma;
Vvnew[NEUTRON]=Vomega+Vrho;
Vvnew[PROTON ]=Vomega-Vrho+Vcoul;if(i<=30||i % iprint == 0){
fprintf(FPO,"%9.5f %.7e %.7e %.7e %.7e %.7e %.7e %.7e %.7e\n",i*drb

72

,Vsnew[NEUTRON],Vvnew[NEUTRON],Vsnew[PROTON],Vvnew[PROTON],Vsigma,Vomega,Vrho,Vcoul);}for(idxt=0;idxt<IDXT;idxt++){diff=fabs(Vsnew[idxt]-Vs[idxt][i]); if(diffVs[idxt]<diff) diffVs[idxt]=diff;
diff=fabs(Vvnew[idxt]-Vv[idxt][i]); if(diffVv[idxt]<diff) diffVv[idxt]=diff;Vs[idxt][i]=(1-p)*Vs[idxt][i]+p*Vsnew[idxt];
Vv[idxt][i]=(1-p)*Vv[idxt][i]+p*Vvnew[idxt];}}
diff = diffVs[NEUTRON] ;diff2 = diffVs[PROTON] ; if ( diff2 > diff ) diff = diff2 ;diff2 = diffVv[NEUTRON] ; if ( diff2 > diff ) diff = diff2 ;diff2 = diffVv[PROTON] ; if ( diff2 > diff ) diff = diff2 ;printf("# max (for r) potential inconsistency(Vs(n),Vv(n),Vs(p),Vv(p),all) (MeV)");
printf("%e %e %e %e %e :PD\n",diffVs[NEUTRON],diffVv[NEUTRON],diffVs[PROTON],diffVv[PROTON],diff);fprintf(FPO,"\n");if(diff<1.0e-8) {printf("# selfconsistency fulfilled: iter %d, potential difference<%e\n",iter,diff);
break;}} // for(iter): loop for selfconsistency : endfor(i=0;i<2*IDXR-1;i++){r=i*drb;
fprintf(FFD,"%9.5f %.7e %.7e %.7e %.7e\n",r,Dens[0][i], Dens[1][i], Dens[2][i], Dens[3][i]);fprintf(FBD,"%9.5f %.7e %.7e %.7e %.7e\n",r,phi[SIGMA][i],phi[OMEGA][i],phi[RHO][i],phi[COULOMB][i]);}fprintf(FFD,"\n");fprintf(FBD,"\n");fclose(FFD); fclose(FBD); fclose(FPO); fclose(FWF);
return 0;}//-------------------------sort energy and calculate density--------------------------int prog2(int idxt, orbital *e, double *etot, int idxp, int idxj, int idxn){
int i,n,imax,fn,idxr,idxrd;
double r, dr, E, *Fsptr, *Gsptr;
idxtGlobal=idxt;73

// energy sweep (for a test of this code)
if(0){
int mxi=200, node ; double maco;
for(i=0;i<mxi;i++){
E=Mass[idxt]-mxi+i;
runge(0,E,idxp,idxj,idxn,&node,&maco);}}// calculation of the energies and wavefunctions of the eigenstates// for all the combinations of the quantum numbers
i=0;
imax=IDXP*IDXJ*NODE-1;
for(idxn=0;idxn<NODE;idxn++){for(idxj=0;idxj<IDXJ;idxj++){for(idxp=0;idxp<IDXP;idxp++){E=solveDirac(idxp, idxj, idxn);
e[i].E=E;
e[i].idxp=idxp;
e[i].idxj=idxj;
e[i].idxn=idxn;i++;} //>for(idxp)} //>for(idxj)} //>for(idxn)// sorting the eigenstates in the asceding order of energysort(e);// calculation of the densities
double f6 = 4.0/3.0;
double f7 = -1.0/3.0;
dr=RadialGridSpacing/2;
fn=NZ[idxt];
for(idxrd=0;idxrd<2*IDXR-1;idxrd++){Dens[idxt][idxrd]=0.0;
Dens[idxt+IDXT][idxrd]=0.0;}*etot=0; // sum of single-particle energies
for(i=0;i<imax;i++){n=2*e[i].idxj+2;
if(fn<n){n=fn;}*etot+=n*e[i].E;
Fsptr=Fs[idxt][e[i].idxp][e[i].idxj][e[i].idxn];

74

Gsptr=Gs[idxt][e[i].idxp][e[i].idxj][e[i].idxn];
for(idxr=1;idxr<IDXR;idxr++){// debug "<=" --> "<"double Ft,Gt;
idxrd=2*idxr;
r=idxrd*dr;
Ft=Fsptr[idxr];
Gt=Gsptr[idxr];
Dens[idxt][idxrd]+=(Ft*Ft+Gt*Gt)*n/(4.0*M_PI*r*r);
Dens[idxt+IDXT][idxrd]+=(Ft*Ft-Gt*Gt)*n/(4.0*M_PI*r*r);}
Dens[idxt][0]=f6*Dens[idxt][2]+f7*Dens[idxt][4];Dens[idxt+IDXT][0]=f6*Dens[idxt+IDXT][2]+f7*Dens[idxt+IDXT][4];fn-=n;
if(fn<=0) break;}//>for(i)*etot-=NZ[idxt]*Mass[idxt];
if(fn>0){// changed from "if(fn<0)" on 07/9/3
fprintf(stderr,"error in the calculating density: insufficient number of levels\n");exit(1);}interpolate(Dens[idxt],2*IDXR-1);
interpolate(Dens[idxt+IDXT],2*IDXR-1);
fprintf(stderr,"%f %f %f %f #3#\n",Dens[0][0],Dens[1][0],Dens[2][0],Dens[3][0]);//return 0;}//------------------------------------------------------------------------------------int writeWavefunction(FILE *FO, int idxt, int idxp, int idxj, int idxn){ // writes the wavefunction to a stream
int i;
double r1, F, G, *Fsptr, *Gsptr;
Fsptr=Fs[idxt][idxp][idxj][idxn]; Gsptr=Gs[idxt][idxp][idxj][idxn];
F=0;G=0;
for(i=0;i<IDXR;i++){if(F<fabs(Fsptr[i]))F=fabs(Fsptr[i]);if(G<fabs(Gsptr[i]))G=fabs(Gsptr[i]);}fprintf(FO,"# (idxp,idxj,idxn)=(%d %d %d) E=%16.8f\n",idxp,idxj,idxn,Es[idxt][idxp][idxj][idxn]);
fprintf(FO,"# r(fm) F/%12e G/%12e Vve Vsc\n",F,G);
F=1.0/F; G=1.0/G;
for(i=0;i<IDXR;i++){r1=RadialGridSpacing*i;

75

fprintf(FO,"%8.4f %16.8e %16.8e %16.8e %16.8e\n",r1,Fsptr[i]*F,Gsptr[i]*G,Vv[idxt][2*i],Vs[idxt][2*i]);}}//------------------------------------------------------------------------------------int swap(orbital *e, int i, int j){ // called only by function "sort"orbital tmp;
tmp = e[i];e[i] = e[j];e[j] = tmp;}//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - int sort(orbital *e){
int idxt, i, j, imax=IDXP*IDXJ*NODE-1;
for(i=0;i<imax;i++){for(j=i+1;j<=imax;j++){if(e[i].E>e[j].E){
swap(e, i, j);}}}}//------------------------------------------------------------------------------------int initialPotential(int sw) {// 0th compoment of vector potential (MeV) and scalar potential (MeV)
int i,idxt;
double r,dr;
double A,R,rhov,rhos,rho3,t1,t2,Vv0[IDXT],Vs0[IDXT],Vc0; // local variable// external variables// r0_ws,a_ws,NZ[]
A=NZ[NEUTRON]+NZ[PROTON];
R=r0_ws*pow(A,1/3.0);
rhov=0.16;
rhos=rhov;
rho3=rhov*(NZ[NEUTRON]-NZ[PROTON])/A;
t1=pow(Coupling[OMEGA]*HbC/BosonMass[OMEGA],2)*rhov;
t2=pow(Coupling[RHO]*HbC/BosonMass[RHO],2)*rho3;
Vv0[NEUTRON]=t1+t2;
Vv0[PROTON ]=t1-t2;
Vs0[NEUTRON]=-pow(Coupling[SIGMA]*HbC/BosonMass[SIGMA],2)*rhos;
Vs0[PROTON ]=Vs0[NEUTRON];
Vc0=NZ[PROTON]*fineStructureConstant*HbC;
dr=RadialGridSpacing/2;
for(idxt=0;idxt<IDXT;idxt++){for(i=0;i<2*IDXR-1;i++){76

r=i*dr;Vv[idxt][i]=Vv0[idxt]/((1.0+exp((r-R)/a_ws))*(1.0+exp((-r-R)/a_ws)));Vs[idxt][i]=Vs0[idxt]/((1.0+exp((r-R)/a_ws))*(1.0+exp((-r-R)/a_ws)));if(idxt == PROTON){
if(r<=R){Vv[idxt][i]+=0.5*Vc0/R*(3-pow(r/R,2));}else {
Vv[idxt][i]+=Vc0/r;}}}
if(sw>0){
fprintf(stderr,"vector potential: t,strength,R,a=%d %f %f %f\n",idxt,VvStr[idxt],R_ws,a_ws);
fprintf(stderr,"scalar potential: t,strength,R,a=%d %f %f %f\n",idxt,VsStr[idxt],R_ws,a_ws);}}
return 0;}//------------------------------------------------------------------------------------// A set of functions to solve the Dirac eigenvalue equation with spherically symmetric// vector and scalar potentials which are regular at radius zero.
//* idx2qn//* setRadialGrid//* solveDirac//* rungeBCi//* rungeBCo//* (rungeFun1 and rungeFun2 :commented out)//
* runge//------------------------------------------------------------------------------------int idx2qn(int sw,int idxp,int idxj, double *J,int *w,int *L,int *kappa){// converts the indices into quantum numbers*J=idxj+0.5; // J = angular momentum
if(idxp == 0){*w=-1;} else {*w=1;} // w = quantum number "varpi"*L=idxj+(1+ *w)/2; // L=J+varpi/2 : orbital angular momentum*kappa=*w *(idxj+1); // kappa = varpi*(J+1/2)
if(sw > 0){
fprintf(stderr,"quantum numbers: J=%3d/2 varpi=%d L=%d kappa=%d\n",(int)(2*(*J)),*w,*L,*kappa);
return 0;}}//------------------------------------------------------------------------------------int setRadialGrid(int sw){ // sets up the parameters of the radial grid to express the radial wavefunctions
double dr;

77

Idxr.i=10; // For radial grid number <= Idxr.i, rungeBCi is used.
Idxr.o=2500; // For radial grid number >= Idxr.o, rungeBCo is used.
Rout=25.0; // (fm) corresponding to Idxr.oif(1 > Idxr.i || Idxr.i > Idxr.o || Idxr.o >= IDXR || Rout <= 0.0){
fprintf(stderr,"setRadialGrid: error: Idxr.i,Idxr.o,IDXR,Rout %d %d %d %f\n",Idxr.i,Idxr.o,IDXR,Rout);exit(1);}
dr=Rout/Idxr.o; RadialGridSpacing=dr;
Rin=Idxr.i*dr;
Idxr.m=(R_ws+3*a_ws)/dr+0.5;
if(Idxr.m < Idxr.i) Idxr.m = Idxr.i;
if(Idxr.m > Idxr.o) Idxr.m = Idxr.o;
Rmatch=Idxr.m*dr;if(sw>0){
printf("R(in,match,out)=%9.5f %9.5f %9.4f (fm)\n",Rin,Rmatch,Rout);
printf("idxr(in,match,out)=%9d %9d %9d dr=%15.10f (fm)\n",Idxr.i,Idxr.m,Idxr.o,dr);}}//------------------------------------------------------------------------------------double solveDirac(int idxp,int idxj,int idxn){ // calculates the energy and the wavefunction of the eigenstate
int i,idxt,node;
double maco,maco1,maco2; // matching condition is "maco=0"
double E,E1,E2,dE, Eprec=1.0e-14; // energy (MeV)
idxt = idxtGlobal;
runge(0,Mass[idxt]-0.001,idxp,idxj,idxn,&node,&maco);//

printf("solveDirac: runge done.\n"); // #1#

if(node < idxn || node == idxn && maco < 0){fprintf(stderr,"no bound states: For E=Mass, node=%d maco=%e\n",node,maco);
E=Ebig;goto fin;}
E1=Mass[idxt]-30.0; dE=10.0;
for(;;){
runge(0,E1,idxp,idxj,idxn,&node,&maco);////

printf("solveDirac: runge done in for(;;). E1=%f dE=%f node=%d maco=%f\n",E1,dE,node,maco); // #2#

if(node < idxn || node == idxn && maco < 0) break;78

E1=E1-dE; dE*=1.25;}
fprintf(stderr,"lower bound (E,node,maco)=(%f %d %e)\n",E1,node,maco);
E2=E1*0.8+Mass[idxt]*0.2;
for(;;){
runge(0,E2,idxp,idxj,idxn,&node,&maco);if(node > idxn || node == idxn && maco > 0) break;
E2=E2*0.8+Mass[idxt]*0.2;
if(Mass[idxt]-E2 < 0.001) { // Actually,this condition was alread checked.
fprintf(stderr,"no bound states(2): For E=Mass, node=%d maco=%e\n",node,maco);
E=Ebig;goto fin;}}
fprintf(stderr,"upper bound (E,node,maco)=(%f %d %e)\n",E2,node,maco);// bisection methodfor(i=0;i<60;i++){if(E2-E1<=Mass[idxt]*Eprec) break;
E=(E1+E2)*0.5;
runge(0,E,idxp,idxj,idxn,&node,&maco);// fprintf(stderr,"%12f %5d %15.12f %12f %12f\n",E,node,maco,E1,E2);
if(node < idxn || node == idxn && maco < 0){E1=E;} else {E2=E;}}
E=(E1+E2)*0.5;
runge(1,E,idxp,idxj,idxn,&node,&maco);Es[idxt][idxp][idxj][idxn]=E;
fprintf(stderr,"iter=%d E=%f\n",i,Es[idxt][idxp][idxj][idxn]);
fin:
Es[idxt][idxp][idxj][idxn]=E;
return E;}//------------------------------------------------------------------------------------int rungeBCi(double r, double E, int idxp, int L, double *F, double *G){// Boundary Condition for r -> 0
double vv,vs,epsilon,mu,rL,rL1,rL2,rL3;
int idxt;
idxt=idxtGlobal; // A global variable is copied to a local variable.
if(r <= 0.0){*F=0; *G=0;
if(r == 0) return 0; else return 1;}

79

vv=Vv[idxt][0]; vs=Vs[idxt][0];epsilon=(E-vv)/HbC;
mu=(Mass[idxt]+vs)/HbC;{//input : int L, double r, output : double rL=pow(r,L)
int pt=L; double xt=r, rt=1.0;
while(pt != 0){if(pt & 1) rt*=xt; xt*=xt; pt>>=1;} rL=rt;}
rL1=rL*r; rL2=rL1*r; rL3=rL2*r;*F=rL1-(epsilon*epsilon-mu*mu)/(4*L+6)*rL3;
if(idxp == 0){ // varpi=-1, L=J-1/2=0,1,2,...
*G=-(epsilon-mu)/(2*L+3)*rL2;}
else { // varpi=1, L=J+1/2=1,2,3,...
*G=(2*L+1)/(epsilon+mu)*rL-(epsilon-mu)/(2*L+3)*rL2;}
return 0;}//------------------------------------------------------------------------------------int rungeBCo(double r, double E, double *F, double *G){// Boundary Condition for r -> infinity
double Mass1;
Mass1=Mass[idxtGlobal];*F=exp(-r*sqrt(Mass1*Mass1-E*E)/HbC);*G=-sqrt((Mass1-E)/(Mass1+E))*exp(-r*sqrt(Mass1*Mass1-E*E)/HbC);
return 0;}//------------------------------------------------------------------------------------#define RUNGE_OPTION 3#if RUNGE_OPTION == 1// - - - - - - - - - - - - - - function version - - - - - - - - - - - - - inline double rungeFun1(double rinv,double F,double G,double E,int kappa,double V){
return (E-V+Mass[idxtGlobal])*G*(1.0/HbC)-kappa*F*rinv;}
inline double rungeFun2(double rinv,double F,double G,double E,int kappa,double V){
return -(E-V-Mass[idxtGlobal])*F*(1.0/HbC)+kappa*G*rinv;}#elif RUNGE_OPTION == 2// - - - - - - - - - - - - - - - macro version - - - - - - - - - - - - - - #define rungeFun1(rinv, F, G, E, kappa, V)\((-(kappa))*(rinv)*(F)+(1.0/HbC)*((E)-(V)+Mass[idxtGlobal])*(G))#define rungeFun2(rinv, F, G, E, kappa, V)\((1.0/HbC)*(-((E)-(V)-Mass[idxtGlobal]))*(F)+(kappa)*(rinv)*(G))#elif RUNGE_OPTION == 3// - - - - - - - - - - - - - - - inline version - - - - - - - - - - - - - - //Inline version uses neither functions nor macros.
#endif80

int runge(int sw, double E, int idxp, int idxj, int idxn, int *node, double*maco){// Runge-Kutta method for radial wave fn.
// sw & 0x01 > 0 --> the wavefunciton iwws stored in the global arrays// sw & 0x02 > 0 --> fprintf(stderr) some information/*
used global variables
Fs, Gs, Mass, idxtGlobal, RadialGridSpacingused functions
rungeBCi, rungeBCo, idx2qn*/
int i,is,ii,i2,i3;
int w,L,kappa;
int forbac; // 0 for forward solution, 1 for backward solution
int idxt;
double *Fsa,*Gsa;
double r1,r1inv;
double k1,l1,k2,l2,k3,l3,k4,l4,k,l;double F,F2,F3,F4,oldF, G,G2,G3,G4;
double dr, pmdr, halfpmdr; // pmdr = plus or minus drdouble J, vv, vs;double Fm[2], Gm[2];double c1,c2,c3,c4,c5,c6,c7;
idxt=idxtGlobal;
Fsa=Fs[idxt][idxp][idxj][idxn];
Gsa=Gs[idxt][idxp][idxj][idxn];
idx2qn(0,idxp,idxj, &J,&w,&L,&kappa);
dr=RadialGridSpacing;if(sw & 0x01){
for(i=0;i<=Idxr.i;i++){r1=i*dr;rungeBCi(r1,E,idxp,L,&Fsa[i],&Gsa[i]);}for(i=IDXR-1;i>=Idxr.o;i--){r1=i*dr;rungeBCo(r1,E,&Fsa[i],&Gsa[i]);}}*node=0;
for(forbac=0; forbac<2; forbac++){81

if(forbac == 0){
rungeBCi(Idxr.i*dr,E,idxp,L,&F,&G);
pmdr=dr; is=Idxr.i; ii=1;}
else {
rungeBCo(Idxr.o*dr,E,&F,&G);
pmdr=-dr; is=Idxr.o; ii=-1;}
c1=pmdr*kappa;
c2=pmdr*(1.0/HbC);
c3=E+Mass[idxt];
c4=E-Mass[idxt];
halfpmdr=pmdr*0.5;
r1=is*dr; r1inv=1.0/r1; c5=c1*r1inv;
i2=is*2; vs=Vs[idxt][i2]; vv=Vv[idxt][i2];for(i=is+ii ; ; i+=ii){i2=i*2;#if ( RUNGE_OPTION == 1 ) || ( RUNGE_OPTION == 2 )//----------------------- function and macro versions ----------------------k1=rungeFun1(r1inv,F,G,E,kappa,vv-vs)*pmdr;
l1=rungeFun2(r1inv,F,G,E,kappa,vv+vs)*pmdr;
r1=r1+halfpmdr; r1inv=1.0/r1;
i3=i2-ii; vs=Vs[idxt][i3]; vv=Vv[idxt][i3];F2=F+0.5*k1; G2=G+0.5*l1;
k2=rungeFun1(r1inv,F2,G2,E,kappa,vv-vs)*pmdr;
l2=rungeFun2(r1inv,F2,G2,E,kappa,vv+vs)*pmdr;
F3=F+0.5*k2; G3=G+0.5*l2;
k3=rungeFun1(r1inv,F3,G3,E,kappa,vv-vs)*pmdr;
l3=rungeFun2(r1inv,F3,G3,E,kappa,vv+vs)*pmdr;
r1=dr*i; r1inv=1.0/r1;
vs=Vs[idxt][i2]; vv=Vv[idxt][i2];F4=F+k3; G4=G+l3;
k4=rungeFun1(r1inv,F4,G4,E,kappa,vv-vs)*pmdr;
l4=rungeFun2(r1inv,F4,G4,E,kappa,vv+vs)*pmdr;#elif RUNGE_OPTION == 3//----------------------------- inline version ----------------------------k1=c2*(c3-vv+vs)*G-c5*F;
l1=c5*G-c2*(c4-vv-vs)*F;
r1=r1+halfpmdr; c5=c1/r1;
i3=i2-ii; vs=Vs[idxt][i3]; vv=Vv[idxt][i3];F2=F+0.5*k1; G2=G+0.5*l1;
c6=c2*(c3-vv+vs);c7=-c2*(c4-vv-vs);k2=c6*G2-c5*F2;
l2=c5*G2+c7*F2;
F3=F+0.5*k2; G3=G+0.5*l2;
k3=c6*G3-c5*F3;
l3=c5*G3+c7*F3;82

r1=dr*i; c5=c1/r1;
vs=Vs[idxt][i2]; vv=Vv[idxt][i2];F4=F+k3; G4=G+l3;
k4=c2*(c3-vv+vs)*G4-c5*F4;
l4=c5*G4-c2*(c4-vv-vs)*F4;#endifk=(k1+2*k2+2*k3+k4)*(1.0/6.0);l=(l1+2*l2+2*l3+l4)*(1.0/6.0);
oldF=F;
F=F+k;
G=G+l;
if(F*oldF<0) (*node)++;if(sw & 0x01){
Fsa[i]=F;Gsa[i]=G;}
if(i==Idxr.m) break;}Fm[forbac]=F; Gm[forbac]=G;} // end of for(forbac)*maco=Fm[0]*Fm[1]*(Fm[0]*Gm[1]-Gm[0]*Fm[1]);//fprintf(stderr,"%12f %5d %15e E,node,maco\n",E,*node,*maco);//---------------- normalization of the wavefunction : begin ---------------if(sw & 0x01){ double s,s1,s2,t1,t2,mf,nf1,nf2;
s1=0.0;for(i=0;i<Idxr.m;i++){t1=Fsa[i];t2=Gsa[i];s1+=t1*t1+t2*t2;}
s2=0.0;for(i=Idxr.m;i<IDXR;i++){t1=Fsa[i];t2=Gsa[i];s2+=t1*t1+t2*t2;}
mf=Fm[0]/Fm[1]; // It is preferable to choose between Fin/Fout and Gin/Gout// the one which suffers from less numeircal error.
s=dr*(s1+s2*mf*mf);
nf1=1/sqrt(s); // sign of the tail for r-> infinity agrees with rungeBCo
if(mf<0.0) nf1*=-1;
nf2=nf1*mf;
for(i=0;i<Idxr.m;i++){
Fsa[i]*=nf1;

83

Gsa[i]*=nf1;}
for(i=Idxr.m;i<IDXR;i++){
Fsa[i]*=nf2;
Gsa[i]*=nf2;}}//---------------- normalization of the wavefunction : end -----------------return 0;}// block 2 : begin//
This is an independently available part.
// =============================================================================//
proper/screened Poisson equation solver//
No global variables are referred.
// =============================================================================//------------------------- Poisson equation solver -------------------------int Poisson(int imax,double h,double *source,double *phi){double r0,r1;
double I1,I2;
double rhol;
int i,j,k,l,a1;
double intcoef[6]={ 11.0/1440.0, -93.0/1440.0, 802.0/1440.0,
802.0/1440.0, -93.0/1440.0, 11.0/1440.0};
//r=0
I2=0;
for(j=0;j<imax-1;j++){for(k=0;k<6;k++){l=j+k-2;r1=l*h;a1=abs(l);if(a1<imax) rhol=source[a1];else rhol=0;
I2+=intcoef[k]*rhol*r1;}}
phi[0]=h*I2;//r<0
for(i=1;i<imax;i++){r0=i*h;I1=0;for(j=0;j<i;j++){for(k=0;k<6;k++){l=j+k-2;84r1=l*h;a1=abs(l);if(a1<imax) rhol=source[a1];else rhol=0;
I1+=intcoef[k]*rhol*r1*r1/r0;
}
}
I2=0;
for(k=0;k<6;k++){for(j=i;j<imax-1;j++){l=j+k-2;r1=l*h;a1=abs(l);if(a1<imax) rhol=source[a1];else rhol=0;
I2+=intcoef[k]*rhol*r1;}}
phi[i]=h*(I1+I2);}
return 0;}//------------------ screened Poisson equation solver ------------------------int screenedPoisson(int imax, double h, double m, double *source, double *phi){double D;// the value of Green’s function
double r0,r1;// variables r and r’ of Green’s function
double s,s1,s2,exp1,exp2,sourcel;
int i,j,k,l,n,al,nexp;
double intcoef[6]={ 11.0/1440.0, -93.0/1440.0, 802.0/1440.0,
802.0/1440.0, -93.0/1440.0, 11.0/1440.0};#define NEXP 100000
double expmmr[NEXP];
int expsft = 3;
nexp=NEXP;
exp1=exp(-m*h);for(k=0;k<NEXP;k++){if(k % 20 == 0) expmmr[k]=exp(-m*h*(k-expsft));
else expmmr[k]=expmmr[k-1]*exp1;
if(expmmr[k]<1.0e-32) {nexp=k+1;
break;}}
if(expmmr[nexp-1]>1.0e-16) {fprintf(stderr,"screenedPoisson: warning: You have to increase NEXP\n");fprintf(stderr,"expmmr[%d]=%e\n",NEXP,expmmr[NEXP-1]);}#undef NEXP

85

// calculation of phi(r) at r=0
s2=0;
for(j=0;j<imax-1;j++){ // integral over [0,(imax-1)*h]for(k=0;k<6;k++){l=j+k-2;r1=l*h;n=l+expsft; if(n<nexp) exp1=expmmr[n]; else exp1=0;
al=abs(l); if(al<imax) sourcel=source[al]; else sourcel=0;
s2+=intcoef[k]*sourcel*r1*exp1;}}
phi[0]=h*s2;//

fprintf(stderr,"phi[0]=%.7f\n",phi[0]);// calculation of phi(r) for r>0
for(i=1;i<imax;i++){r0=i*h;s1=0;for(j=0;j<i;j++){ // integral over [j*h,(j+1)*h], (j+1)*h<=r0
for(k=0;k<6;k++){l=j+k-2;r1=l*h;n=i-l+expsft; if(n<nexp) exp1=expmmr[n]; else exp1=0;
n=i+l+expsft; if(n<nexp) exp2=expmmr[n]; else exp2=0;
al=abs(l); if(al<imax) sourcel=source[al]; else sourcel=0;
s1+=intcoef[k]*sourcel*r1*(exp1-exp2);
}
}
s2=0;
for(j=i;j<imax-1;j++){ // integral over [j*h,(j+1)*h], j*h>=r0
for(k=0;k<6;k++){l=j+k-2;r1=l*h;n=l-i+expsft; if(n<nexp) exp1=expmmr[n]; else exp1=0;n=l+i+expsft; if(n<nexp) exp2=expmmr[n]; else exp2=0;
al=abs(l); if(al<imax) sourcel=source[al]; else sourcel=0;
s2+=intcoef[k]*sourcel*r1*(exp1-exp2);}}
phi[i]=h*(s1+s2)/(2*m*r0);}
return 0;}// block 2 : end// block 3 : start//
No global variables are referred.

86

//----------------------------------------------------------------------------int interpolate(double *y, int n){/*double y[n],y[0],y[2],y[4],...,y[n-5],y[n-3],y[n-1] : given.
y[1],y[3],y[5],...,y[n-6],y[n-4],y[n-2] : to be interpolated.
n must be an odd integer not less than 7.
y[-i]=y[i] assumed.
y[i]=0 assumed for i >= n.
*/
int i;
double f0= 3.0/256.0;
double f1=-25.0/256.0;
double f2= 75.0/128.0;
double f3= 75.0/128.0;
double f4=-25.0/256.0;
double f5= 3.0/256.0;
if(n%2==0||n<7){exit(1);
}
y[1]
y[3]

=f0*y[4]=f0*y[2]+f1*y[2]+f1*y[0]+f2*y[0]+f2*y[2]+f3*y[2]+f3*y[4]+f4*y[4]+f4*y[6]+f5*y[6]+f5*y[8];;for(i=5;i<=n-6;i+=2){
y[i]=f0*y[i-5]+f1*y[i-3]+f2*y[i-1]+f3*y[i+1]+f4*y[i+3]+f5*y[i+5];}// NB) density beyond the maximum radius is approximated with zero.
// A better treatment is to construct 5-pt and 4-pt formulae// and use them for these points.
y[n-4]=f0*y[n-9]+f1*y[n-7]+f2*y[n-5]+f3*y[n-3]+f4*y[n-1];y[n-2]=f0*y[n-7]+f1*y[n-5]+f2*y[n-3]+f3*y[n-1];}// block 3 : end

87