原子核の相対論的平均場模型

2007 年 2 月
福井大学 工学部 物理工学科

15 年度入学 48 番 山田昌平

目次
序章

2

第 1 章 原子核の成り立ち

1.1

原子核

1.2
1.3

原子核の密度分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

アイソスピン . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 2 章 Dirac 方程式

4
4
4
5
6

2.1
2.2

Dirac 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
中心力ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
7

2.3
2.4

固有値問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 中間子

3.1

中間子場

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10
10

3.2
3.3

Green 関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Gauss-Seidel 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10
12

第 4 章 計算結果

4.1
4.2

13

核子の波動関数の例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
中間子場. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

中間子の古典場の Green を利用した解

4.2.1
4.2.2
4.2.3

. . . . . . . . . . . . . . . . . .
中間子の古典場の反復解法 . . . . . . . . . . . . . . . . . . . . . . . . .
動径グリッド の間隔を変えて比較する . . . . . . . . . . . . . . . . . . .

第 5 章 まとめ

13
14
15
15
17
19

関連図書

20

謝辞

21

付録

Program List

22

1

序章
今日, 私たちを取り巻く世界は原子から出来ていることは知られている. 原子は物質の最小構成単位の粒子である。現在発見されているものだけでも約 3000 種類存在し 、数え方によっては約 6000 種類に達する。しかしながら、電子の数が等しいものを同じ原子と考えると、約 110 種類の元素にまとまる。その原子の構造が次第に明らかになり, 原子核の存在も物理学の中で大きな役割りをもつことがわ
かった.
現在の最先端の知識の一つである原子核物理学は , 素粒子や宇宙, 物性など幅広くい分野の土台となる学問であり, また原子力やアイソト ープ生産のような近代的な応用分野の基礎となるものでもある.
本論文では相対論的平均場模型を用いて, 原子核の特性について研究を進める.
原子核は陽子と中性子から構成されてれており, 粒子同士が中間子を交換することによって引力が働く, これを核力という. 核力は非常に近い距離において陽子間のクーロン力を上回るために , 原子核が束縛される. ちなみに原子核の半径はウラン 238 では半径 R は約 7.0 fm になる.
次に原子核のエネルギーは大きいので , 核子の速度は光速に近く相対論的運動学は一見重要と思える.
しかし , 原子核の基底状態付近では核子の運動エネルギーは平均すると約 22 MeV であり, 原子核から１個の核子を分離するのに必要な分離エネルギーは 5∼10 MeV ある. これらは核子の静止エネルギー
約 940MeV と比べて十分小さいので , 相対論的運動学はあまり重要ではないと考えられる.
しかし , 相対論を用いることにより大きな利点が二つある. まず非相対論ではスピン軌道力は後から付
け加えなければならないが , 相対論では自然に導かれる. 第二に, 原子核の状態が高温 · 高圧, つまり運動エネルギーがきわめて高い時にも外挿できるのも利点である.
非相対論と相対論とを比較すると, 非相対論では普通, 核子同士が直接相互作用する. その様子は Schrö

dinger 方程式で現される. 一方, 相対論では Dirac 方程式で現される核子が中間子を放出 · 吸収する. この様子は場の理論で扱うが , その解を求めることは非常に困難なので , 以下の二つの仮定を置いて平均場模型という枠組に置き換える. 第一の仮定として, 中間子を古典場に置き換える. 第二に、Dirac 方程式
から負のエネルギーの解が導き出される. 普通, 励起状態にあるとき光子を放出してエネルギーを失い ,低い励起状態に落ちていき , 最後に基底状態に落ち着く. ところが負のエネルギーの解があると , 基底状
態はエネルギーが最低の状態でなくなる. 負のエネルギー状態のエネルギーには下限がないので限りな
く光子を出し続けることになる. 原子の基底状態は安定なので , 基底状態からの遷移が起こってはならな
い.

1930 年に Dirac は ,「負エネルギーの状態がすべて占有されていれば , パウリ原理のために負エ

ネルギーの状態への遷移は禁止される」ことに注意を払った.Dirac によれば , 粒子が存在しないと考えられる真空は , すべての負エネルギーの状態が核子によって完全に占有されている状態なのである. しかし 負のエネルギーの核子の寄与を無視する. これを no-sea 近似と呼ぶ.
ラグランジアン密度から出発し ,Dirac 方程式の固有値方程式が求まる. その方程式の波動関数を動径波動関数と球面調和関数に分離すると、原子核の動径関数の従う二成分一階の方程式となって表される.
この方程式の中に二つのポテンシャルがあり, 二つの値を代入することにより, 波動関数が出てくる.
そして Green 関数を用いて計算する方法と ,Gauss-Seidel 法を用いて計算する, 二つの方法で核子密度から中間子場を求めることを試みる。
最後に自己無撞着に波動関数, 核子密度, 中間子場のそれぞれの解を求めることができるので、波動関数から核子密度を求め、次に核子密度から中間子場を求めて、最後に中間子場から波動関数を求める。
この三つのサイクルを変化が無くなるまで繰り返すと、波動関数, 核子密度, 中間子場のそれぞれの解を2

求めることができる。

3

第1章

1.1

原子核の成り立ち

原子核

原子核は陽子と中性子からつくられている. 陽子と中性子を総称して核子と呼ぶ. 原子核から核子がとび出さないのは , 核子の間に引力が働くからであろう. この力を核力と呼ぶ. 核力はクーロン力ではな
い . 中性子にはクーロン力は働かないし , 陽子の間のクーロン力は斥力である. 核力は万有引力でもない .
陽子の間の万有引力はクーロン力の約 1/1036 にすぎないからである.
原子核の内部での核子の間のような非常に短い距離にでは , 核力は陽子の間のクーロン力よりはるかに強い力でなければならない . この強い力を強い相互作用とも呼ぶ. 陽子と陽子の衝突は , 陽子間の距離が約 2 × 10−15 m 以上ではクーロン斥力による散乱として説明されるので , 核力はその到達距離が約

2 × 10−15 m という非常に短い距離の力である.

核力の原因を研究し , 核力は質量が電子の質量の約200 倍の新しい粒子を核子同士が交換されるために生じる, という新しい理論を 1935 年に提唱したのは湯川秀樹である. この新しい粒子を今日は中間子 (meson) と呼ぶ。
a という種類の質量 ma の中間子が二つの核子の間に交換されると生じる核力の到達距離は , 量子力→
学の不確定性関係から , 簡単に次のように理解される. 図 2-1 のように,−
r1 にある核子 1 が中間子を放出→−
し , r2 にある核子 2 がそれを吸収する過程を考える. 中間子の放出に伴うエネルギーの不確定の程度は
ma c2 で , 対応する時間の程度は ∆t ∼ h̄/ma c2 である. 中間子の伝播の速さは光速 c 以下だから , 中間子が核子 1 から離れる距離は c∆t ∼ h̄/ma c ∼ 0.379 fm 以下である. その範囲に核子 2 があれば , 中
間子を吸収して, 核子が生じる.

p

n
meson

Dirac eq.
pn

図 1.1: 2 核子間で中間子が交換されるダ イアグラム

1.2

原子核の密度分布

電子や陽子と原子核との散乱実験によって, 原子核内の核子の密度分布 ρ(r) が調べられる.r を原子核の重心から測った距離である. 電子散乱では電磁相互作用により電荷分布 ρc (r) が , 陽子散乱は強い相互作用によって陽子と中性子の寄与の和である密度分布が調べられる. 原子核の体積は構成核子の個数 Aに比例すること, 中心部の密度 ρ(0) は A によらずほぼ一定値 ρ(0) ' 0.17 fm−3 であることが知られ4ている. ちなみに 1 fm = 10−13 cm である。これを密度の飽和性という。
核子密度 ρ(0) の 90 ％から 10 ％に減少する距離と定義される表面の厚さは 2∼3 fm であること等が明らかになった.A が小さい軽い核を除けば , 表面の厚さが原子核全体の広がりに比べて小さいので ,ρ(r)が ρ(0) の 1/2 となる距離を核半径 R と定義すると ,1

R = r0 A 3 (r0 = 1.2 f m)

(1.1)という関係がよい近似で成り立っている.

1.3

アイソスピン

陽子も中性子もスピン 1/2 で、質量はほぼ等しいので、核子の二つの状態を考える。スピン J = 1/2の粒子の電子にはスピン上向き Jz = 1/2 と下向き Jz = −1/2 の状態があるように、核子はアイソス
ピン I = 1/2 の粒子で、陽子はアイソスピンが上向き I3 = 1/2、中性子は下向き I3 = −1/2 の状態だと考える。
アイソスピンは仮想的な荷電空間での角運動量で、核力は荷電空間での回転でほぼ不変である。これ
を核力の荷電独立性という。この核力の荷電独立性のために、核反応ではアイソスピンが保存する。また、スピン J の原子・分子のエネルギー準位は Jz = J, J − 1, · · · , −J − 1 の合計 2J + 1 個の状態が縮退しているように、アイソスピンが I の準位は I3 = I, I − 1, · · · , −I − 1 の合計 2I + 1 個の状態が縮退している。この質量数とスピンが等し く、電荷が異なる 2I + 1 個の原子核の準位の組を荷電
多重項または荷電 2I + 1 重項という。

5

第 2 章 Dirac 方程式

2.1

Dirac 方程式

質量 m の自由粒子に対する Schrödinger 方程式

ih̄

∂ψ
∂t

=−

h̄2
2m

∇2 ψ

(2.1)

は非相対論的な波動方程式である. 一方, 特殊相対論に於いて, 光速 c に近い速度で動く粒子のエネルギー

E と運動量 p の間には ,E 2 = (pc)2 + (mc2 )2(2.2)

が成り立つ. ここで m は静止質量である. この式は Einstein の関係式という. この方程式を (2.1) 式に代入すれば ,(∇2 −∂2

)φ = µ2 φ , µ =

mc

(2.3)
h̄
となる. これを Klein-Gordon 方程式という. この方程式は相対論的な波動方程式だが , スピンが入っていない . 速度が遅いときは Schrödinger 方程式と近似的に一致するスピン 1/2 の粒子の相対論的な波動方程式を導く. そこで Einstein の関係式を,E=

c2 ∂t2p2 c2 + m2 c4
E

2

mc
→
−
=−
v ·→
p +(
mc2 )
E

(2.4)

と書き直して
2

∂
mc
−
→
→
→
p → −ih̄ , −v →−α c , E → ih̄
,
→β
∂t
E

(2.5)

と置き換えると , 新しい波動方程式

ih̄

∂ψ
∂t

+ ih̄cα · ∇ψ − βmc2 ψ = 0(2.6)

が得られる. この方程式はイギリスの Dirac によって導かれたので ,Dirac 方程式という.
ここで α と β は関係はαk β + βαk = 0, β 2 = 1, αk αj + αj αk = 2δkj

(2.7)

を満たす行列だとすると , 式 (2.6) に微分演算子

ih̄

∂ψ
∂t

− ih̄cα · ∇ + βmc2(2.8)

を作用させると ,Ψ は Klein-Gordon 方程式を満たす,

h̄2

∂2ψ
∂2t

+ h̄2 c2 ∇2 ψ − (mc2 )2 ψ = 0

(2.9)

→
関係 (2.8) を満たす −α と β は少なくとも 4 行 4 列の行列でなければならない . この理由としては 2 行
2 列や 3 行 3 列では (9) を満たすものがないからである. たとえば ,
[
]
[
]
0 σk
I 0
αk =
,β =
σk 0
0 I
6

(2.10)

は関係 (9) を満たす. このときの σk はパウリ行列で ,I は 2 行 2 列の単位行列である,
[
]
[
[
[
]
]
]
0 1
0 −i
1 0
1 0
σ1 =
, σ2 =, σ3 =
,I =
1 0
i 0
0 −1
0 1
このように,α,β は 4 行 4 列の行列式なので ,Ψ は 4 成分を持ち,
Ψ1 (x, t) Ψ2 (x, t) Ψ(x, t) =  Ψ3 (x, t) 
Ψ4 (x, t)(2.11)(2.12)のようになる.

2.2

中心力ポテンシャル

中心力による静ポテンシャル V (r) の中の粒子を記述する Dirac 方程式は

ih̄

∂ψ
∂t

→
+ i−
α · ∇ − βm − V (r)ψ = 0(2.13)と書ける. この場合の Dirac のハミルト ニアンはHD = α · p + βm + V (r).

(2.14)

である. これは回転と反転によって不変である. すなわち→−[HD , J ] = 0, [HD , P ] = 0.

(2.15)

それゆえ決まった角運動量と偶奇性をもつ固有値 Ψ を求めることができる.

HD Ψ = EΨ.
このときの解 Ψ を次の形に書くと都合がよい

Ψ=

(

φ
χ

(2.16)).

(2.17)→
ただし ,φ は Ψ の 1,2 成分,χ は Ψ の 3,4 成分である.φ と χ は −r および z 軸方向のスピン成分 µ の関
数である.
−
→
Ψ は J 2 , Jz , P に共通な関数である. 角運動量を固定する関数を (J M ) で表す. 空間反転に結びついた演算子 P に対する偶奇性については , 量子数 $(バルパイ) を{1(偶奇性 (−)J +1/2 の状態に対して)$=0(偶奇性 (−)J −1/2の状態に対して)(2.18)

Mと定め, 全角運動量 (J M ) の固有関数を YLJ(θ, ϕ) とする. これは軌道偶奇性の演算子 P (0) に対する偶奇性 (−)L をもつ関数である. 角運動量の結合規則により,L は二つの値, すなわち

L=l≡J+

1
2$, L = l0 ≡ J −1
2$(2.19)

M
しか取りえない . そして二つの関数 YlJと YlM
0 J は反対の偶奇性をもつ .
M
要約すれば ,Y$J
によって角運動量 (JM ) と偶奇性 (−)$+1/2 をもつ状態を表わすとき,()
M
1
F (r)ylJM
Ψ$J =riG(r)ylM0Jという形に書き表わすことができる. ここで F と G は r の任意の関数である.

7(2.20)

2.3

固有値問題

固有値問題
M
HD ΨM
$J = EΨ$J(2.21)を解いてみる.
動径運動量 pr と動径速度 αr をつぎのように導入する.

pr ≡ −i

1 ∂
r ∂r

r

(2.22)→
αr ≡ (−
α · r̂) = ρ1

−
→
→
σ ·−
r

(2.23)r

そうすると[αr (pr −i$(J + (1/2))rという形に書き表わされる. このとき

Mβ) + mβ + V (r)]ΨM$J = EΨ$J

M→(−
σ · r̂)ylJ= −ylM
0J
M
M→−( σ · r̂)y 0 = −y(2.25)(2.26)lJ

l J(2.24)

を使えば , 上の方程式から動径関数 F (r) と G(r) に関する二つの微分方程式が一組に成って導き出される. すなわち

[
−
[

d
dr
d
dr

+

$(J + 12 )+$(J + 12 )rr]
G(r) = (E − m − V )F (r)(2.27)

F (r) = (E + m + V )G(r)(2.28)]と書ける.
またラグランジアン密度から出発して, 原子核の連立一階微分動径方程式を求めると[]$(J + 21 )d−+G(r) = (E + V0 − m − Vs )F (r)drr][$(J + 21 )d+
F (r) = (E + V0 + m + Vs )G(r)drr(2.29)(2.30)と二成分一階の微分方程式となって現される。この連立方程式の F (r) または G(r) を消去し 、一つの式に書き直すと Schrödinger 方程式に似た式が導き出される。
ここで V0 はベクト ルポテンシャル ,Vs はスカラーポテンシャルである. ベクト ルポテンシャル , スカ
ラーポテンシャルはVs (r) = gσ φσ (r)(2.31)V0 (r) = gω φ0ω (r) + gρ τ3 φ0ρ3 (r) + e

1 − τ3A0 (r)(2.32)
2
φ は中間子場を表わし , その上付きの文字は各中間子のスピン , 下付きの文字は中間子のアイソスピンを3表す.e 1−τA0 (r) はクーロン場である. 上の式では用いられている σ 中間子,ω 中間子,ρ 中間子につい2てまとめたものが表１である.
表１. 中間子の種類
アイソスピン t 質量 (MeV/c2 )

中間子

スピン s

σ

0

0

520

φσ

ω

1

0

783

φµ
ω

ρ

1

1

770

φµ
ρ 3

8

場

2.4

Runge-Kutta 法

一段階法の特徴を生かしつつ,x0 x1 の間の, 中間の x における解の導関数, すなわち f (x, y) の関数
値を何回か計算し , それらの重みつき平均を用いることによって, 次数を達成するのが ,Runge-Kutta
法の原理である.
T解を求める範囲の上限値を T とし , 分割する個数 N (N は自然数) を定めて ∆ = Nとする. 初期条件として Y0 = a を用いる.tj の情報から Runge-Kutta 法で Yj+1 を求める方法は , 以下のとおりである.

k1 = f (tj , Yj )∆t∆t
, Yj +
k1 )k2 = f (tj +
2
2
∆t
∆t
k3 = f (tj +, Yj +
k2 )
2
2
∆t
∆t
k4 = f (tj +, Yj +k3 )1
1∆t
Yj+1 = Yj +(k1 + 2k2 + 2k3 + k4 )6
このようにして, 格子点 j の情報から次点の J + 1 を導き出す. この式の内容を説明する.

(2.33)(2.34)(2.35)(2.36)(2.37)
k1 は格子点 j での微分解 (傾き) である. この情報を使い k2 を求める.

k2 も微分解の式であるが , 代入される t の成分は tj + ∆tである. これは , 格子点 j と j + 1 との中間
2点での差分解を求めることに他ならない . 代入される Y 成分は Yj + ∆tk であるから , こらは t 成分の2 1に合うように k1 の値を使って Y の値を出したものである. k3 は tj + ∆tの Y 成分を出す
tj + ∆t2
2
のに k2 を使う. k4 も tj + ∆t(=tj+1 ) で Y 成分を出すのに k3 を使い , そこでの微分解を求める. この様にしてもとめて k は各地点での微分解であるから , これを重みを付けて平均をとる.(4.5) 格子点の端である k2 , k3 は重みを 2 として計算する.
この様にして,Yj での差分解を求める. さらにこの差分解から求めた Yj+1 を Yj として考え, 次点のYj+1 を求めていく. これが Runge-Kutta 法のアルゴリズムである.

9

第3章

3.1

中間子

中間子場

中間子と光子の古典場の従う方程式は(−∆ − m2σ )φσ = −gσ ρs (r) − g2 σ 2 (r) − g3 σ 3 (r),(3.1)(−∆ − m2ω )φ0ω = −gω ρv (r),(3.2)(−∆ − m2ρ )φ0ρ 3 = −gρ ρ3 (r),(3.3)−∆A0 (r) = eρp (r),(3.4)

である。それぞれ σ 中間子,ω 中間子,ρ 中間子と光子の従う方程式である.mσ ,mω ,mρ は σ 中間子,ω 中間子,ρ 中間子の質量、A0 はクーロンポテンシャルである。gσ ,gω ,gρ は結合定数であり、ρs (r),ρv (r),ρ3 (r),ρp (r)は各種の密度である. 次に核子密度は

ρs =
ρv =
ρ3 =

1 ∑
L3

Ψ†i γ0 Ψi ,1 ∑L3

Ψ†i Ψi ,(3.6)Ψ†i τ3 Ψi ,

(3.7)

i

1 ∑
L3

(3.5)

i

i

と表す。Ψ は Dirac 方程式の動径波動関数の lower component であり、τ3 は荷電空間のベクト ルの第三成分である。L は中間子を箱の中に入れて考えた時の箱の一辺である。

3.2

Green 関数

たとえば , 中間子の古典場の従う方程式

d2
dr 2

φ(r) +

2 dr drφ(r) − m2 φ(r) = −gρ(r)を解くために、まず次の方程式を考える。
)
( 2
2 d
d
2
+
− m D(r, r 0 ) = δ(r − r 0 )
dr 2
r dr

(3.8)

(3.9)

この解は中間子の source term、と Green 関数を用いて書き現すことができる. その Green 関数の具体形がD(r, r 0 ) = −

1
mrr 0

sinh(mr< ) exp(−mr> ),r< = min(r, r 0 ),r> = max(r, r 0 ).
10

(3.10)

この Green 関数の右辺を 0 に置き換えた方程式で、原点で正則な解は 1r sinh(mr) の式になり、無限−mr

遠で正則な解は , e rの式になるような境界条件を持つ. 図 3.1 は縦軸 −Drr 0 , 横軸 r の Green 関数

を現したものである。図 3.1 から r = r 0 の時, 一階微分が不連続になり, 二階微分すると δ 関数が出て
くる.

0.6

sinh(mr<)exp(-mr>)

0.5

-Drr’

0.4
0.3
0.2
0.1
0
0

r’

1

2

3

4

5

r
図 3.1: Green 関数

もし 、方程式 (3.9) の解 D(r, r 0 ) が得られたならば 、一般の場合の ρ(r) に対して次の計算をすると∫ ∞φ(r) =
r 02 dr 0 [−gρ(r 0 )]D(r, r 0 ).
(3.11)
02

d2 d2
ここで f (r) = −gρ(r),£ = dr
2 + r dr − m とおいて、得られた関数 φ(r) は∫ ∞{£D(r, r 0 )} f (r 0 )r 02 dr 0 .
I(r) =(3.12)

0

I(r) = f (r) なので、£D(r, r 0 ) = 0 とは考えずに
∫ r−ε
∫ r+ε
∫ ∞
0
0
02
0
0
0
02
0
I(r) ={£D(r, r )} f (r )r dr +{£D(r, r )} f (r )r dr +{£D(r, r 0 )} f (r 0 )r 02 dr 0
0

r−ε

∫ r+ε

r+ε

{£D(r, r 0 )} f (r 0 )r 02 dr 0

=
r−ε

∫ r+ε
I(r) = lim

ε→+0

r−ε
∫ r+ε

= lim

ε→+0

r−ε

{£D(r, r 0 )} f (r 0 )r 02 dr 0
f (r 0 )r 02

∂2
∂r 2

D(r, r 0 )dr 0

11∫ r+ε= f (r)r 2 lim

ε→+0

[

= f (r)r

2

lim

ε→+0

r−ε

∂
∂r

∂

∂

∂r 0 ∂r
0

D(r, r 0 )dr 0]r−εD(r, r )r+ε式の導出が不完全なので書き直すとI(r) = f (r)(3.13)となって、求める方程式の解になっている。このような関数 D(r, r 0 ) を Green 関数という。

3.3

Gauss-Seidel 法

方程式 (3.9) を解くことは連立一次方程式を解くことに帰着される。その方程式は

d
1
φi =
(−φi+2 + 8φi+1 − 8φi−1 + φi−2 ) + O(h4 )
dr
12h
d2
1
φi =
(−φi+2 + 16φi+1 − 30φi + 16φi−1 + φi−2 ) + O(h4 )
2
dr
12h2

(3.14)(3.15)

方程式 (3.14),(3.15) を (3.9) に代入し 、両辺に h2 を掛け、φ̃i = φ(i, φi+2 , φi+1 , φi−1 , φi−2 , ρi )とおくと、

φ̃ = c2 (1 +2i)φi+2 + c1 (2 +2i)φi+1 + c1 (2 −2i)φi−1 + c2 (1 −

1
2
− 12
3
,
c
=
,
c
=
1
2
5
5
+ m2 h2+ m2 h2+ m2 h2222c0 = 5h22i)φi−2 + c0 ρi + O(h5(3.16))(3.17)

連立１次方程式 Ax = b を解く一方法である。係数行列 A を対角成分が 0 の行列 A0 と対角行列 Dとの和に分解すれば , 方程式は (A0 + D)x = b, すなわち x = D 1 (b − A0 x) と書ける。したがって, 置
き換えれば x ← D 1 (b − A0 x) を繰り返し , 収束すればその x が解である。いつも収束するわけではな
いが、行列 A が正定置 (固有値がすべて正) であれば収束する。特に、対角成分が非対角成分より大きい行列なら収束する。一般には Gauss 法が速いことが多いが、この方法は係数の書き換えがないので、成分の多くが 0 であるような行列 (疎行列) で便利である。Gauss ではせっかくの疎行列でなくなってしまう.

Gauss-Seidel 法の置換え
x ← x + (D 1 (b − A0 x) − x)(3.18)= (1 − c)x + cD 1 (b − A0 x)(3.19)の c = 1 の場合と考えられるが , この c を 1 より少し 大きくすると収束が速くなることがある. このような方法を SOR(successive over-relaxation) 法という.
この Gauss-Seidel 法を用いて

φ0i = φi + (φ̃ − φi )(3.20)= (1 − c)φi + cφ̃

(3.21)

この方程式を用いて反復代入ループを使い中間子場を求める。

12

第4章

4.1

計算結果

核子の波動関数の例

鉛 (A=208) の原子核での場合、各ポテンシャルは図 4.1 のようにとると適当である。これらのポテンシャルを方程式 (2.29),(2.30) に代入すると F (r) と G(r) の従う微分方程式が求まる。波動関数を求めた軌道の量子数を書け図 4.2 を見ると、原点付近では F (r) と G(r) の傾きに違いがわかる。これは Dirac 方程式の特徴として、軌道角運動量 l, l0 の二つある。F (r) は s 軌道、G(r) は p 軌道の違いからである。また振動しているのもわかり、0 に近づくところでは指数関数的な特徴も見てとれる。

potential energy
300

vecter potential
scalar potential

Energy(MeV)

200
100
0
-100
-200
-300
-400
0

5

10
radius(fm)

図 4.1: ポテンシャル V0 ,Vs

13

15

20

wave function
1

F(r)
G(r)

0.5
0
-0.5
-1
0

5

10

15

20

図 4.2: 動径波動関数 F (r),G(r)

4.2

中間子場

核子密度から中間子場を計算機で解いた結果を記す。解いた原子核は鉛 208 である。鉛 208 は二重閉 k 殻なので球形であるなど開発したプログラムに適した原子核である。
グリッド 上でしか、ρ(r) の値が与えられないとうい問題点に対し 、中間子の古典場の Green 関数を利用いた解と ,Gauss-Seidel 法を用いた中間子の古典場の反復解法から計算し , 二つの方法で中間子場を導き出す。
鉛 208 の核子密度を次のように仮定した。

φ(r) =

1(4.1)(1 + exp r−R)(1 + exp −r−R)aa

この式は偶関数化のため通常の Woods-Saxon 関数に修正を加えたものである。ここで R = 7.0 fm

density(fm^-3)

であり、a = 0.7 fm である (図 4.1)

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

rho

0

5

10

15
r(fm)

図 4.3: Green 関数

14

20

25

30

4.2.1

中間子の古典場の Green を利用した解

まず、中間子の古典場の Green 関数を利用いた解についての計算法を考えてみる。π 中間子の質量140 MeV に等し くとる。動径グリッド の間隔を 0.01 fm とした。
方程式 (3.8) の解は、(3.10),(3.11) となり、φ(0) の表式を求めると∫ ∞0φ(0) =
e−mr ρ(r 0 )r 0 dr 0

(4.2)

0

となる。ここでは台形則による積分をする。直線で補間するため台形則は精度が低い。但し 、滑らかな関数の一周期にわたる積分で、また r = ±∞ で滑らかにゼロに収束する関数の積分の場合では最高精∫b度になる。積分 S = a f (x)dx を台形則で近似すると

S = h(f0
2+ f1 + · · · +fn
2) + O(h)(4.3)となる。また Simpson 則は二次関数で補間する近似法があるが、この近似で先程の S を表すと

S=

h
3

(f0 + 4f1 + 2f2 + 4f3 + · · · + fn ) + O(h4 )

(4.4)

となる。
方程式 (4.2) の r = 0 での値では彼積分関数が原点について対称で遠方で速やかにゼロに収束するため、台形則による積分で高い精度が得られる。一方、r > 0 での値では彼積分関数が r = r 0 で一階以上の微分が不連続となるため、台形則による積分結果は精度が低く、h の１乗に比例する。そのように積分して中間子場 φ(r) が求まる。

4.2.2

中間子の古典場の反復解法

次に中間子の古典場の反復解法についての計算法について考える。同じように π 中間子の質量 140MeVに等し くとる。グリッド の上の点しか分からないので、動径グリッド の間隔を 0.01 fm とした。
まず核子密度に比例する中間子場の湧き出し ρ(r) の値と中間子場 φ(r) の値を同じとして Gauss-

Seidel 法の反復のための初期値の設定する。そして Gauss-Seidel 法の減衰因子 c = 1.99 とし 、反復
代入ループで収束するまで繰り返す。その時、更新後の中間子場は方程式 (3.21) となり、収束すれば中間子場 φ(r) が求まる。この Gauss-Seidel 法を用いる場合は、不連続点に煩わされないのが長所であるが、反復解法なので、収束の判定がわかりにくいこと、収束が非常に遅くなることが短所である。

15

図 4.4 は二通りの計算法で導いた中間子場である。rho は核子密度、pho(GF) は中間子の古典場の Green 関数を利用した解から求めた中間子場であり、pho(GS) は中間子の古典場の反復解法から求めた中間子場である。図 4.5 の縦軸は常用対数表示である。
2rho
phi(GF)
phi(GS)

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0

5

10

15

20

25

30

r (fm)

図 4.4: 核子密度と中間子場

100

rho
phi(GF)
phi(GS)

1

0.01

1e-04

1e-06

1e-08

1e-10

1e-12
0

5

10

15

20

25

r (fm)

図 4.5: logscale y の核子密度と中間子場

1630

4.2.3動径グリッド の間隔を変えて比較する

図 4.6, 図 4.7 は二つの計算法で動径グリッド の間隔を h=0.01,0.5,0.1 fm と変えてあらわしたものである。
2

h=0.01
h=0.5
h=0.1

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0

2

4

6

8

10

12

14

r (fm)

図 4.6: 中間子の古典場の Green 関数を利用した解2

h=0.01
h=0.5
h=0.1

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0

2

4

6

8

10

r (fm)

図 4.7: 中間子の古典場の反復解法

17

12

14

図 4.6 は中間子の古典場の Green 関数を利用した解から求めた中間子場であり、, 図 4.7 は中間子の古典場の反復解法から求めた中間子場であるが、グリッド の間隔を h=0.01,0.5,0.1 fm に変えたときに少しばらつきが見てとれるのは、中間子の古典場の Green を利用した解の方である。それは Green関数を用いて計算する場合、台形則による積分なので精度が低いことが短所である。0 < r 0 < r で rが小さい場合は低い精度になるのでこのような結果になったのだと推測される。

18

第 5 章 まとめ
自己無撞着に波動関数, 核子密度, 中間子場のそれぞれの解を求めることができることが知られているので、はじめにラグランジアン密度から出発し 、中間子を古典場に置き換え, さらに Dirac 方程式から導かれる負のエネルギーの核子の寄与を無視する, という二つの仮定して上で、Dirac 方程式の波動関数を動径波動関数と球面調和関数に分離して、動径関数を求めた。これにより波動関数から核子密度が得られる。そして、Green 関数を用いる計算法と、Gauss-Seidel 法を使う計算法で中間子場を求めた。この二つの計算法にはそれぞれ長所と短所がある。そこで、二つの方法を組み合わせることでより良い手法を考案した。すなわち、自己無撞着解を求めるための反復毎の変化の大きい初期の段階では、計算量の固定した Green 関数による方法を用い、反復毎の変化が小さい最終段階では収束が速くなっているはずの Gauss-Seidel 法を用いて計算するというものである。
今後の課題は、上記の中間子場を求めるための複合型手法を完成させ、さらに、自己無撞着に核子の波動関数と中間子場を求めるプログラムを完成させることである。

19

関連図書
[1] 有馬 朗人:「原子と原子核−量子力学の世界−」 朝倉書店, 1991

[2] 原 康夫 :「朝倉現代物理学講座 11 −素粒子−」 朝倉書店, 1981

[3] 中村 誠太郎:「大学院原子核物理」 講談社, 1996

[4] メシア著:小出昭一郎, 田村二朗 訳 :「メシア 量子力学 3」 東京図書株式会社, 1976

[5] 高橋 大輔:「数値計算」 岩波書店, 1996

[6] Brian D. Serot, John Dirk Walecka ”ADVANCES IN NUCLEAR PHYSICS”
PlenumPress, p38-84, 1986

[7] Peter Ring, Georgios A. Lalazissis, Dario Vretenar ”Relativistic Hartree-Bogoliubov
theory and application far from the valley of stability”
RIKEN Review, NO. 26, p10-16, January, 2000

[8] 奥村 晴彦: 「 C 言語による最新アルゴリズム事典」 技術評論社, 199120

謝辞
本論文を作成するにあたり，終始変わらぬ懇切丁寧なご指導と助言頂いた 鈴木敏男 先生, 林明久 先生, 田嶋直樹 先生には誠に感謝し ，お礼申し 上げます。田嶋直樹 先生には本研究に対し 、貴重な御時間を割いてまで根気強く丁寧に御指導頂いたこと、また本論文作成に関しても、多大な御協力を頂きました。また物理工学科の諸先生方に対しても、日頃の厚い御指導に深く感謝しています。また同じ研究室の院生の三和先輩、石田先輩やゼミのみんなにも、本当に感謝し 、謝辞の言葉とさせて頂きます。

2007 年 2 月
山田昌平

21

付録

Program List//mesonGF.c//Green’s function method to obtain the classical meson field//Poisson1: グリーン関数を使って核子密度から中間子場を求めるプログラム//rhocheck: 求めた中間子場を微分して核子密度を計算し , 入力された核子密度に//一致することを確かめる.
//2007 年 1 月 29 日

山田 昌平

作成

#include<stdio.h>#include<math.h>#defineN 3001//-------------------- 関数プロトタイプ宣言 ---------------------int Poisson1(int imax,double h,double m,double *rho,double *phi);double density(double r);
double rhocheck(double h,double m,double rho[],double phi[]);//--------------------------- main 関数-----------------------------

int main(){double m=0.7;
double h=0.01;//中間子の質量 (Mev/c^2)//動径グリッド 間隔 (fm)double r;
double rho[N];//動径座標//中間子の湧き出しdouble phi[N];int i,j;//中間子場for(i=0;i<N;i++){r=i*h;
rho[i]=density(r);}Poisson1(N,h,m,rho,phi);rhocheck(h,m,rho,phi);}//----------------------- Poisson1 関数 ------------------------------int Poisson1(int imax,double h,double m,double *rho,double *phi){double D;
double r0,r1;//Green 関数の値//Green 関数の変数 r,r’

double s;
int i,j;// r=0 での値
i=0;
r0=i*h;s=0;22//彼積分関数が原点について対称で遠方で速やかにゼロに収束するため、//台形則による積分で高い精度が得られる。
for(j=1;j<imax;j++){ // Point r1=0 does not contribute.
r1=j*h;s+=exp(-m*r1)*r1*rho[j];}phi[0]=s*h;// r>0 での値
for(i=1;i<imax;i++){r0=i*h;s=0;for(j=1;j<imax;j++){//彼積分関数が r=r’ で一階以上の微分が不連続となるため//台形則による積分結果は精度が低く、h の１乗に比例する。
r1=j*h;if(r1>=r0)
else

D=sinh(m*r0)*exp(-m*r1)/(m*r0*r1);
D=sinh(m*r1)*exp(-m*r0)/(m*r0*r1);

s+=r1*r1*rho[j]*D;}phi[i]=s*h;}
return 0;}//------------------------ 密度を求める式-----------------------------

double density(double r){
double R=7.0; //鉛の原子核 (A=208) の半径は 7.11fm
double a=0.7; //原子核の表面の領域の厚さのパラメータ (fm)return 1/((1+exp(r-R)/a)*(1+exp(-r-R)/a));}//--------------------

得られた密度の値をチェック --------------------

double rhocheck(double h,double m,double rho[],double phi[]){double ST;//中間子の souce term
double Fdphi;
double Sdphi;//中間子場 phi の一階微分//中間子場 phi の二階微分

double r;
int i;//【注】phi[0] が他の点と異なる (高精度の) 計算方法によるため。
//
phi[0] を含む有限差分は微分値とは大きく異なる値となる。
for(i=0;i<N-2;i++){r=i*h;
if(i==0){
Sdphi=(-phi[2]+16*phi[1]-30*phi[0]+16*phi[1]-phi[2])/(12*h*h);ST=-(3*Sdphi-m*m*phi[0]);}else if(i==1){
Fdphi=(phi[2]-phi[0])/(2*h);Sdphi=(-phi[2]+16*phi[1]-30*phi[0]+16*phi[1]-phi[2])/(12*h*h);ST=-(Sdphi+(2/r)*Fdphi-m*m*phi[i]);}else{Fdphi=(phi[i+1]-phi[i-1])/(2*h);Sdphi=(-phi[i+2]+16*phi[i+1]-30*phi[i]+16*phi[i-1]-phi[i-2])/(12*h*h);ST=-(Sdphi+(2/r)*Fdphi-m*m*phi[i]);}23printf("%f %e %e %e\n",r,rho[i],ST,phi[i]);}}//mesonGS.c//The Gauss-Seidel method the classical meson field//Poisson2: ガウス・ザイデル法により核子密度から中間子場を求めるプログラ厶
//2007 年 2 月 2-7 日 山田 昌平 作成#include<stdio.h>#include<math.h>#define N 3001//動径グリッド の点の個数。第 0 点が r=0 に対応する。

//-------------------- 関数プロトタイプ宣言 ---------------------int Poisson2(int imax, double h, double m, double *rho, double *phi );double density(double r);//--------------------------- main 関数
int main(){-----------------------------

double m=0.7;
double h=0.01;//π中間子の質量 (139.57MeV = 0.71 fm^(-1))//動径グリッド の間隔 (fm)double phi[N];double rho[N];//中間子場 φ (r)//中間子場の湧き出し (核子密度に比例する) ρ (r)double r;
int i;//動径座標 (r)for(i=0;i<N;i++){r=i*h;
rho[i]=density(r);phi[i]=rho[i];//Gauss-Seidel 法の反復のための初期値の設定}
Poisson2(N,h,m,rho,phi);for(i=1;i<N;i++){r=i*h;
printf("%f %e %e\n",r,rho[i],phi[i]);}}//--------------------- Poisson2 関数 -------------------------int Poisson2(int imax, double h, double m, double *rho, double *phi ){//Screened Poisson 方程式の解を Gauss-Seidel 法により求める。
//m
中間子の質量 (139.57MeV = 0.71 fm^(-1))//h=0.01;//phi[N];

動径グリッド の間隔 (fm)中間子場 φ (r)//rho[N];
中間子場の湧き出し (核子密度に比例する) ρ (r)double newphi; //更新後の中間子場 (作業用)double r;
double c=1.99;//動径座標 (fm)//Gauss-Seidel 法の減衰因子 (0<c<2,c>1 では加速)

double s,t,err;
double c0,c1,c2,z0,z1,z2;
int i,j;t=m*h; t=1/(t*t+2.5);c2=(-1.0/12.0)*t;24

c1=( 2.0/ 3.0)*t;
c0=

h*h*t;

t=m*h; t=1/(t*t+7.5);z2=-0.5*t;
z1= 8.0*t;
z0= h*h*t;// Gauss-Seidel 法の反復代入ループ

for(j=0;j<1000;j++){
s=0;

i=0;
newphi= z2*phi[2]+z1*phi[1]+z0*rho[0];t=newphi-phi[i];s+=t*t;phi[i]+=c*t;i=1;
t=2.0/i;newphi= c2*((1+t)*phi[i+2]+(1-t)*phi[1])+c1*((2+t)*phi[i+1]+(2-t)*phi[i-1])-c0*rho[i];t=newphi-phi[i]; s+=t*t; phi[i]+=c*t;for(i=2;i<imax-2;i++){t=2.0/i;newphi = c2*((1+t)*phi[i+2]+(1-t)*phi[i-2])+c1*((2+t)*phi[i+1]+(2-t)*phi[i-1])-c0*rho[i];t=newphi-phi[i]; s+=t*t; phi[i]+=c*t;}i=imax-2;
t=2.0/i;newphi= c2*((1+t)*0+(1-t)*phi[i+1])+c1*((2+t)*phi[i+1]+(2-t)*phi[i-1]) -rho[i]*c0;
t=newphi-phi[i];s+=t*t;phi[i]+=c*t;i=imax-1;
t=2.0/i;newphi= c2*((1+t)*0+(1-t)*phi[i+1])+c1*((2+t)*phi[i+1]+(2-t)*phi[i-1])-c0*rho[i];t=newphi-phi[i];s+=t*t;phi[i]+=c*t;

s=s*h;
if(s<1.0e-16) break;}if(s>=1.0e-16) fprintf(stderr,"warning: not converged: error=%e\n",s);err=sqrt(s*h);fprintf(stderr,"# #iteraction=%dreturn 0;

error=%e\n",j,err);}//----------------------double density(double r){
double R=7.0;
double a=0.7;

密度を求める式--------------------------//鉛の原子核 (A=208) の半径は 7.11 fm//原子核の表面の領域の厚さのパラメータ (fm)return 1/((1+exp(r-R)/a)*(1+exp(-r-R)/a));}25