核子系のアイソスカラー対相関

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

鈴木 紀史

目 次

1

序章

1 原子核の構造

31.1

固有値と固有関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31.2

殻模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.3 Nilsson 模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.3.1

4 重極変形 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.3.2

16 重極変形 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

1.3.3

スピン軌道結合 . . . . . . . . . . . . . . . . . . . . . . . . . . .
~l2 ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

1.3.4

2 アイソスピン

12
13

2.1

アイソスピン変数の導入 . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2

アイソスピン変数を導入した固有状態. . . . . . . . . . . . . . . . . .

15

2.3

核力の性質

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

16

2.3.1

π 中間子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.3.2

核力ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.3.3

荷電対称性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.3.4

荷電不変性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3 BCS 理論21

3.1

BCS 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.2

純粋な対相関力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.3

BCS 解を求めるプログラムの説明 . . . . . . . . . . . . . . . . . . . . .

26

4 異種核子間にも働くアイソスカラー対相関
4.1

アイソスピンの取り扱い方 . . . . . . . . . . . . . . . . . . . . . . . . .

1

34
34

4.2

4.1.1

アイソスピンという量子数 . . . . . . . . . . . . . . . . . . . . .

34

4.1.2

相互作用のアイソスピン依存性 . . . . . . . . . . . . . . . . . .

35

4.1.3

スピン, アイソスピン射影演算子の行列表現 . . . . . . . . . . .

36

4.1.4

相互作用の軌道運動部分の行列要素 . . . . . . . . . . . . . . . .

39

Hartree-Fock-Bogoliubov 法 . . . . . . . . . . . . . . . . . . . . . . . .

40

4.2.1

相互作用のアイソスピン依存性 . . . . . . . . . . . . . . . . . .

40

4.2.2

Hartree-Fock-Bogoliubov 方程式 . . . . . . . . . . . . . . . . . .

41

5 アイソスピン形式の HFB 法による対相関の計算

43

5.1

エネルギースペクトルの作製. . . . . . . . . . . . . . . . . . . . . . .

43

5.2

HFB 方程式の各行列要素の説明 . . . . . . . . . . . . . . . . . . . . . .

44

5.3

Jacobi 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

5.4

HFB 方程式の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

6 結果と考察

47

謝辞

56

参考文献

57

付録

Program List

58

i

序章
原子核における対相関は、質量公式における対項の存在をはじめ、偶々核のスピンが零であることや、回転準位から導かれた慣性能率が剛体値の数分の一しかないことなどの原因としてきわめて重要である。
原子核の対相関は、固体の超伝導現象の微視的理論である BCS 理論 [1] と同じ理論形式で扱うことができる。超伝導体中の電子対と同様に、陽子同士、中性子同士の対相関のみを考慮すれば十分であると考えられることが多いが、陽子と中性子が対を組むという原子核に特有の組み方もありうる。このような可能性を考える場合、アイソ
スピンの概念が重要である。
アイソスピン1 とは陽子と中性子の性質が電荷を除いてきわめてよく似ていることから、陽子と中性子を核子という粒子のアイソスピンという量子数の異なる状態とみるために導入されたものであるが、核力はアイソスピンを保存することがわかっているので、単なる便宜上の概念でなく本質的なものである。
陽子と陽子、中性子と中性子の対はアイソスピンが１であり、アイソベクトル対と呼ぶことができる。一方、陽子と中性子の対はアイソスピンが 0 の場合と 1 の場合があり、前者をアイソスカラー対と呼ぶ。一般に核力はアイソスカラー対間のほうが強く働く。

1970 年代までの陽子・中性子間の対相関に関する研究は、文献 [2] に数ページの簡潔で明解な review がある。それらの多くは軽い核の殻模型の相互作用にもとづいた研究であり、学ぶべき発見は多いが、その後 20 数年の計算機の進歩により、現在では研究で用いる計算の規模がまったく違っており、現在の理論的道具だてのなかで研究しなおす価値が年々たかまっていると思われる。
1

Isospin, isotopic spin, 荷電スピン。素粒子物理学では、アイソスピンを記号 I で表す習慣である

が、原子核物理学においてはアイソスピンは 記号 T で表す習慣である。これは記号 I が原子核の全角運動量を表す記号として使われるからである。角運動量が I である理由は、原子物理学で全角運動量 Jの一部として原子核の角運動量に I が使われているからである。また、原子核物理学ではアイソスピンの第 3 成分 が 12 のものを中性子に、− 12 のものを陽子に対応させるが、この対応は素粒子物理学の場
合と正反対である。

1

陽子・中性子間の対相関の研究は、その後しばらく研究が低調な期間をはさんで、1990 年代前後からの不安定核物理の実験的研究の進展にともない、質量公式におけるWigner 項 [3] をアイソスカラー対相関に帰すことができることが一部の人の注目を集めるようになった [4]。2000 年以降 Satula と Wyss らにより [5, 6]、原子核の回転状態の解析のために高度に発達してきたクランキング模型2 の理論を、荷電空間でのクランキング (強制回転) に翻訳して N = Z 核近傍の不安定核の性質の理解に役立てる論文が目を引いた。N = Z の核から出発して、中性子のフェルミ準位を高くし陽子のフェルミ準位を低していくことが、アイソ空間での回転の角速度を増大させていくことだと解釈できるのである。
最も一般的にアイソスピンの自由度を取り入れて対相関を扱う方法は、Hartree-Fock-

Bogoliubov(HFB) 法をアイソスピン形式で記述することである。そこで、本研究では、
この形式の HFB 法の計算プログラムを作成することを目的とする。そして完成したプ
ログラムによる計算例として、N = Z = 64 核をアイソ空間で回転させた場合に T = 1
対相関状態と T = 0 対相関状態がどのように移りかわるかを示すことにする。

2

Cranking model 「強制回転模型」または rotating shell model 「回転系での殻模型」と呼ばれる。

2

第1章
原子核の構造原子核の特徴として、密度の飽和性と結合エネルギーの飽和性がある。密度の飽和性とは質量数によらず密度はどんな原子核でもほぼ等しいということを意味し、結合エネルギーの飽和性とは質量数によらず 1 核子が持つ結合エネルギーは質量数にはほとんどよらないことを意味する。
原子核の構造に関する理論的模型として、液滴模型と殻模型があげられる。
構成粒子の個性を塗りつぶす考え方の液滴模型は、上記の 2 つの飽和性を示す最も簡単な系である。これは、平均的な量、例えば質量公式などを考えるときには理解しやすい模型である。しかし、原子核の魔法数を説明できないという欠点がある。
これに対し、原子核の性質を 1 粒子運動から出発して記述し説明しようというのが殻模型である。この模型は、魔法数の存在、閉核付近の原子核の諸説明の説明に驚くばかりの成功をおさめた。この模型は、核構造解析の有力かつ標準的な武器として広く使われている。

1.1 固有値と固有関数
まずこの節では、核子の固有値と固有状態がどのように表すことができるのかについて記していく。
ある１次元空間に粒子が存在しているとする。原点 O からの距離に比例する引力

Fx = −kx

(1.1)

k 2
x
2

(1.2)

はポテンシャル

V (x) =から導くことができる。力 −kx による単振動を古典的に扱うと、x = Acos(ωt + δ) を3

得るがこのとき角振動数は

s

ω=

k
m

(1.3)

で与えられる。ふつう −k の代わりに得られる k = mω 2 を用いることが多い。したがって、このような力を受けている粒子に対するシュレーディンガー方程式は!

h̄2 ∂ 2mω 2 2−+x ϕ(x) = εϕ(x)2m ∂x22(1.4)と表すことができる。このときのハミルトニアンは、

H=−h̄2 ∂ 2mω 2 2+x
2m ∂x22(1.5)

である。こうして、x 軸上で原点からの距離に比例する引力 − mω 2 x を受けて単振動している質点に対するシュレーディンガー方程式を得ることができる。これは、1 次元調和振動子に体するシュレーディンガー方程式とも言える。
また、物理的に意味のある解を与えるような ε の値をε1 , ε2 , ε3 , · · ·εn , · · ·(1.6)とし、これらそれぞれに対応する解 ϕ(x) を

ϕ1 (x), ϕ2 (x), ϕ3 (x), · · ·, ϕn (x), · · ·(1.7)と記すことにする。これらはHϕn (x) = εn ϕn (x)(1.8)

を満たす。一般に、演算子を関数に作用させると別の関数が得られるわけであるが、特にそれが元の関数の定数倍になる場合に、その関数をその演算子の固有関数、その定数を固有値という。シュレーディンガー方程式を解くことは、ハミルトニアン H の固有値（エネルギー固有値）ε1 ,ε2 ,··· と固有関数 ϕ1 (x),ϕ2 (x),··· を求める固有値問題である。また、ｎ は正整数であるから、(1.8) 式の ε はあきらかにとびとびなものに限られる。これはちょうどスペクトルのように見えることから、固有値スペクトルとも呼ばれる。このようにとびとびの固有状態に番号づけする数は量子数と呼ばれる。
これまで、1 次元空間についての固有値と固有関数を記述してきた。それでは、取り扱う空間を 1 次元から 3 次元に拡張してみることにする。この場合、シュレーディンガー方程式は
(

h̄2
−
2m

∂2
∂2
∂2
+
+
∂x2 ∂y 2 ∂z 2!

)mω 2 2+(x + y 2 + z 2 ) ϕ(r) = εϕ(r)
2
4

(1.9)

である。

r 2 = x2 + y 2 + z 2(1.10)と書くと、この方程式はそれぞれ、x, y, z 方向の運動に関する独立な 3 つの部分の和として書ける。そのため変数分離ができて、φ(xyz) = φnx (x)φny (y)φnz (z)(1.11)

3h̄ω
2n0 = nx + ny + nz
ε = n0 +(1.12)(1.13)であることがわかる。ただし、ε=E+V0 である。V0 は原子核の中心でのポテンシャルの深さを表す定数を示す。最もエネルギーの低い状態は、量子数が (nx , ny , nz ) = (0, 0, 0)の状態であり、3
ε = h̄ω(1.14)
2を持つ [7]。この場合でも nx ,ny ,nz は正整数であるから、(1.12) 式の ε はあきらかにとびとびなものに限られる。
ところで、中心力ポテンシャルをもつシュレーディンガー方程式の固有状態は球座
標を用いて表すこともできる。この場合量子数は、動径方向の波動関数の節点の数 n、軌道角運動量 l、その z 成分 m である。このとき l は、

l = 0, 1, 2, 3, · · ·(1.15)と整数のみを取る。また m は z 成分であるから、絶対値が l が超えることはないので、

m = −l, −l + 1, −l + 2, · · ·, l − 1, l

(1.16)というとびとびの値をとる。こうすることにより粒子の状態は n, l, m の 3 つの量子数でも記述することが可能である。なお、

nx + ny + nz = 2n + l(1.17)という関係がある。

1.2 殻模型
多数の実験事実により、原子核の陽子数 Z 、中性子数 N のいずれかが、


Z 

N 

= 2, 8, 20, 28, 50, 82, 126

5(1.18)

をもつ原子核は、結合エネルギーが大きく、安定していることが結論づけられている。
これらの数を、原子核の魔法数 (magic number) という [7]。
ところで、原子にも魔法数が存在し、原子番号 2,10,18,36,54,86, の希ガス元素は電子軌道の閉殻にあたり安定になる。このことから魔法数の存在は、原子核にも閉殻が存
在し、核を構成している核子が 1 つの共通のポテンシャル場の中でそれぞれの軌道を描き、独立の運動をしている殻模型が成立することを示唆していると考えられる。

1.1 節で述べた調和振動子ポテンシャルにより導出した式 (1.14) の場合を考えてみる。この状態にスピン上向き、下向きで、2 個の陽子または中性子が入り、この軌道を満員にする。これが、Z ＝ 2、または、N ＝ 2 が、魔法数である理由である。このように、一つの軌道が満員になることを閉核になると言う。次に、N ＝ 1 の核では、(1, 0, 0), (0, 1, 0), (0, 0, 1)、と言う 3 つの状態が縮退している。スピンの上, 下を考慮にいれると、この軌道には陽子または中性子が 6 個入れる。そこで、2+6 ＝ 8 だから、Z ＝ 8
または N ＝ 8 が魔法数になる。全く同様に、Z(N ) ＝ 8+12 = 20、Z(N ) ＝ 20+20 ＝ 40、さらに、70、112 が魔法数になることが予想される [7]。ところが、こうして導かれた
魔法数は 40 以上については実際の数とは異なっている。
では次に、平均ポテンシャルとして調和振動子型でなく、原子核の密度分布に良く似た関数形をもつ密度に比例したポテンシャルである Woods-Saxon ポテンシャルを導入してみる。このポテンシャルは下記の関数形をもつ。

V (r) = −V0
1 + {exp(r − R)/b}(1.19)

ここで b は、核表面のぼやけを表すパラメータである。このポテンシャルでは、調和振動子ポテンシャルでおきたような縮退は解け、エネルギー準位は分裂する (図 1.1 参照)。しかし、このポテンシャルを用いても 28 以降の魔法数を説明することはできない。
なお、軌道角運動量の大きな軌道は、遠心力で外側に密度が多くなるため、ポテン
シャルが調和振動子から Woods-Saxon 型に変わることで大きくエネルギーが下がる (図1.1 参照)。このことから、Woods-Saxon ポテンシャルは、調和振動子ポテンシャルに
軌道角運動量の 2 乗に比例するポテンシャルを加えたもので近似できることがわかる。

結局、もっともらしいどのような形の平均ポテンシャル V (r) を用いても原子核の魔法数を再現することはできない。 この矛盾は、1949 年、Mayer と Jensen によって独立に解決された。彼らは、核子の内部角運動量であるスピンと軌道角運動量を結合させるスピン軌道力が、原子核内の平均ポテンシャルでは非常に強いと考えると、魔法

6

数が説明できることを発見したのである。スピン軌道力は、

ξ(~l · ~s)(1.20)と書ける。~l は軌道角運動量を、~s はスピン角運動量を表し、ξ は定数であるとする。式(1.9) に、Woods-Saxon ポテンシャルを近似する ~l2 項とこの項を付け加えると、)(h̄ ~ 2 mω 2 2∇ +r + a~l2 + ξ(~l · ~s) ϕ = εϕ−2m2(1.21)となる。この式の固有状態は、n, l, j = l ± 12 , m、という 4 つの固有値で指定される。
図 1.1 は、一粒子エネルギー準位が、それぞれのポテンシャルによって分裂していく様子を表している。記号 s,p,d,f ,g,h,i は、l=0,1,2,3,4,5,6 を表し、例えば 1h11/2 は、n =1,l = 5,j = 11を表す。
2
このとき、ξ < 0 とすると、~l · ~s > 0 の方がエネルギーが低くなる。~l · ~s > 0 のときは、j = l + 1 、~l · ~s < 0 のときは、j = l − 1 になるので図のように分裂する。ここで、
2

2

N ＝ 3 に属する準位は 1f7/2 ,2p3/2 ,2p1/2 ,1f5/2 、の 4 本の準位に分裂する。このうち一番低いのは 1f7/2 である。j ＝ 27 の状態は、m ＝ ± 72 ,± 25 ,± 32 ,± 12 をもつから、8 個の陽子または中性子で閉殻となり、魔法数 20+8=28 が得られるのである [7]。同様に 50、82、126 の魔法数も説明できる。

7

図 1.1: 1 粒子エネルギー準位の分裂 (文献 [7] より転載)

8

図 1.2 のように、エネルギーの低い準位から順番に P auli 原理にしたがって核子を詰めていけば、原子核全体としての最低エネルギー配位ができる。これを基底状態配位という。これを基準にして、非占有状態からなる 1 粒子状態を粒子空間、占有準位からなる空間を空孔空間と定義する。基底状態配位は、場の理論での真空に対し粒子も空孔も存在しない状態とみなせる。他方励起状態は、図 1.2 のように粒子空孔ペア − が何個か生成された多粒子多空孔状態として記述できる [9]。

図 1.2: 殻模型における粒子-空孔励起 (概念図)(文献 [7] より転載)

9

1.3 Nilsson 模型
原子の場合には、中心にある重い原子核が強い力で電子の群の集団運動を統率しているが、原子核にはこのような中心体がない。そのことからも予想できることである
が、原子核は容易に変形することができる [9]。原子核の変形を考えれば、当然一粒子運動として球対称でない力の場の中の運動を考えることになる。つまり、他の核子からの相互作用を球対称でない平均ポテンシャルで近似し、個々の核子がこのポテンシャ
ル内で独立に運動しているという描像を仮定することになる。このように。球対称でないポテンシャル中での運動へ殻模型を拡張したのが、Nilsson 模型である [8][10]。
この模型では、4 重極変形、16 重極変形、スピン軌道結合、~l2 ポテンシャルで構成される Nilsson ポテンシャルが用いられる。このポテンシャルは、

1
2 2
2 2
U (r) = m(ω⊥x + ω⊥y + ωk2 z 2 )2 s4π+ 2h̄ω0 rt2
ε4 Y40 (r̂)
9+ 2κh̄ω0~lt · ~s22
− κµh̄ω0 (~lt − < ~lt >N )(1.22)と、表される。1.3.1∼1.3.4 節で、上式右辺の各項について説明する。
なお、Nilsson 模型の与えるのは、個々の核子のエネルギーである。個々の核子のエネルギーをすべての核子について足し合わせると、相互作用のエネルギーが二重に取り込まれるため、原子核全体のエネルギーにはならない。しかし、原子核に核子一個を加えたときの原子核全体のエネルギーの変動は、個々の核子のエネルギーの和の変動によって良く近似できる。

1.3.1 4 重極変形
Nilsson はまず、4 重極変形を持つ核の１粒子状態をよく表しているものとして、原子核の変形が軸対称であるとし、軸対称の調和振動子型ポテンシャルで粒子軌道を考えた。このポテンシャルは通常の殻模型のときの振動子型ポテンシャルに対応している。

1
U (r) = m(ω12 x2 + ω22 y 2 + ω32 z 2 )2(1.23)

対称軸が第 3 軸であるとすれば、ω1 = ω2 = ω⊥ である。また、ω3 = ωk とする。
変形のパラメータ ε2 を導入し、

1
ω⊥ = ω0 (1 + ε2 )3
10

(1.24)

となり、ε2 が小さいときは2ωk = ω0 (1 − ε2 )3のように ω⊥ ,ωk を変化させることにする。

(1.25)

式 (1.24), 式 (1.25) から、変形パラメータ ε2 > 0 のときは偏長変形、ε2 < 0 のときは偏平変形、ε2 = 0 のときは球形を表すことが分かる。ここで ω0 は、h̄(ω1 ω2 ω3 )1/3 =

N −Z
41
1
∓
A−1/3
3A



[MeV]

になるように、各 ε2 の値に対して決める。ただし、陽子の運動を考える場合は + であり、中性子の運動を考える場合には − である。このように ω0 を決めると原子核の密度の飽和性を満たすように、変形によらず原子核の体積を一定にすることができる。
図 1.3 は、変形殻モデルの 1 粒子エネルギー準位で、50 < Z < 82 領域の陽子に対する計算例である。スペクトルの疎らな場所が安定と考えられるので、ε2 = 0 の場所では、50,82 の陽子数で安定になっている。しかし、ε2 = 0.2 あたりでは、66 の陽子数で安定になる。また、ポテンシャルの変形によっても 1 粒子エネルギーが変わる事も分かる。

図 1.3: 変形核モデルの 1 粒子エネルギー (概念図)(文献 [7] より転載)

11

1.3.2 16 重極変形
このポテンシャルは、4 重極変形の補正として導入されているものであり、変形度 ε4は 4 重極変形に比べて 1 桁程度小さくなる。ここで Y40 (r̂) は、球面調和関数であり、ε4は、各 ε2 の値に対しエネルギー最小値を与えるように決めるのが一般的であるが、本研究では簡単のため ε4 = 0 とする。

rt は、stretched coordinate での動径長を表し、小さな変形では、rt ' r と考えてよい。

1.3.3 スピン軌道結合
球形殻模型の場合と同様に、スピン軌道力は必要である。κ は、ポテンシャルの強度を表すパラメータである。本研究では κ = 0 とする。

1.3.4 ~l2 ポテンシャル
1 粒子準位の順序や準位間隔を、Woods-Saxon ポテンシャルのものに近づけるために、~l2 に比例する項を更に付け加えている。µ はポテンシャル強度パラメータである。

12

第2章
アイソスピン
原子核は主に、陽子と中性子から構成されている。陽子と中性子はひとまとめにして核子と総称する。この 2 種類の粒子は、非常に似た性質を持っていることが知られている。すなわち質量についていえば、陽子および中性子の質量 Mp , Mn はそれぞれ、

Mp = 1.67239 × 10−27 (kg)(2.1)

Mn = 1.67470 × 10−27 (kg)(2.2)

であり、質量はほぼ等しい。またこの 2 種類の核子は、電荷が異なるという特徴を持っている。つまり、同一粒子の異なる電荷状態であると言える。このため、電子のスピ
ン z 成分を区別するスピン演算子にならって、陽子と中性子を z 成分の異なる固有状
態として区別しようと考えられた。こうして、スピンと区別してこの自由度を特にア
イソスピン（荷電スピン）と呼んでいる。またこの理論は核子以外の素粒子にも拡張され、素粒子を特徴づける基本量のひとつにもなっている。

2.1 アイソスピン変数の導入
先ほども述べたように、陽子と中性子は質量がほとんど等しく電荷が異なる粒子である。この両者を荷電状態が異なる１種類の粒子として取り扱うためには、荷電状態を示す新しい変数を導入する必要がある。陽子の電荷 e を単位にとれば、中性子と陽子の電荷は 0 および 1 である。よってこれを固有値とする変数を考えてもよいが、両者を対称的に扱うため 1, −1 を固有値とする変数τz を導入し便宜上 − 1 を陽子に + 1 を中性子にあてはめることにする。このような変数の導入により核子は 2 成分の波動関数で表される。したがって、この 2 つの成分に関係する演算子は 2 行 2 列の行列で表され13

る [9]。

τx = 


0 1 





1 0



 0 −i τy = 

i 0


(2.3)



1

0 

0 −1

τz = 


式 (2.3) を τ スピンまたはアイソスピン行列と呼び、この演算子で表される物理量をア
イソスピンと呼ぶ。
これらの行列は、3 次元のアイソスピン空間内でのベクトルの x,y,z で示される軸方向の成分とみなされる。対角である z 成分は、上に導入したアイソスピン変数に対応する。
行列 (2.3) からアイソスピン演算子

1
t= τ
2

(2.4)

が得られる。この演算子の成分は角運動量ベクトルの交換関係を満たす。(t)2 = 34 であるから、核子は全アイソスピン t = 12 、z 成分 mt = tz = + 21 (中性子に対し) 、および

− 21 (陽子に対し) を持つ。
またスピンにおいても、変数σ を定義しその x, y, z 成分はτ と同じ行列で表される。


 0 1 

σx = 
σy = 

1 0

0 −i 

i 0

(2.5)



 1

σz = 

0 

0 −1

スピン演算子は σ と次式で関連付けられる。

1
s= σ
2
が得られる。
式 (2.3), 式 (2.5) は Pauli 行列とも呼ばれている。

14

(2.6)

2.2 アイソスピン変数を導入した固有状態
2 つの核子 1,2 があり、それぞれの軌道角運動量を l1 , l2 、z 成分を m1 , m2 とする。2 核子系の軌道角運動量の和を L と書くと、

L = l1 + l2

(2.7)という関係が成り立つ。そのとき z 成分は

M = m1 + m2

(2.8)

である。古典力学であれば合成した角運動量 L = |L| は、l1 + l2 から |l1 − l2 | まで自由な値を取れる。しかし、量子力学では様子が違い L は

L = l1 + l2 , l1 + l2 − 1, l1 + l2 − 2, · · ·, |l1 − l2 |(2.9)

である。スピン角運動量についても同じであり、2 電子のスピンの合成の結果は

S=

1 1
1 1
1 1
+ , S = + − 1, · · ·, S = −
2 2
2 2
2 2

(2.10)

S = 1, 0(2.11)

すなわちの 2 つの値を取る。軌道角運動量 L とスピン角運動量 S の和 J を全角運動量と呼ぶ。J は

J =L+S(2.12)

で与えられ、ふたたび式 (2.9) 式の場合のように

J = L + S, L + S − 1, · · ·, |L − S|(2.13)というとびとびの値を取る [7]。
このようにして 2 核子系の状態は、軌道角運動量 L、全スピン S 、全角運動量 J を使って分類できる。それではさらに、アイソスピンをも変数として考えてみる。このとき 2 核子系の全アイソスピンは、

T = t1 + t 2 =τ1 + τ22(2.14)となる。t1 , t2 はそれぞれ核子 1, 2 のアイソスピンである。またT は代数的には全スピ
ンS と同じ構造を持っている。例えばT の各成分は互いに可換ではないがT 2 とは可換である。したがってT 2 とT の成分の 1 つ、例えば Tz との同時固有状態を作ることがで15

きる。このときの Tz の固有値 Mt の最大値を T とすると、T 2 の固有値は T (T + 1) である。この T をT の大きさというのが普通である。そこで、状態を全アイソスピン T の値で分類すると、2 核子系では T = 1 または 0 ということになる。陽子 − 陽子系、中性
子 − 中性子系は T = 1 の状態であり、陽子 − 中性子系は T = 0 と T ＝ 1 の両方の状態がありうる。

2.3 核力の性質
原子核のような多体系を扱う場合に、多体力が存在するか、存在するとすればどの程度の寄与を現象に与えるかは検討しなければならない問題である。しかしながら、これに関する知識は現在のところ実験的にも理論的にも不確かな点が多い。したがって、核子間の相互作用の問題を考える場合、まず 2 体の相互作用だけで考えてみることにする。

2.3.1 π 中間子
原子核内での 2 核子間の距離は、多少の重なりを考慮に入れて 2fm 程度である。陽子や中性子をこのくらいの距離に保ちつつ束縛させる核力の起源は、いったい何であろうか。この問いかけに対し 1934 年、湯川秀樹は 2 粒子間である未知の粒子が交換され核力が生じると考え、この粒子を中間子と名付けた。現在では湯川中間子は π 中間子と呼ばれ、他にも多数の中間子が発見されている。

π 中間子には電荷が±e のもの (π ± ) と中性なもの (π 0 ) があり、それぞれの質量は、mπ± c2 = 139.57

mπ0 c2 = 134.96[MeV],[MeV](2.15)である。

2.3.2 核力ポテンシャル
2 つの核子 1, 2 があり、核子間距離を r とする。このとき核子の間に働くポテンシャルの概形を図 2.1 に示す。

16

図 2.1: 核子と核子の相互作用のポテンシャルの依存性
この図におけるいくつかの基本的な特徴の現れを次に示す。

• 相互作用は約 1 fm 程度の短距離力である• この距離においては約 40 MeV の深さで引き付けられる• 約 0.5 fm より短距離になると反発力が強くなる
湯川理論によると r≥2fm の遠距離では 1 個の π 中間子の交換によるポテンシャル

OPEP(One-Pion Exchange Potential) により、
(

!

)

1 f2
3
3
e−µr
VOP EP =mπ c2 (τ1 · τ2 ) (σ1 · σ2 ) + 1 +
+
S
12
3 h̄c
µr (µr2 )µr(2.16)となる。添字 1, 2 は核子の番号を示す。S12 は、

S12 = 3(σ1 · r)(σ2 · r)
− σ1 · σ2r2

17(2.17)

で定義される。mπ はπ 中間子の質量 (mπ c2 '140MeV)、µはπ 中間子のCompton 波長 mh̄π c '1.4fmの逆数である。また f はπ 中間子と核子の結合定数で、π 中間子の核子による散乱実験か2ら fh̄c '0.08 である。VOP EP の形から分かるように、核力は相互作用している 2 核子の波動関数がどのようなスピン, アイソスピン状態にあるかに依存する。また、相対運動の
軌道角運動量にも強く依存する。このような「強い状態依存性」は核力の重要な特徴である。
原子核という多体系の中での核力は、さまざまな多体効果のために 2 核子系での核力と異なる。これを核内での有効相互作用という。核力の強い状態依存性を反映して、
有効相互作用の性質は複雑になる。1960 年代以降、核力と核構造の関連についてさま
ざまな研究が積み重ねられ、かなりの性質についてかなりの知見が得られた。しかし、今なお多くの未解決の問題を残している。

核子は、その電気的な性質を別にすれば、陽子と中性子を同等と考えてよい。また、核子間の力も、電磁相互作用によるものを別にすれば、核子の荷電状態に無関係であるらしい。そこで、原子核系で電磁相互作用を取り去ると、系はその荷電状態に無関係である、すなわち、荷電空間での回転に対して不変であると仮定することができる。
核力は電磁気的な力よりずっと強いから、原子核系は主としてこの荷電不変な相互作用に支配され、この不変性が原子核の性質に現れることになる。2.3.3 節、2.3.4 節では、
荷電対称性と荷電不変性という核力の性質について述べる。

2.3.3 荷電対称性
π 中間子にはπ ± , π 0 の 3 種類があり、核子には陽子と中性子がある。それぞれ電荷が違うにもかかわらず質量の差はわずかである。とすると、核力が荷電状態によらないことが予想される。
まず 2 陽子間, 2 中性子間の核力はπ 0 の交換で生じるが、核子とπ 中間子の結合の強さが電荷によらなければまったく同じ力になる。このように、中性子と陽子をまったく入れかえたとき核力が同じであることを荷電対称性という。
荷電対称性がどの程度成り立っているかの例として、15 N と 15 Oの準位構造を図 2.2 に示す。

18

図 2.2: 15 Nと15 O の準位構造図 (文献 [9] より転載)
これによれば、両方の核の準位は互いに対応するものが大体同じエネルギーのところに現れていることがわかる。エネルギー値が少しずつ食い違っているのは、状態による核の Coulomb エネルギーの差などの高次の補正項によるものと考えられてる [9]。

2.3.4 荷電不変性
中性子と陽子間ではπ ± が交換される。中性子と陽子からなる系には、

n(1)p(2) + p(1)n(2)√2(2.18)

n(1)p(2) − p(1)n(2)√(2.19)
2の 2 状態がある。ここで n(i), p(i) は、i 番目の核子がそれぞれ中性子, 陽子の状態にあることを意味している。π ± を交換して中性子は陽子にまたはその逆になるから、1 つの核子を時間を追ってみていると陽子になったり中性子になったりする。そこで、このように n(1)p(2) と n(2)p(1) の 2 状態を混ぜておくのである。2 陽子系 p(1)p(2), 中性子系 n(1)n(2) が番号 1 と 2 の入れかえに対して不変なことから、式 (2.18) の状態
は nn, pp の状態とよく似ていて、その状態のとき np 間の核力が nn, pp の場合と同じになることが予想され、実験でも確かめられた。これを核力の荷電不変性と呼んでいる。
例として,14 C(MT = 1),14 N(MT = 0),14 O(MT = −1) の 3 つの組を考えてみる。図 2.3 にこの 3 つの組の核のスペクトルを示す。

19

図 2.3: 14 C(MT = 1),14 N(MT = 0),14 O(MT = −1) の準位構造図(文献 [9] より転載)
14

Nの基底状態は T = 0 であり、2.31Mev の O+ 準位は、14 C,14 Oの基底状態と共に T =1 の準位であると考えられる。14 Nの基底状態をエネルギーの原点に取ると,14 C,14 Oの基底状態はそれぞれ0.156,5.15MeV のところにある。これに、中性子-陽子質量差 1.29MeVを考慮に入れると, これらと 14 N の 2.31MeV 準位との差は、それぞれ 3.44MeV および

4.13MeV である。一方、Coulomb エネルギーの差は, 上の式では, それぞれ、3.8MeV お
よび 4.3MeV となる。定性的な傾向としては大体一致している。定量的には,Coulomb
エネルギーの式に核の大きさの実験値を用いるとか、核子間相関を考慮に入れるなどにより, もっとよい一致を得ることができる [9]。

なお、この節の執筆にあたっては主として参考文献 [9] の第 2 章を参考にした。

20

第3章

BCS 理論ある種の金属または合金の電気抵抗が、その物質に固有な温度以下で 0 になる現象を、
超伝導、その状態を超伝導状態という。これは、1911 年、H.Kamerlingh-Onnes によって 4.2[K] 以下の水銀について発見され以来、近年の高温超伝導ブームに続く長い実験的探求により様々な物質がなりうる状態であることがわかってきた。
ところで、軽い核と閉殻近傍の核を除いて、ほとんどの原子核の基底状態は超伝導状態になっていると思われる。すなわち、フェルミ面近傍の中性子と陽子がそれぞれ、
クーパー対と呼ばれる対を形成している。原子核のクーパー対は、運動量が逆向きではなく、角運動量の方向が逆向きの核子どうしが強く相関することにより形成される。
結合エネルギーの偶奇効果、偶々核と偶奇核の間の励起スぺクトルの相違、変形核の慣性モーメントが剛体直の 1/2〜1/3 に減少しているなどの実験事実は、核が超伝導状態になっていることの反映である。

BCS 理論とは、1957 年、J.Bardeen、L.N.Cooper、および J.R.Schrieffer によって提唱された超伝導の微視的理論である。ここでは、彼らの推測による超伝導体での基底状態の決定にならい、以下の式で偶々核の場合の波動関数がよく表されるとする。

|BCS >=

Y+(uk + vk a+k a−k )|− >

(3.1)

k>0

式 (3.1) の中の uk ,vk が変分パラメータである。k は一粒子状態のラベルであり、ゼロでない整数値をとるものとし、k でラベルされる状態があれば、−k でラベルされる状態もあるとする。時間反転不変なハミルトニアンに対しては、状態 k と 状態 −k は時
間反転操作によってお互いに変換される状態であるようにラベルを割り振る。例えば、個体中の電子のような第ゼロ近似として空間的に一様で並進対称なポテンシャルのもとで運動する粒子を扱う場合は、k を波数ベクトル、sz をスピンの z 成分として、状態 |k, sz i と状態 | − k, −sz i とが時間反転状態対となる。一方、原子核のような第ゼロ近似として空間的に局在した球対称なポテンシャルをもつ状態の場合は、全角運動量21

を j 、その z 成分を m として状態 |j, mi と |j, −mi とが時間反転状態対となる。ただし、対を組む状態間の相対位相を首尾一貫してつけることに注意を払わなければならない。式 (3.1) の乗積記号

Qの添字の k は、正の値を網羅する。u2k ,vk2 は、状態の対での占有確率を表す。これらに状態のノルムが依存しないように、|u2k | + |vk2 | = 1(3.2)という条件を課す。この理論は原子核の全ての核子を扱うことを許し、そして異なったタイプの相互作用の場合に、容易に一般化することができる。

3.1 BCS 方程式
多体系がハミルトニアン

H=X

tk1 k2 a+
k 1 ak 2 +k1 k2 6=0

X1+
v k 1 k 2 k 3 k 4 a+
k 1 ak 2 ak 4 ak 3
4 k1 k2 k3 k4 6=0(3.3)

よって記述されると仮定する。試行波動関数 (3.1) のパラメータ、u と v は、エネルギーの変分によって決定される。しかしながらこの変分は、粒子数の期待値が目標値 N に一致すつべしという、拘束条件によって制限される。

< BCS|N̂ |BCS >= 2X

vk2 = N

(3.4)

k>0

この条件は、変分ハミルトニアンに項 −λN を加えることにより満たすことができる。

H 0 = H − λN̂(3.5)

ここで、ラグランジュ乗数 λ は、条件 (3.4) によって固定される。E =< BCS|H|BCS >
とすると、

dE
(3.6)
dN
が成り立つことが示せるので [11]、λ は粒子数の変化に伴うエネルギーの増加を表す
λ=

ことが分かる。したがって λ は化学ポテンシャル、あるいはフェルミ準位と呼ばれる。

(3.1) 式、そして (3.3) 式から、H 0 の BCS 期待値として、



X


X
1 X
2 2
2
0
0
v
+
v̄
v
v̄kk̄k0 k¯0 uk vk uk0 vk0< BCS|H 0 |BCS >=
+
(t
−
λ)v
0
kk kk k kkkk
2 k0 6=0
k6=0
kk0 >0(3.7)

式を得る。BCS 波動関数はパラメータ vk そして (3.2) 式によって完全に決定されるので、変分

δ < BCS|H 0 |BCS >= 0
22

(3.8)

は、

(

∂
duk ∂
+
) < BCS|H 0 |BCS >= 0∂vk
dvk ∂uk(3.9)

をもたらす。これを微分することにより、BCS 方程式は、

2˜k uk vk + ∆k (vk2 − u2k ) = 0 ,k>0(3.10)となる。ここで、˜k 、そしてギャップパラメータ ∆k は、





X
1
˜k =
tkk + tk̄k̄ +(v̄kk0 kk0 + v̄k̄k0 k̄k0 )vk20 − λ

2
k0 >0
X

∆k = −v̄kk̄k0 k¯0 uk0 vk0(3.11)(3.12)k0 >0である。˜k と ∆k の値を固定すれば、(3.2) 式と (3.10) 式により、u2k と vk2 がそれぞれ以下の解を持つ。



vk2 =



1
˜k

1± q
2
2
˜ + ∆2k(3.13)

k





1
˜k

u2k = 1 ± q
222
˜ + ∆
k(3.14)k

相互作用がない場合は、∆ = 0 だから占められた軌道 (˜ < 0) は、vk2 = 1,u2k = 0 を持つ。それゆえ式 (3.13) と式 (3.14) において唯一の可能な解は、




1
˜k

vk2 = 1 − q
2˜2 + ∆2k

(3.15)

k





1
˜k

u2k = 1 + q
2
22
˜ + ∆
k(3.16)k式である。従って、変分原理 (3.8) 式は、方程式 (3.11)、(3.12) そして、(3.15)、(3.16)をもたらす。粒子数条件、

2X

vk2 = N

(3.17)k>0

式と共に、それらは BCS パラメータ uk 、vk の計算を可能にする。一般に、これらの方程式は、非線形であり、そして反復的解法によって解かなければならない。これらの方程式での性質の議論の場合、式 (3.15), 式 (3.16) を式 (3.12) に挿入して得られるいわ
ゆるギャップ方程式∆k = −∆k0
1 Xv̄kk̄k0 k¯0 q
2 k0 >0˜2 0 + ∆2 0kがしばしば有用である。

23

k(3.18)

3.2 純粋な対相関力
純粋な対相関力は、原子核の対相関の説明において、非常に単純だが強力なモデルを提供する。以下では、この特別な相互作用の場合に BCS 理論の式がどのような形をとるかを述べる。まずハミルトニアンは、

H=

X
k>0

X

+
a+
k ak̄ ak̄ ak0(3.19)

1∆2(˜k vk2 + Gvk4 ) −2
G
k>0(3.20)+
k (a+k ak + ak̄ ak̄ ) − G

kk0 >0となる。H 0 の期待値式 (3.7) は、
X< BCS|H 0 |BCS >= 2である。この場合、ギャップパラメータ ∆ は、k によらない。

∆=GX

uk v k 0(3.21)

˜k = k − λ − Gvk2

(3.22)k>0

また、˜k が、によって与えられる。(3.20) 式での項 Gvk4 と (3.22) 式での項 Gvk2 は、しばしば無視される。なぜなら、その効果は一粒子エネルギーをシフトさせることであるが、このような図式的な相互作用を使用する場合、既に対相関を導入する前の段階で一粒子エネル
ギーは現実の原子核を記述するのに適切と思われる値に設定されており、その値をさらに変更することは原子核の記述精度をむしろ低下させることになるからである。したがって、より正確に言えば、これらの項は無視されるのでなく、一粒子エネルギーの決定時にすでに考慮されているので二重には入れないことにするということである。
この場合、




u2k 





1
k − λ= 1 ± q
2 
2
2
2

vk
(k − λ) + ∆(3.23)

式を得ることができる。また、この式は、極限 G → 0、すなわち ∆ → 0 において、占拠されるべきレベルは vk2 → 1 となり、空席であるべきレベルは vk2 → 0 となっているべしという物理的な要請を満足していることが容易に見て取れる。この場合、vk2 はステップ関数であるのに対して、相互作用しているケース (∆ 6= 0) はステップ関数が、幾分ぼやけたものになる。相互作用のために、粒子はフェルミ準位の下から上へ散乱していく。これはフェルミレベルの下に状態の部分的な減損を、上に状態の部分的な満足をもたらす。最期に、ギャップ方程式 (3.18) は、単純な形

∆=

GX
∆
q
2 k>0 (k − λ)2 + ∆2
24

(3.24)

をとることになる。次節では、BCS のギャップ方程式を反復法を用いて解き、ギャップ
パラメータやフェルミ準位の様子を観察してみることにする。

25

3.3 BCS 解を求めるプログラムの説明

図 3.1: 一粒子スペクトルの設定

26

図 3.1 のような原子核の殻模型の特徴に似せた 1 粒子エネルギースペクトルを考え、BCS のギャップ方程式を解くことにより、ギャップパラメータやフェルミ準位の様子を観測する。ギャップ方程式は以下の式で与えられる。各準位は 2 重に縮退していて 1 粒子状態ラベル i は i = ±1, ±2, · · ·, ±Ni をとる。




Ni
G X
1

q
∆=−
2 i=1 (i − λ)2 + ∆2

(3.25)

下記は計算のために用いた C 言語プログラミングのうち、エネルギースペクトル iを作製するための部分と、(3.25) 式を反復代入法を用いて計算するための部分である。

++++++++++++++++++++++++++++++++++++++++++++++++++++
エネルギースペクトルを作製++++++++++++++++++++++++++++++++++++++++++++++++++++
void eps_calculation(double eps[] , int nl , int ns , double e_bace ,double d1 , double d2){
int i;
eps[0]=e_bace;for(i=0;i<=nl-1;i++){if(((i+1)%ns)==0){
eps[i+1]=eps[i]+d2;
printf("%d %lf\n",i+1,eps[i]);}else{
eps[i+1]=eps[i]+d1;
printf("%d %lf\n",i+1,eps[i]);}}}++++++++++++++++++++++++++++++++++++++++++++++++++++
式 (3.25) を反復代入法を用いて計算++++++++++++++++++++++++++++++++++++++++++++++++++++
double delta_calculation(double eps[] , int nl , double* delta ,double g_int , double* lamda){
int i,k=0;
double wa,sol=0,sa_2=1.0,sol_2,delta_2;
while(1.0e-10<=fabs(sa_2) && k<1000){
k++;
wa=0;
for (i=0;i<=nl-1;i++){wa+=1.0/(sqrt(pow((eps[i]-lamda),2)+pow(delta,2)));}sol=(-g_int/2.0)*wa*delta ;sol_2=fabs(sol);delta_2=fabs(delta);
sa_2=delta_2-sol_2;
delta=sol;27

}if(k>=1000){delta=0.0;
puts("delta couldn’t be calculatied!!\n");}return delta;}+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Program 1
上記のプログラムの各変数は、次のような意味を持つ。

eps[i]= εi(エネルギースペクトル)(3.26)

ebace =

(エネルギー基準値)(3.27)

d1=(主殻内の準位間隔)(3.28)

d2=(主殻間エネルギーギャップ)(3.29)

delta=∆(対相関ギャップ)(3.30)

gint=G(対相関相互作用強度)(3.31)(フェルミ準位)(3.32)

lamda = λi の値を 0 から nl まで求め、N で番号づけする。i = 0 から N まで、ループで回し、
wa =

N
X

q

i=1

sol = −

1

(i − λ)2 + ∆2G
∗ wa ∗ ∆2.0(3.33)(3.34)を求める。sol と ∆ の絶対値の差が、1.0 × 10−10 より小さくなるまで繰り返す。このようにして ∆ を算出する。
ギャップ方程式 (3.25) が右辺を計算して得られた ∆ の値を再び右辺に代入することを繰り返すという「反復代入法」で解ける理由は、右辺を ∆ の関数として見たときの傾き d(右辺) が正で 1 より小さいからである。これが反復代入法の収束条件を満たし
d∆ていることについては例えば文献 [12] を参照されたい。

28

図 3.2: 一粒子スペクトル
図 3.2 は、 NL = 12、NS = 4、D1 = 0.5、D2 = 2 と与え、一粒子スペクトルを設定したものである。この一粒子スペクトルをもとに G = −0.3 とし、ギャップ方程式 (3.25)を計算した。

29

図 3.3: フェルミ準位と対相関ギャップエネルギーの関係

図 3.4: フェルミ準位と粒子数の関係
図 3.3 は、横軸にフェルミ準位 λ をとり縦軸に対相関ギャップエネルギー ∆ を、図

3.4 は、横軸にフェルミ準位 λ をとり縦軸に粒子数の期待値 N をとったものである。図3.3 から図 3.2 で設定したスペクトルの主殻内にフェルミ準位のあるとき、対相関ギャッ
プエネルギーはゼロでない値をとり、主殻間のスペクトルの空隙にフェルミ準位のあるとき対相関ギャップエネルギーはゼロになることが分かる。一つのスペクトル準位には 2 個の粒子が入ることが可能である。今の計算では、1 つの殻内に 4 本の準位を設定したため、図 3.4 では 8 個づつ粒子が詰まっては、殻間のギャップのためしばらく休むといういく様子が見てとれる。

30

次に算出された ∆ を (3.23) 式に代入し vk2 を求め、式 (3.17) を用いることにより粒子数 N を求める。下記は、粒子数 N を求めるプログラムの部分と λ を 2 分法を用いて求めるプログラムの部分である。

++++++++++++++++++++++++++++++++++++++++++++++++++++
粒子数 N を求めるプログラム++++++++++++++++++++++++++++++++++++++++++++++++++++
double n_calculation(double eps[] ,int nl , double delta , double g_int, double lamda){
int i;
double v_1,wa,n_1,du;
wa=0;
for (i=0;i<=nl-1;i++){du=1.0-(((eps[i])-lamda)/(sqrt(pow((eps[i]-lamda),2)+pow(delta,2))));
v_1=(1.0/2.0)*du;
wa+=v_1;}
n_1=2*wa;return n_1;}++++++++++++++++++++++++++++++++++++++++++++++++++++
2 分法でフェルミ準位を決定するプログラム++++++++++++++++++++++++++++++++++++++++++++++++++++
double lamda_calculation(double eps[] , int nl , double n_0 , double* delta ,int k=0;
double m_n=100.0,l_n,h_n;double l_delta,m_delta,h_delta;
double l_lamda=-100.0,m_lamda,h_lamda=100.0;
while(fabs(m_n-n_0)>1.0e-8 && k<1000){k++;m_lamda=(l_lamda+h_lamda)/2.0;
m_delta=delta_calculation(eps , nl , &delta , g_int , &m_lamda);
m_n=n_calculation(eps , nl , m_delta , g_int , m_lamda);
if(m_n<n_0){
h_lamda=m_lamda;}
else{
l_lamda=m_lamda;}}
return m_lamda;}++++++++++++++++++++++++++++++++++++++++++++++++++++
Program 2

ギャップ方程式プログラム

i = 0 から N までループで回し、
i − λdu = 1 − q(i − λ)2 + ∆2
31

(3.35)

v1 =

1.0
∗ du
2.0

wa =

nl−1
X

(3.36)

v1(3.37)

n1 = 2 ∗ wa

(3.38)

i=0を求める。
先ほど算出された ∆、そして N により、2 分法によって λ が算出される。図 3.2 の一粒子スペクトルをもとに、粒子数 N フェルミ準位 λ を計算した。

32

図 3.5: 粒子数と対相関ギャップエネルギーの関係

図 3.6: 粒子数とフェルミ準位の関係
図 3.5 は、横軸に粒子数 N そして縦軸にギャップエネルギー ∆ を、図 3.6 は、横軸に粒子数 N そして縦軸にフェルミ準位 λ をとったものである。先ほども述べたが、今計算では一つのバンド内に 4 本のスペクトルを設定したため、１つのバンド内には 8 個の粒子が入ることが可能である。また図 3.6 からは、8 個の粒子が増えるごとにフェルミ準位が不連続に増加する様子が見てとれる。

33

第4章
異種核子間にも働くアイソスカラー対相関系を構成している粒子の運動状態を考えるとき、それが他の粒子の状態と無関係に記述できるならば、その粒子は他粒子と相関がないといえる。粒子間に相互作用が働いている系では粒子間相関が存在するのが普通である。直接に相互作用がなくても、Pauli
原理のように互いに他の運動を制約する法則があるときには、やはり相関が現われる。
相互作用のある系ではごく一般的にいえば、すべての粒子は互いに相関を持っているから、系が N 個の粒子から成り立っていれば N 粒子相関まで考えなければならない。
しかし、作用している力の性質、系全体のエネルギー値などによって相関の程度は変わってくる [9]。
核子間の相互作用は短距離力でしかも強いから、2 体相互作用が重要であることは容易に推測される。原子核内の核子間の平均距離は 2[fm] 程度で、核力のかなり弱い領域に相当しているから、3 粒子以上の相関は 2 体相互作用に比べてずっと弱いと考えても不当ではない。以下では 2 体相互作用を考えていくことにする。

4.1 アイソスピンの取り扱い方
4.1.1 アイソスピンという量子数
原子核内における核子、陽子と中性子の区別は電荷によってすることができた。つまり、陽子は電荷 + e を持ち、中性子は電荷をもたない。この両者を荷電状態が異なる 1 種類の粒子として取り扱うために、荷電状態を示す新しい変数τ を導入することにした。便宜上、核子が全アイソスピン t = 21 のとき、その z 成分が陽子の場合 tz = + 21
、中性子の場合 tz = − 21 とおく。こうすることにより、原子核内における 2 種類の核子をはっきりと区別して示すことが可能になった。

1 粒子の固有状態を記述するために、電子の場合には軌道運動 n, l, m、そしてスピ34

ン s の 4 つの量子数を与える必要があった。これと同様に核子の場合には、陽子と中性子を区別する必要性から、n, l, m, s, t という 5 つの量子数を与えることにより、1 粒子の固有状態を指定できる。
このように考え、ある核子 i の波動関数は第 2 量子化を使って、|i >= |n, l, m, s, t >

(4.1)

と書く。

4.1.2 相互作用のアイソスピン依存性
2 核子系について核子間に相互作用がある場合を考えてみる。2 核子系の状態 1, 2 が相互作用により 3, 4 になったとする。このとき核力ポテンシャル V1234 は、

V1234 = < n1 l1 m1 s1 t1 , n2 l2 m2 s2 t2 |V̂ |n3 m3 l3 s3 t3 , n4 l4 m4 s4 t4 >(4.2)

と書くことができる。上式は、2 核子の始状態が相互作用により終状態に変化することを表す式に相違ない。
では次に、相互作用を 2 核子の全スピン S 、全アイソスピン T 、全軌道角運動量 L の値で分類してみる。テンソル力を無視すれば、S, T, L は保存されるからである。このときの V を VT SL とすると、

V̂ =とおける。ここで
Q

Q

Y

Y
T

Y
S

L VT SL(4.3)

T は 2 核子のアイソスピンが T である部分空間への射影演算子でQある。また、 S はスピン S の空間への、 L は相対角運動量 L の空間への射影演算子である。それでは全角運動量 L を、

L = l1 + l2

(4.4)

とおき、さしあたり強く働く L = 0 のみ (s 波) を考える。全スピン S は L によらず常に、S = |s1 + s2 | = 1, 0(4.5)

の値しかとりえないので、全角運動量 J は、J = |L + S| であり空間の等方性により当然保存される。よって、L = 0 のときは J = S であり、S は保存される。即ちテンソル
力を無視しなくとも、L = 0 のみを考えるときは S は保存される。また全アイソスピ
ン T も、

T = t1 + t2 = 1, 0
35

(4.6)

の値をとる。よって、式 (4.3) は、

V̂

=

Y

+

Y
T =0

S=0 V̂T =0,S=0,L=0 +

Y

Y
T =0

Y

S=1 V̂T =0,S=1,L=0 +Y
T =1

S=0 V̂T =1,S=0,L=0

Y

Y
T =1

S=1 V̂T =1,S=1,L=0

(4.7)

となる。
核子はフェルミ粒子である。このことから、2 核子の状態は 2 核子の入れ換え操作に対して反対称でなければならない。

表 4.1: 2 核子の入れ換え操作に対する対称性
対称

T =1 S=1

L = 0, 2, 4, ...

反対称

T = 0 S = 0 L = 1, 3, 5, ....

このように考えると全角運動量 L = 0 のとき可能な組合せは、T = 1, S = 0 また
は T = 0, S = 1 の 2 チャンネルしか取りえない。結局、L = 0 の場合相互作用を書き下
すと、

V̂ =

Y

Y
T =0

S=1 V̂T =0,S=1,L=0 +YY

S=0 V̂T =1,S=0,L=0

T =1

(4.8)

となる。式 (4.2) と式 (4.8) より 2 核子系の状態を表す相互作用行列要素 V1234 は、

V1234 = < t1 t2 |Y

T =0 |t3 t4 >< s1 s2 |Y

S=1 |s3 s4 >< n1 l1 m1 , n2 l2 m2 |V̂010 |n3 l3 m3 , n4 l4 m4 >+ < t1 t2 |Y

T =1 |t3 t4 >< s1 s2 |Y

S=0 |s3 s4 >< n1 l1 m1 , n2 l2 m2 |V̂100 |n3 l3 m3 , n4 l4 m4 >(4.9)

と記述することができるのである。4.3,4.4 節では、上式をさらに詳しく考察していくことにする。

4.1.3 スピン, アイソスピン射影演算子の行列表現
前節では、波動関数のアイソスピン部分を |t > と記述することにより、原子核における 2 種類の核子、陽子そして中性子の区別をすることができた。本節では、陽子に対する波動関数のアイソスピン部分を |p > 、中性子に対する部分を |n > とし、2 種類の核子を区別してみる。このように陽子と中性子を区別することにより、核子が対を組むケースを以下の 4 通りに場合分けすることが可能になる。

|p1 > |p2 >,|p1 > |n2 >,|n1 > |p2 >,36|n1 > |n2 >(4.10)

前節で述べたように、2 核子系の状態は全アイソスピン T の値で分類すると、T =0 と T = 1 の両方の状態がありうる。また、T の角運動量を表す成分 Tz は、T = 0 の場
合 Tz = 0、T = 1 の場合 Tz = 0, −1, 1 である。このことと、式 (4.10) より、以下の式を導くことができる。

1|t = 0, tz = 0 > = √ (|p1 > |n2 > −|n1 > |p2 >)21|t = 1, tz = 0 > = √ (|p1 > |n2 > +|n1 > |p2 >)2|t = 1, tz = 1 > = |p1 > |p2 >(4.11)(4.12)(4.13)|t = 1, tz = −1 > = |n1 > |n2 >(4.14)

それでは、上式を用いることにより式 (4.9) の内部の < t1 t2 |Q

T =0 |t3 t4 > と t1 t2 |Q

T =1 |t3 t4 >部分を計算してみる。まず、アイソスピン交換演算子 Pτ を定義する。Pτ は、Y

1(1 − Pτ )2
Y1
T =1 = (1 + Pτ )2
YY
T =0 +
T =1 = 1(4.15)

T =0 =(4.16)(4.17)という条件を持つ。また、Pτ と、式 (4.10) の 2 核子の状態との乗算はそれぞれ、

Pτ |p1 > |n2 >= |n1 > |p2 >(4.18)

Pτ |n1 > |p2 >= |p1 > |n2 >(4.19)

Pτ |p1 > |p2 >= |p1 > |p2 >(4.20)

Pτ |n1 > |n2 >= |n1 > |n2 >(4.21)となる。式 (4.14) と、アイソスピン交換演算子 Pτ から与えられた条件式、式 (4.17), 式 (4.21) よ
り、
Y
Y
Y
Y

T =0 |T = 1, Tz = 0, ±1 >= 0

(4.22)

T =1 |T = 1, Tz = 0, ±1 >= |T = 1, Tz = 0, ±1 >(4.23)

T =0 |T = 0, Tz = 0 >= |T = 0, Tz = 0 >(4.24)

T =1 |T = 0, Tz = 0 >= 0

(4.25)という式を導くことができる。式 (4.9) の内部の < t1 t2 |Q

T =0 |t3 t4 > と t1 t2 |Q

T =1 |t3 t4 >は、上式より求めることができる。表 4.1, 表 4.2 に 2 つの式の、おのおので求まる値を記す。

37

表 4.2: < t1 t2 |表 4.3: < t1 t2 |Q

Q

T =0 |t3 t4 > の値

T =1 |t3 t4 > の値

この表は、縦軸に t1 t2 を横軸に t3 t4 をとったものである。また、式 (4.9) の内部の

< s 1 s2 |Q

S=0 |s3 s4 >, < s1 s2 |Q

S=1 |s3 s4 > は 2 核子系の状態のスピンに対応する部分であり、これはアイソスピンと同様に振る舞うため、表 4.1, 表 4.2 の値と同様なものになると考えて良い。
ここまでで、式 (4.9) の 2 つの項の内部における、前半部分の係数の値を求めることができた。次節では、式 (4.9) の後半部分の係数を求めてみることにする。

38

4.1.4 相互作用の軌道運動部分の行列要素
あらためて、式 (4.9) の第 1 項の後半部分の係数を書き出してみる。

< n1 l1 m1 , n2 l2 m2 |V̂010 |n3 l3 m3 , n4 l4 m4 >(4.26)

上式は、Clebsch-Gordan 係数を用いることにより、

σl0 m0 |n1 n2 l0 m0 >< l1 m1 l2 m2 |l0 m0 >σlm |n3 n4 lm >< l3 m3 l4 m4 |lm >(4.27)と記述することができる。上式をまとめて書くと、X< l1 m1 l2 m2 |l0 m0 >< l3 m3 l4 m4 |lm >< n1 n2 l0 m0 |V̂T S0 |n3 n4 lm >(4.28)

lml0 m0となる。VT SL は、全スピン S 、全アイソスピン T 、全軌道角運動量 L の値で分類していることを考慮すると、V̂100 では L = 0 なので軌道角運動量 l = l0 = 0 のときのみ有効であり、その他では 0 になる。また軌道角運動量が 0 ということは、その z 成分も 0 になるので m = m0 = 0 である。よって、式 (4.9) の 2 つの項の後半部分 をまとめて書くと、< n1 l1 m1 , n2 l2 m2 |V̂T S0 |n3 l3 m3 , n4 l4 m4 >= < l1 m1 l2 m2 |00 >< l3 m3 l4 m4 |00 >< n1 n2 00|V̂T S0 |n3 n4 00 >(4.29)となる。
また、< l1 m1 l2 − m2 |00 > という Clebsch − Gordan 係数は、(−1)l1 −m1< l1 m1 l2 − m2 |00 >= δl1 l2 δm1 −m2 √2l1 + 1(4.30)という簡単な形をとることができる。ただし、この式は l1 = l2 かつ m1 = −m2 が必要条件であり、これ以外の値をとると 0 になってしまう。よって、クロネッカーデルタを用いることにより、δl1 l2 δm1 −m2 という因子がかけられている。< l3 m3 l4 − m4 |00 > も同様の形をとる。以上のことをまとめることにより、式 (4.29) は、(−1)l1 −m1 +l3 −m3
V̂T S0 (l1 n1 n2 l3 n3 n4 )
δl1 l2 δm1 −m2 δl3 l4 δm4 −m4 q(2l1 + 1)(2l3 + 1)

(4.31)となる。またV̂T S0 (l1 n1 n2 l3 n3 n4 ) は厳密に言えばすべての l1 , n1 , n2 , l3 , n3 , n4 の組合せについて計算しなければならないが、本研究では、研究の第一歩として、近似的対相関相互作用としてよく利用される seniority force と同様の考え方で、同じ主殻内の２39

粒子が対を組み、その対のまま別の主殻内に散乱されていくという形の行列要素のみを零でないと仮定することにする。すなわち、

V̂T S0 δn1 n2 δn3 n4(4.32)に近似するのである。
後述するように、主殻の異なる軌道にある２核子が対を組めればアイソスカラー対相関が N 6 Z の核でも強く残存する可能性が高まるので、研究の次のステップでは、厳密に行列要素を計算することが大切である。したがって、本研究で開発する HFB 法のプログラムは、行列要素に関する seniority force 型の仮定にほとんど依存しないように作られている。
以上、3, 4 節で式 (4.9) の内部の説明をしてきた。これらより、相互作用行列要素 V1234は
、

V1234

(−1)l1 −m1 +l3 −m3= δl1 l2 δm1 −m2 δl3 l4 δm4 −m4 δn1 n2 δn3 n4 q
(V 010 < t1 t2 |Y

+V 100 < t1 t2 |

(2l1 + 1)(2l3 + 1)

T =0 |t3 t4 >< s1 s2 |YY

T =1 |t3 t4 >< s1 s2 |

S=1 |s3 s4 >Y

S=0 |s3 s4 >)(4.33)と書くことができる。

4.2 Hartree-Fock-Bogoliubov 法
原子核にはすべての核子によって作られている 1 体ポテンシャルの中を、核子が独立に運動しているという描像が近似的に成立している。A 個からなる原子核の中の核子は他の A − 1 個の作る平均 1 体場内を運動しているが、この核子は同時に他の核子の感じる平均 1 体場を作り出している。有限系に特有なこの自己無撞着性を満たすな 1 体場の微視的基礎づけを与えるのが、Hartree-Fock(以後、HF) 理論である。しかし HF 理論
では、原子核の比較的近い励起状態で重要な対相関を取り扱えない [13]。したがって、

HF 理論を拡張した Hartree-Fock-Bogoliubov(以後、HFB) 理論について述べる必要がある。

4.2.1 相互作用のアイソスピン依存性
前節まで議論していた Vijkl を 2 体行列要素として< ij|V̂ |kl >40(4.34)

とする。またこのとき、反対称化 2 体行列要素を、< ij|V̂ |kl >a =< ij|V̂ |kl > − < ij|V̂ |lk >(4.35)と定義する。HFB 理論では準粒子 a+
i , ai がBogoliubov 変換

a+
i =

X

(Uij c+j + Vij cj )(4.36)(Uij∗ cj + Vij∗ c+
j )(4.37)

j

ai =

X
j

のように定義される。この a+i , ai が反交換関係を満たすためには、式 (4.36), 式 (4.37) のユニタリ性が必要であり、これを行列 U, V の満たす関係として与えれば

UU+ + V V + = 1(4.38)

UV + + V U+ = 0(4.39)を得る。ここで (U + )ij = Uij+ とする。本研究では実数しか現れないので転置と同じことである。粒子 − 空孔密度行列ρij と対テンソル tij とを、

ρij = (V + V )ij(4.40)

tij = (V + U )ij(4.41)

で導入する。ρij は粒子 − 空孔相関を、tij は粒子 − 粒子相関を平均 1 体場に取り込むのに基本的な役割を果たす。このとき 2 体相互作用 < ij|V̂ |kl > に含まれる粒子 − 空孔
相関と粒子 − 粒子相関とを各々 1 体ポテンシャルとして取り込んだHartree-Fock ポテ
ンシャル, 対ポテンシャルは

Γij =

X

< ik|V̂ |jl >a ρlk(4.42)

1X< ij|V̂ |kl >a tkl
2 kl(4.43)

kl∆ij =で与えられる。

4.2.2 Hartree-Fock-Bogoliubov 方程式
この多体系がハミルトニアン

Ĥ =

X
ij6=0

Tij c+
i cj +

1 X+
Vijkl c+i c j c k al
4 ijkl6=041(4.44)

によって記述されたとする。Tij は、1 粒子の運動エネルギーである。Ĥ は粒子数を保存する。すなわちĤ は粒子数演算子
X

c+
i ci

(4.45)< N̂ >= N̄

(4.46)N̂ ≡iと可換であるといえる。拘束条件を満たすために式 (5.10) にフェルミエネルギーλを導入することでハミルトニアンはĤ 0 = Ĥ − λN̂(4.47)となる。ところで、陽子数をN̂p 、中性子数をN̂n とすると

N̂p =

X

c+
i ci

(4.48)

c+i ci(4.49)

i

N̂n =

X
i

とおける。N̂p は 1 粒子状態のアイソスピン成分 t ＝ −1 のときだけ値が入り、N̂n は t ＝ 1 のときだけ値が入る行列である。よって、陽子のフェルミエネルギーをλp 、中性子のフェル
ミエネルギーをλn として、

Ĉ = λp N̂p + λn N̂n(4.50)

Ĥ 0 = Ĥ − Ĉ

(4.51)

とおく。これより式 (4.47) はとおくことができる。また、Ĥ は 1 粒子エネルギーをi とおくと

Hij = (Tij + Γij )ij = δij i(4.52)と近似する。i には Nilsson ポテンシャルのエネルギー準位を用いる。

HFB 方程式は U, V と準粒子エネルギー Ei を定める固有方程式であり Ĥ − Ĉ


ˆ∗
−∆





ˆ
∆
−(Ĥ − Ĉ)∗

 Ui 
  Ui  = Ei 


Vj

Vj

(4.53)

で与えられる。上式から明らかなように、結果として得られる固有ベクトルを HF ポ
テンシャル, 対ポテンシャルが含んでいるという非線形な固有値問題となっており、実際には反復計算で数値的に求めることができる。

42

第5章
アイソスピン形式の HFB 法による対相関の計算
HFB 法を行列の個有値問題として解くにあたり、まず 1 粒子状態のエネルギースペク
トルの作製方法について説明する。

5.1 エネルギースペクトルの作製
核子の固有状態は、主量子数 n、軌道角運動量 l、その z 成分 m、スピンの z 成分 s
、アイソスピンの z 成分 t で表すことができることはこれまでに述べてきた。では、実際の 1 粒子固有状態がどのように表されるか解説していく。
例えば 1 粒子軌道量子数の最大値を Nmax = 1 とする。これは n = 0 と n = 1 の場合があることを意味する。n = 0 の場合は、l = 0、m = 0、s = −1, 1、t = −1, 1 であり全部で 4 通りの固有状態がありえる。次に n = 1 の場合は、l = 1、m = −1, 0, 1、s = −1, 1、t = −1, 1 であり全部で 12 通りの固有状態がありえる。よって、1 粒子軌道量子数の最大値が Nmax = 1 の場合は、16 通りの固有状態がありえることになる。ここで、s, t は、

s = 2szt = 2tz(5.1)という条件を持っている。このときの n, l, m, s, t で表される核子は 1 粒子エネルギーei (1 <= i <= 16) を持っている。このエネルギー ei を!

n(n + 3)−ei ＝h̄ω (n + 1.5) − κµ l(l + 1) −

2.0

s



5
l(l + 1) − 3m2 β(n + 1.5)4π(2l − 1)(2l + 3) 
(5.2)という式で与える。上式の各変数は、次のような値を入れる。

h̄ω = 7.0 [MeV] (調和振動子主殻間のエネルギー間隔)
43

(5.3)

κµ = 0.025(ポテンシャル強度パラメータ)

β = 0.22

(変形度)

(5.4)
(5.5)

計算のために用いた C 言語プログラミングのうち、エネルギースペクトル  を作製するための部分は、付録に載せたプログラムに含まれている。

5.2 HFB 方程式の各行列要素の説明
HFB 方程式は、
 Ĥ − Ĉ


ˆ∗
−∆



ˆ
∆
−(Ĥ − Ĉ)∗







  Ui  Ui 
 = Ei 


Vj

Vj

(5.6)

で与えられる。(本研究で用いた値からは虚数は現れないため、∗ は考えないものとする
ˆ は 1 粒子軌道量子数の最大値を Nmax とすれば、Nmax ×Nmax の行列になって)Ĥ, Ĉ, ∆,ˆ の作製方法について説明していくことにする。
いる。以下ではĈ, Ĥ, ∆•行列Ĉ の説明
4.2 節で述べたようにĈ = λp N̂p + λn N̂n(5.7)

である。N̂p は、ti = −1 のとき 1 を, ti = 1 のとき 0 を、N̂n は ti = −1 のとき 0 を, ti = 1 のとき 1 を対角要素に入れていく対角行列である。よって、Ĉ は ti = −1 のときλp を、ti =1 のときλn を対角要素に入れていく対角行列であると言える。
•行列Ĥ の説明
4.2 節で述べたようにĤ = δij i(5.8)

である。i は 1 粒子エネルギースペクトル作製時にできるのi 値を使用する。δij よりĤ は、対角要素にi の値が入る対角行列であると言える。

ˆ の説明
•行列∆
4.2 節で述べたように
∆ij =

1X
< ij|V̂ |kl >a tkl
2 kl

44

(5.9)

である。< ij|V̂ |kl > は 4 章で求めた

V1234(−1)l1 −m1 +l3 −m3= δl1 l2 δm1 −m2 δl3 l4 δm4 −m4 δn1 n2 δn3 n4 q
(V010 < t1 t2 |Y

+V100 < t1 t2 |

(2l1 + 1)(2l3 + 1)

T =0 |t3 t4 >< s1 s2 |YY

T =1 |t3 t4 >< s1 s2 |

S=1 |s3 s4 >Y

S=0 |s3 s4 >)(5.10)

である。本研究では V010 = −3.0, V100 = −2.0 として計算した。また、U, V は




ˆ
∆

 Ĥ − Ĉ


∗
∗
ˆ
−∆
−(Ĥ − Ĉ)(5.11)の 固有ベクトルより求まる値である。U + , V +は U, V の転置である。計算初期では U, V は値がないため、これらより求まる行列要素 tkl は乱数を使い 0.5 ∼ 0.6 として作製した。
計算のために用いた C 言語プログラミングは、付録に転載した。

5.3 Jacobi 法
行列要素の固有値と固有ベクトルを求める方法として、本研究では Jacobi 法を採用した。Jacobi 法とは、実数の対称行列の固有値、固有ベクトルを同時に全て求める簡単な反復法である。
実数の対称行列 A を適当な直行行列 X1 により A − X1T AX1 と変換し、A のどれか 1 つの非対角成分を 0 にする。次にまた適当な直行行列 X2 により別の非対角成分を 0 にする。その際、以前に 0 にした成分は 0 でなくなるが、最初よりは 0 に近いことが多い。このようなもぐら叩きを、A のすべての非対角成分がほぼ 0 になるまで繰り返すと
、対角化
T
Xm
...X2T X1T AX1 X2 ...Xm = A(5.12)が完成する。A の非対角成分 apq を 0 にする直行行列 X は pq 平面内の回転

xpp = cos θ,
xpq = − sin θ,

xpq = sin θ
xqq = cos θ(5.13)とし、変換後の apq が 0 になるという条件から、θを定める。A の非対角成分を消す順序
については、(1) 絶対値の大きい非対角成分から順に消していく。
(2)0 でない非対角成分を片っ端から消していく。
(3) 絶対値がある値 K 以上の非対角成分を片っ端から消し、K を次第に 0 に近づける。などの方法が考えられる。本研究では文献 [14] のものを利用した。

45

5.4 HFB 方程式の計算
それでは、HFB 方程式を Jacobi 法と反復代入法を用いて解いてみることにする。5.2
ˆ を用いて節で説明した Ĥ, Ĉ, ∆


f =


ˆ
∆



Ĥ − Ĉ


∗
∗
ˆ
−∆
−(Ĥ − Ĉ)(5.14)

を作製する。この行列から Jacobi 法を使用することにより、固有値と固有ベクトルを求める。
このとき求まった固有値を仮に Ei1 とする。固有ベクトルから U, V を求めることが出来る。しかし、1 回目の計算ではの初期値は任意に決めていたため求まった値は真値とˆ を新たに作製することによって f をは言いがたい。それゆえ求まった V + , U を使って∆作り、再びJacobi法により固有値と固有ベクトルを求める。Ĥ とĈ の行列は、一度作製されればその後、値の変更はない。

2 回目のjacobi法より求まった固有値を、仮に Ei2 とする。こうして固有値 Ei1 と Ei2の絶対値の差を求める。このときの差が 1.0 × 10−4 MeVより小さくなるまで計算を繰り返す。このようにして固有値 Ei と固有ベクトル、そして U, V を算出する。

46

第6章
結果と考察
この章では、付録に載録した自作のプログラム isopair.c を用いて行なった数値計算の結果の一例を解説する。
計算では、一粒子状態として、３次元調和振動子の主量子数 N が 0 から 5 までの 56状態に、スピン・アイソスピン自由度 4 を乗じた 224 個の状態を取り入れた。主殻間のエネルギー間隔 を h̄ωosc = 7 MeV、Nilsson ポテンシャルの l2 項の係数を κµ = 0.025、
四重極変形度を β2 = 0.22 とした。取り入れた状態のエネルギーのスペクトルを図 6.1に示す。

図 6.1: 1 粒子エネルギースペクトル
中性子および陽子のフェルミ準位を

λ n = λs +ωτ,2

λ p = λs −

ωτ
,
2

(6.1)

と表し、λs = 40 MeV に固定し、アイソ空間のクランキング (角) 周波数 ωτ を 0 MeVから 10 MeV まで掃引して、各々の ωτ の値に対して HFB 解を求めた。ωτ = 0 MeV

47

の場合は λn = λp であり、N = Z 核に対応する。ωτ の増大に伴って、λn は上昇し λpは下降するので、N > Z の核に対応するようになる。
対相関相互作用の強度 VT SL は V010 = −1.0 MeV、V100 = −0.66 MeV とした。２つの強度の値の比は文献 [5, 6] を参考にして決めた。強度の大きさは対相関の働く空間との兼ね合いで適当な値に決めた。
このようにして計算した HFB 法の基底状態での対相関エネルギー

Epair =

1X
∆ij t∗ij
2 ij

を ωτ の関数として描いたのが図 6.2 である。

48

(6.2)

図 6.2: 基底状態での対相関エネルギーωτ < 2.4 MeV では Epair = −40.0 MeV という一定値をとっている。この値を T = 0の相互作用からの寄与と T = 1 の相互作用からの寄与に分けて計算したところ、100%が T = 0 の寄与であることがわかった。ωτ = 2.4 MeV で不連続に減少して、それよ
り右側では、100% が T = 1 の寄与からなることがわかった。
なお、転移点 ωτ = 2.4 MeV の左側の部分で一定値であるのは、アイソクランキン
グという外場の与え方の特殊性によるものである。また、この転移の位置には履歴性はなく、左側から計算しても、右側から計算しても同じ位置で転移が起きる。
図中の点線は T = 1 の相互作用のみを考慮して解を計算した場合の Epair である。
転移点の右側の領域では T = 1 と T = 0 の両方の相互作用を取り入れて求めた解と一致している。左側では、完全な相互作用を用いた場合と違い、T = 1 対相関は消えることなく存続している。
逆に、T = 0 の相互作用のみを用いて解を求めた場合は、転移点の左側では完全な相互作用を用いた場合と一致し、右側では、対相関は完全にゼロになってしまう。

49

図 6.3: 中性子および陽子の個数の期待値
図 6.3 は、中性子および陽子の期待値を ωτ の関数として描いたものである。T = 1の相互作用のみを用いて解を求めた場合はなめらかに変化するのに対し、完全な相互作用を用いた場合は、T = 0 対相関の区間では中性子と陽子とは同じ個数であり、転移点で不連続に変化して個数が違ってくることがわかる。

50

図 6.4: 全エネルギー
図 6.4 は対相関エネルギー Epair に平均場エネルギー

Emf =

X

ρii i(6.3)

iを加えた全エネルギー Etot = Emf + Epair である。やはり転移点において不連続な変化をしていることがみてとれる。
なお、全エネルギーの順位が逆転する点と転移点が一致しないのは、解を求める基
準が全エネルギーでなく、拘束条件のための付加項を含んだもの (Routhian と呼ばれる) だからである。
また、図 6.4 では、転移以外にも曲線のうねりが気にかかるかもしれないが、これは、アイソクランキングに伴うフェルミ準位と１粒子準位との交差に起因するもので、とりたてて興味深い情報は含んでいない。

51

図 6.5: 正常状態の全エネルギー
図 6.5 には対相関力を全く入れずに求めた解の全エネルギー (=平均場エネルギー)が描いてあるが、この変化が対相関力のある場合のうねりの元であるとみることができる。

52

図 6.6: Routhian の期待値
このようなフェルミ準位の変動によるエネルギーのうねりを除去して、対相関によるエネルギーの変化がどれだけか明確に見てとれるようにするためには、Routhian と呼ばれる下式で定義される演算子

R̂ = Ĥ − λp N̂p − λn N̂n

(6.4)

0
+ Epair ,R = hR̂i = Emf

(6.5)

0
= Emf − λp hN̂p i − λn hN̂n i
Emf

(6.6)

の期待値

をプロットするとよい。というのは、HFB 解を求めるという変分操作の対象が Etot =hĤi ではなく R = hR̂i であるため、図 6.5 の Etot (に含まれる Emf ) の示すような不連
0
続な変化は、R(に含まれる Emf) には存在しないからである。

図 6.6 は R を ωτ に対してプロットしたものである。実線は対相関力を考慮しない場合に得られる正常状態解の R である。T =1 対間にのみ相関力が働くようにすると、

53

破線が得られる。実線と破線との間隔が、アイソベクトル対相関によるエネルギーの利得に相当する。ダイヤ印 (◇) は T =0 対の間にのみ対相関力が働くようにして得られた解である。ダイヤ印をつなぐことで、T =0 対相関解は R = −1213 MeV の水平な直線で表されることがわかる1 。実線で表された正常状態とこの水平線とが交差するωτ = 3.1 MeV を境にして HFB 解は T =1 対相関状態から正常状態に転移することがわかる2 。プラス印 (＋) は T =0 対と T =1 対の両方に対相関力が働くようにした場合に得られた解である。ωτ = 2.8 MeV での R の値の逆転により、T =0 対相関状態から

T =1 対相関状態への転移が起きることが明確に見てとれる3 。この場合は、正常状態の他に T =0 対相関状態と T =1 対相関状態という２つの状態があり、HFB の解法で変分の対象となっている R の値の逆転にともなって片方の対相関状態からもう片方の対相関状態への不連続転移がおきているわけである。

1

水平な直線になるのは、一粒子スペクトルがアイソスピンに依存しないこととアイソクランキングという外場の与え方の特殊性に起因するものである。
2交差点付近では、転移点の前から解の R 値が散乱し始め、正常相への転移は ωτ = 2.6 MeV ですでに起きてしまっている。これは、HFB 解の解法に採用した反復法が、必ずしも R が最低の解に収束しないからであると考えられる。反復ごとの減衰因子を大きくとるなどすれば、改善できると期待される。
3交差点より前の ωτ = 2.5 MeV ですでに T = 0 対相関状態から T = 1 対相関状態への移行が起きている。この理由もまた、反復解法は必ずしも最低 R の状態に収束するとは限らないためであると思われる。

54

まとめると、本研究では、陽子と中性子それぞれに対して別個に対相関を考えるのでなく、他のアイソスピンの組み方の対相関をも完全に含めて扱うことのできるアイ
ソスピン形式の HFB 法の計算プログラムを作成し、そのプログラムの使用例として、

N = Z = 64 核をアイソ空間で回転させた場合に T = 0 対相関状態から T = 1 対相関状態へ不連続な転移が起きることを見た。
今後、このプログラムを利用し、また改良して、下記のような研究をすることが望
ましいと思う。

1. 主殻の異なる軌道にある核子が対を組むことができるように、相互作用行列要素
がそのような対の散乱については零であるという現在の近似を修正する。このことにより N と Z の差の大きい核でもより強く T = 0 対相関が残存するようになると思われる。

2. このために、相互作用行列要素を実際に積分を実行して計算し、主量子数の違いがどの程度行列要素を小さくするかを調べる。

3. スピン軌道結合ポテンシャルを導入することで、スピン軌道混合がどれくらいT = 0 対相関を変化させるかを調べる。特に興味深いのは ls partner 軌道がそれぞれの核子のフェルミ準位付近にあるとき、T = 0 対相関が N 6= Z 核でも強く残る可能性があることである。一方、N = Z 核では、T = 0 対相関が弱まる可能性がある。

4. 変形や l2 ポテンシャルを抜くことで、縮退が対相関のモードにどのような効果があるかを調べる。

55

謝辞
本研究を行うにあたり、御指導を頂いた田嶋直樹先生、鈴木敏男先生、林明久先生に厚く御礼申し上げます。特に田嶋先生には最初から最後まで懇切丁寧な御指導を賜わり感謝の念に堪えません。また、研究を進めるにあたって長谷川健二先生、高木丈夫先生には多くの助言を頂き、有難うございました。
最後になりましたが、応用物理学科の諸先生方に対しましても、日頃の御指導に感謝し、謝辞の言葉と致します。

2002 年 2 月
鈴木 紀史

56

参考文献
[1] J. Bardeen, L.N. Cooper, and J.R. Schrieffer, Phys. Rev. 108, 1175 (1957).
[2] A.L. Goodman, adv. in nucl. phys. 11, 263 (1979).
[3] W.D. Myers and W.J. Swiatecki, Nucl. Phys. 81, 1 (1966).
[4] W. Satula and R. Wyss, Phys. Lett. B393, 1 (1997).
[5] W. Satula and R. Wyss, Phys. Rev. Lett. 86, 4488 (2001).
[6] W. Satula and R. Wyss, Phys. Rev. Lett. 87, 052504 (2001).
[7] 有馬朗人、原子と原子核 (朝倉書店, 東京,1982).
[8] 高木修二 他、岩波講座 現代物理学の基礎 (第 2 版)9 巻 (岩波書店, 東京,1993).
[9] 高木修二, 丸森寿夫、「原子核論」(岩波書店, 東京,1978).
[10] T. Bengtsson, I. Ragnarsson and S. Åberg, “The Carnked Nilsson Model”,
in “Computational Nuclear Physics 1”, ed. K. Langanke, J.A. Maruhn and
S.E. Koonin, (Springer-Verlag, Berlin, 1991) 51.
[11] P. Ring and P. Schuck, The nuclear many-body problem (Springer, New York,1980).
[12] 山本哲朗,「数値解析入門」, p.9, サイエンス社 (1976).
[13] 市村宗武, 坂田文彦, 松柳研一、「原子核の理論」(岩波書店, 東京,1993).
[14] 奥村晴彦、「C 言語による最新アルゴリズム事典」, 技術評論社 (1991).

57

付録

Program List

HFB 方程式を Jacobi 法と反復代入法を用いて解くプログラム/*
isopair.c (ver. 2002/02/11b)
A Hartree-Fock-Bogoliubov program to describe the pairing correlation
taking intofull accout the isospin degree of freedom, i.e. T=0 as well as T=1
pairing,or, p-n as well as p-p and n-n pairing,
written by Norifumi Suzuki and Naoki Tajima, Fukui University【補注】-*- エネルギーの単位には一貫して MeV を想定して作成。
-*- コメント中で A^T は 行列 A の転置行列を表す。
-*- 相互作用チャンネルの説明では T はアイソスピン, S はスピン
L は軌道角運動量の、２核子に関する合成値を表す。
*/#include "matutil.c" /* 行列用小道具集 (奥村晴彦著「アルゴリズム事典」
より)*/#include "lib_jacobi2.c" /* 実対称行列の対角化 (奥村晴彦著同書のものを改変)*/#include <stdio.h>#include <math.h>#include <string.h>#include <stdlib.h>#include <limits.h>#include <sys/time.h> /* Linux での計算時間の計測のために使用 */#include <time.h> /* Linux での計算時間の計測 のために使用 */#define PL 0 /* print level */#define DBG 0 /* debug level */#define I_max 400 /* 取り入れる 1 粒子軌道の最大個数 */const int N_min=0;
値 */const int N_max=5;
値 */const int Max_count=1000;/* 取り入れる 1 粒子軌道の主量子数の最小/* 取り入れる 1 粒子軌道の主量子数の最大/* P1 での自己無動着場を求めるための反復回数の上限 */const double hbar_omega=7.0; /* 調和振動子主殻間のエネルギー間隔 h-bar ω */const double kappa_mu=0.025; /* Nilsson potential parameter κ・μ */const double beta=0.22;/* 軸対称四重極変形度 β_2 */const double lambda_s=40.0; /* アイソクランキングしないときのフェルミ
58

準位 */
const double iomg_beg=0.0;
const double iomg_end=10.0;
const double iomg_inc=0.1;
const double V010=-1.00;

const double V100=-0.66;/* isospace cranking angular velocity *//* _beg 開始値 (begin), _end 終了値, *//* _inc 増分 (increment) *//* T=0,S=1,L=0 channel 相互作用強度(仮標準値 -1.00) *//* T=1,S=0,L=0 channel 相互作用強度(仮標準値 -0.66) *//* Linux での計算時間計測のための変数 : begin */#define NTMPU 20 /* number of time-measured program units */struct tms tbuff;
double cputime[NTMPU];
int calls[NTMPU];
enum{fun1,fun2,fun3,fun7,fun4,fun5,fun6,fun9};/* Linux での計算時間計測のための変数 : end *//* splabel 一粒子状態の量子数をまとめて保持する構造体の型名 */typedef struct{
int n; /* 調和振動子の主量子数 */int l; /* 軌道角運動量量子数 */int m; /* 磁気量子数 (軌道角運動量の z 成分) */int s; /* スピンの z 成分、 s=0 -> -1/2, s=1 -> 1/2 */int t; /* アイソスピンの z 成分 t=0 -> neutron, t=1 -> proton */double e; /* energy */} splabel; /* single particle state label */enum {NEUTRON=0, PROTON=1, DOWN=0, UP=1}; /* splabel の定義はここを参照 *//* ２核子のスピン・アイソスピンの１重項・３重項への射影演算子の行列表現 *//* Proj_...2[s1][s2][s3][s4] = 2 <s1,s2|Projector|s3,s4> *//* s=-1/2 : array index=0,
s=1/2 : array index=1 */int Proj_singlet2[2][2][2][2] = {0,0,0,0, 0,1,-1,0, 0,-1,1,0, 0,0,0,0};
int Proj_triplet2[2][2][2][2] = {2,0,0,0, 0,1, 1,0, 0, 1,1,0, 0,0,0,2};/* 関数プロトタイプ宣言 : begin */int P1(int imax, splabel label[],matrix dnsm, matrix prts, matrix mfhc,
matrix prpt, matrix eigu, matrix eigv ,matrix qpha, matrix qpha1,matrix qpve, vector qpeg , double *num_n_ptr, double *num_p_ptr);
void P2(int imax, matrix mfhc, matrix prpt, matrix qpha);
void P3(int imax, matrix qpha, matrix qpve, vector qpeg, matrix eigu,
matrix eigv);
double P4(int imax, splabel label[], matrix eigu, matrix eigv, matrix dnsm,
matrix prts, double *num_n_ptr, double *num_p_ptr);
int GeneSpBas(int Nmin, int Nmax, int Imax, splabel label[]);void MakeMfhc(double lambda_n, double lambda_p, int imax, splabel label[],matrix mfhc);
void InitPrts(double lambda_n, double lambda_p, int imax, splabel label[],matrix prts);
void CalcPairPot(int imax, splabel label[], matrix prts, matrix prpt,
int channel_t);void CalcEner(int imax, splabel label[], matrix dnsm, matrix prts,
matrix prpt, double *emf_ptr, double *epair_ptr);/* 関数プロトタイプ宣言 : end */59

/****************************************************************************
int main(int argc , char **argv){
int ii,count,loop_iomg=0;
int imax; /* the number of single-particle states */int channel_t; /* interaction is restricted to only T=0 (if channel_t=0),
to only T=1 (if =1), to both T=0 and 1 (if =2) */double num_n,num_p; /* expectation values of the number of neutrons andprotons, <N_n> and <N_p> */double emf,epair,etot; /* Mean-field energy,Pairing energy,total energy */double emf2[3],epair2[3]; /* Mean-field energy and Pairing energy */double iomg;/* λ n = λ s + ω/2, λ p = λ s - ω/2
ω = isospace cranking angular velocity */double lambda_n,lambda_p; /* Fermi levels of neutrons and protons */matrix dnsm,prts,mfhc,prpt,eigu,eigv,qpha,qpha1,qpve; /* 後述の説
明を参照 */
vector qpeg; /* 後述の説明を参照 */splabel label[I_max]; /* 一粒子状態の量子数と準位を記憶しておく配列 */double tm1,tm2; /* variables for CPU time measurement *//*【行列変数の説明】 qp: quasi-particle states,
sp: single-particle basis states
eigu[positive energy qp index][sp index] :Bogoliubov transformation matrix U
eigv[positive energy qp index][sp index] :Bogoliubov transformation matrix V
dnsm[sp index][sp index] : density matrix 密度行列, ρ= V^T V
prts[sp index][sp index] : pairing tensor, t = V^T U
mfhc[sp index][sp index] :mean-field hamiltonian with temrs for number constraint, {curly H} - c
prpt[sp index][sp index] :pairing potential, Δ_{ij}=(1/2) Σ_{kl}<ij|v|kl>t_{kl}qpha[doubled sp index][doubled sp index] :quasi-particle hamiltonian,or HFB super matrix
qpha1[doubled sp index][doubled sp index]:
memory space allocated to preserve previous-step values of qpha
qpve[positive and negative energy qp index][doubled sp index] :quasi-particle vector (eigen vectors of the HFB super matrix)【ベクトル変数の説明】
qpeg[positive and negative energy qp index] :quasi-particle energies, i.e., eigenvalues of qpha*/times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;/* initialization for CPU time measurement : begin */for(ii=0;ii<NTMPU;ii++){
cputime[ii]=0.0;
calls[ii]=0;}/* initialization for CPU time measurement : end */imax = GeneSpBas(N_min, N_max, I_max, label);/* 1 粒子基底の生成, imax:状態の個数 */if(PL>-1){
60

fputs("Program = isopair.c (ver.2002/02/11a)\n",stdout);
printf("Osc.Quant.Num. %d<=N<=%d #single-particle states=%d\n",N_min,N_max,imax);
printf("Mean field: hbar omega=%6.3f kappa_mu=%6.4f beta=%6.4f
lambda_s=%6.3f\n",hbar_omega,kappa_mu,beta,lambda_s);
printf("Pairing interaction strengths V(TSL):
V010=%f V100=%f\n",V010,V100);printf("Isospace cranking: do iso.omega=%f,%f,%f\n",iomg_beg,iomg_end,iomg_inc); fflush(stdout); }if(PL>-1){ /* 一粒子状態を特徴付ける量子数とエネルギー準位の表示 */printf("[ single-paticle states ]\n");printf("%4s %2s %2s %3s %1s %1s %8s\n","No.","n","l","m","s","t","energy");for(ii=0;ii<imax;ii++) printf("%4d %2d %2d %3d %1d %1d %8.4f\n",ii,label[ii].n,label[ii].l,label[ii].m,label[ii].s,label[ii].t,label[ii].e)
fflush(stdout);}/* 関数 malloc の呼び出しによる行列とベクトルの領域確保 */dnsm=new_matrix(imax, imax);
prts=new_matrix(imax, imax);
mfhc=new_matrix(imax, imax);
prpt=new_matrix(imax, imax);
eigu=new_matrix(imax, imax);
eigv=new_matrix(imax, imax);qpha=new_matrix(2*imax, 2*imax);
qpha1=new_matrix(2*imax, 2*imax);
qpve=new_matrix(2*imax, 2*imax);
qpeg=new_vector(2*imax);/* The outer-most loop begins here as for the iso-spacecranking frequency */for(iomg=iomg_beg;
iomg_inc>0 && iomg<iomg_end+1.0e-9 ||iomg_inc<0 && iomg>iomg_end-1.0e-9;
iomg+=iomg_inc){loop_iomg++;lambda_n=lambda_s+iomg/2.0; /* 中性子のフェルミ準位 */lambda_p=lambda_s-iomg/2.0; /* 陽子のフェルミ準位 */if(PL>0) printf("===== iso.omega=%f,iomg,lambda_n,lambda_p);

Fermi levels n=%f

p=%f =====\n"fflush(stdout); fflush(stderr);/* mfhc を作製する。 mfhc は次の iomg までここで作成した値を保つ */MakeMfhc(lambda_n, lambda_p, imax, label, mfhc);/* prts の初期状態を設定する。prts は P1 内のループで自己無撞着性を満たすものに書き換えられる。prpt は prts が変更されるたびに計算しなおされる*/
InitPrts(lambda_n, lambda_p, imax, label, prts);61

/* HFB 方程式の自己無撞着解を求める。count は自己無撞着達成までの反復回数 */
count=P1(imax, label, dnsm, prts, mfhc, prpt, eigu, eigv, qpha, qpha1,
qpve ,qpeg, &num_n, &num_p);if(PL>1){ /* 必要ならば、ここで、一粒子基底状態の占拠確率を表示する */for(ii=0;ii<imax;ii++){printf("%f %f %f %d :occ prb\n",dnsm[ii][ii],label[ii].e,
mfhc[ii][ii],label[ii].t);}}/* 解が求まったあと、epair を T=0 の寄与と T=1 の寄与に分けて計算する */for(channel_t=2;channel_t>=0;channel_t--){/*T=0,T=1 の場合分け*//* 解として得られた対相関テンソル prts をそのまま用い、相互作用を変えて対相関ポテンシャル prpt を計算する */CalcPairPot(imax, label, prts, prpt, channel_t);/* 変更された対相関ポテンシャル prpt を用いて対相関エネルギー epair を計算する */CalcEner(imax, label, dnsm, prts, prpt, &emf2[channel_t],&epair2[channel_t]);}emf=emf2[2]; epair=epair2[2]; etot=emf+epair;/* 平均場エネルギーが変わったらプログラムにバグがある */if(fabs(emf2[2]-emf2[1])>1.0e-8||fabs(emf2[2]-emf2[0])>1.0e-8)
fprintf(stderr,"Error: it must hold emf= %f = %f= %f\n",emf2[0],emf2[1],emf2[2]);/* 対相関エネルギーが各チャンネルの対相関エネルギーの和に一致しなければバグがある */if(fabs(epair2[2]-epair2[0]-epair2[1]) > 1.0e-8)
fprintf(stderr,"Error: it must hold epair %f = %f + %f = %f\n",epair2[2],epair2[0],epair2[1],epair2[0]+epair2[1]);/* 主要出力 */if(PL > -1) {printf("%f %e %e %e %f %f %d iomg,emf,epair(T=0,1),N,Z,iter\n",iomg , emf2[2], epair2[0], epair2[1], num_n, num_p ,count );}} /* end of the outer-most loop for the iso-space cranking frequency
iomg */
times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun1]+=1; cputime[fun1]+=tm2-tm1;/* 計算時間を出力 */printf("total %12.2f sec / %d\n",cputime[fun1],calls[fun1]);62}

printf("GeneSpBas: %e %d %e\n",cputime[fun2],calls[fun2],
cputime[fun2]/calls[fun2]);printf("P3: %e %d %e\n",cputime[fun3],calls[fun3],
cputime[fun3]/calls[fun3]);
printf("P1: %e %d %e\n",cputime[fun7],calls[fun7],
cputime[fun7]/calls[fun7]);printf("P2: %e %d %e\n",cputime[fun5],calls[fun5],
cputime[fun5]/calls[fun5]);
printf("P3: %e %d %e\n",cputime[fun6],calls[fun6],
cputime[fun6]/calls[fun6]);
printf("P4: %e %d %e\n",cputime[fun9],calls[fun9],
cputime[fun9]/calls[fun9]);printf("CalcPairPot: %e %d %e\n",cputime[fun4],calls[fun4],
cputime[fun4]/calls[fun4]);/*******************************************************//* Solve the HFB equation and obtain a self-consistent solution */int P1(int imax, splabel label[] ,matrix dnsm, matrix prts, matrix mfhc,
matrix prpt , matrix eigu, matrix eigv ,matrix qpha, matrix qpha1,matrix qpve, vector qpeg ,double *num_n_ptr, double *num_p_ptr) {int ii,jj,nn,count=0;
double err, errmax, emf, epair, num_n, num_p, qpeg2[2*I_max];
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
while(count<Max_count){
CalcPairPot(imax, label, prts, prpt, 2);/* 対相関ポテンシャルの行列 prpt を作製 */P2(imax, mfhc, prpt, qpha);/* 準粒子ハミルトニアンの行列 qpha を作製 *//* damping of the convergence to suppress oscillations */if(count>0){for(ii=0;ii<2*imax;ii++){for(jj=0;jj<2*imax;jj++){qpha[ii][jj]=0.80*qpha1[ii][jj]+0.20*qpha[ii][jj];}}}for(ii=0;ii<2*imax;ii++){for(jj=0;jj<2*imax;jj++){qpha1[ii][jj]=qpha[ii][jj];}}/* ヤコビ法による対角化で Boguliubov 準粒子解を求める */P3(imax, qpha, qpve, qpeg, eigu, eigv);/* calculation of density matrix, pairing tensor,
expectation values of the numbers of neutrons and protons */P4(imax, label, eigu, eigv, dnsm, prts, &num_n, &num_p);63

if(PL>0){
CalcEner(imax, label, dnsm, prts, prpt, &emf, &epair);printf("%5d %e %e %f %f\n",count+1,emf,epair, num_n, num_p);}if(count>0) {errmax=0.0;
for(ii=0;ii<2*imax;ii++){err=fabs(qpeg2[ii]-qpeg[ii]);
if(err > errmax) errmax=err;}
if(PL>0) printf("iter=%d max change of eigenvalues=%e\n",count+1,errmax);
if(errmax < 1.0e-4) break;}for(ii=0;ii<2*imax;ii++) qpeg2[ii]=qpeg[ii];/* 収束判定のために保存 */count++;}*num_n_ptr=num_n; *num_p_ptr=num_p;

times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun7]+=1; cputime[fun7]+=tm2-tm1;}

return ++count;/*******************************************************//* 準粒子ハミルトニアンの行列 qpha を作製 */void P2(int imax, matrix mfhc, matrix prpt, matrix qpha){
int ii,ii2,jj,jj2;
double err,errmax,dgmin,dgmax;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
for(ii=0,ii2=imax;ii<imax;ii++,ii2++){for(jj=0,jj2=imax;jj<imax;jj++,jj2++){qpha[ii][jj] = mfhc[ii][jj];qpha[ii2][jj2]=-mfhc[ii][jj];qpha[ii][jj2] = prpt[ii][jj];qpha[ii2][jj] =-prpt[ii][jj];}}if(DBG>9) {/* check of the hermiticity of the created matrix */errmax=0;
for(ii=0;ii<imax;ii++){for(jj=0;jj<imax;jj++){err=fabs(qpha[ii][jj]-qpha[jj][ii]);
if(err>errmax) errmax=err;}}
64

if(PL>5||errmax > 1.0e-12)
fprintf(stderr,"error of hermiticity of q.p. hamiltonian= %e\n",errmax);} /* end of if DBG>9 */if(PL>8) {/* show the range of the diagonal elements of the created matrix*/
dgmin=dgmax=qpha[0][0];for(ii=1;ii<imax;ii++){if(qpha[ii][ii] > dgmax) dgmax=qpha[ii][ii];if(qpha[ii][ii] < dgmin) dgmin=qpha[ii][ii];}printf("P2: upper-block diagonal elements are in [%e %e]\n",dgmin,dgmax);
dgmin=dgmax=qpha[imax][imax];
for(ii=imax+1;ii<2*imax;ii++){if(qpha[ii][ii] > dgmax) dgmax=qpha[ii][ii];if(qpha[ii][ii] < dgmin) dgmin=qpha[ii][ii];}printf("P2: lower-block diagonal elements are in [%e %e]\n",dgmin,dgmax);} /* end of if PL>8 */}times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun5]+=1; cputime[fun5]+=tm2-tm1;/*******************************************************//* 準粒子ハミルトニアンの対角化 */void P3(int imax, matrix qpha, matrix qpve, vector qpeg, matrix eigu,
matrix eigv){
int ii,ii2,jj,jj2,kk;
double err,errmax,s1,s2,s3,s4;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;jacobi(2*imax, qpha, qpve); /* ヤコビ法による対角化 */for(ii=0;ii<2*imax;ii++){ /* 固有値を配列にコピーする */qpeg[ii]=qpha[ii][ii];if(PL>5) printf("qp energy[%4d]= %10.5f\n",ii,qpeg[ii]);}/* 正エネルギー準粒子解の固有ベクトルを Bogolliubov 変換行列 U,Vにコピーする */for(ii=0;ii<imax;ii++){ /* ii: quasi-particle index */for(jj=0,jj2=imax;jj<imax;jj++,jj2++){/* jj,jj2-imax:single-particle index */eigu[ii][jj]=qpve[ii][jj];eigv[ii][jj]=qpve[ii][jj2];}}/* 固有値のチェック : begin */if(DBG>-1){
65

errmax=0.0;
for(ii=0,ii2=2*imax-1;ii<imax;ii++,ii2--){if(qpeg[ii] <-1.0e-6) fprintf(stderr,"P3:error:upper blocknegative eigen_value %d %f\n",ii,qpeg[ii]);
if(qpeg[ii2]> 1.0e-6) fprintf(stderr,"P3:error:lower block
positive eigen_value %d %f\n",ii2,qpeg[ii2]);
err=fabs(qpeg[ii]+qpeg[ii2]);
if(err>0.1){
fprintf(stderr,"Eigenvalues are NOT pairwise %d %f %d %f\n",ii,qpeg[ii],ii2,qpeg[ii2]);}
if(err>errmax) errmax=err;}
if(PL>1||errmax>1.0e-2){
printf("Max error of eigenvalue pairwiseness = %e\n",errmax);
if(errmax>1.0e-2) fputs("* Large Error *",stdout);
fputs("\n",stdout);}}/* 固有値のチェック : end *//* check of the orthonomality of U and V : begin *//* U U^T + V V^T = I, U^T U + V^T V = I,
U V^T + V U^T = 0, U^T V + V^T U = 0 */if(DBG>1){
errmax=0;
for(ii=0;ii<imax;ii++){for(jj=0;jj<imax;jj++){if(ii==jj)s1=-1;else s1=0;
s2=s1;s3=0;s4=0;
for(kk=0;kk<imax;kk++){s1+=eigu[ii][kk]*eigu[jj][kk]+eigv[ii][kk]*eigv[jj][kk];s2+=eigu[kk][ii]*eigu[kk][jj]+eigv[kk][ii]*eigv[kk][jj];s3+=eigu[ii][kk]*eigv[jj][kk]+eigv[ii][kk]*eigu[jj][kk];s4+=eigu[kk][ii]*eigv[kk][jj]+eigv[kk][ii]*eigu[kk][jj];}err=fabs(s1); if(err>errmax)errmax=err;
err=fabs(s2); if(err>errmax)errmax=err;
err=fabs(s3); if(err>errmax)errmax=err;
err=fabs(s4); if(err>errmax)errmax=err;}}
if(PL>0||errmax>1.0e-6){
printf("U,V orthonormality error= %e %e %e %e max=%e",s1,s2,s3,s4,errmax);
if(errmax>1.0e-6) fputs("* Large Error *",stdout);
fputs("\n",stdout);}}/* check of the orthonomality of U and V : end */}times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun6]+=1; cputime[fun6]+=tm2-tm1;/*******************************************************/
66

/*密度行列 ρ, 対相関テンソル t, 粒子数期待値 Np, Nn の計算 */double P4(int imax, splabel label[], matrix eigu, matrix eigv,
matrix dnsm, matrix prts, double *num_n_ptr, double *num_p_ptr){
int ii,jj,kk;
double s,num_n,num_p;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;/* Calculation of density matrix dnsm = V^T V */for(ii=0;ii<imax;ii++){for(jj=0;jj<=ii;jj++){s=0;
for(kk=0;kk<imax;kk++) s+=eigv[kk][ii]*eigv[kk][jj];dnsm[ii][jj]=s;if(ii != jj) dnsm[jj][ii]=s; /* because dnsm^T = dnsm */}}/* Calculation of pairing tensor t = V^T U */for(ii=0;ii<imax;ii++){prts[ii][ii]=0;
for(jj=0;jj<=ii;jj++){for(s=0.0,kk=0;kk<imax;kk++) s+=eigv[kk][ii]*eigu[kk][jj];prts[ii][jj]= s;prts[jj][ii]=-s; /* because t^T = -t */}}/* Calculation of expectation values of proton and neutron numbers */for(num_n=0.0, num_p=0.0, ii=0;ii<imax;ii++){if(label[ii].t==NEUTRON) num_n+=dnsm[ii][ii]; else num_p+=dnsm[ii][ii];}*num_n_ptr=num_n; *num_p_ptr=num_p;}

times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun9]+=1; cputime[fun9]+=tm2-tm1;/*******************************************************//* 一粒子状態を生成し量子数とエネルギースペクトルを配列に記憶 *//* Generate single-particle basis */int GeneSpBas(int Nmin, int Nmax, int Imax, splabel label[]){int n,l,m,s,t,i=0,absm,dm;
double e;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
if(Nmin<0||Nmax<Nmin||Imax<1){fprintf(stderr,"GeneSpBas:error:%d %d %d\n",Nmin,Nmax,Imax);exit(1);}for(n=Nmin;n<=Nmax;n++){for(l=n;l>=0;l-=2){
for(absm=0;absm<=l;absm++){dm=absm*2; if(dm==0) dm=1;
for(m=-absm;m<=absm;m+=dm){
for(s=0;s<=1;s++){
67

}

}

}

}

}

for(t=0;t<=1;t++){if(i>Imax){
fprintf(stderr,"GeneSpBas:error i>I_max=%d\n",Imax);exit(1);}e=hbar_omega*((n+1.5)-kappa_mu*(l*(l+1.0)-((n*(n+3.0))/2.0))-beta*sqrt(5.0/(4.0*M_PI))*(l*(l+1)-3*m*m)*(n+1.5)/((2*l-1)*(2*l+3)));label[i].n=n; label[i].l=l; label[i].m=m;label[i].s=s; label[i].t=t; label[i].e=e;
i++;
}

}
times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun2]+=1; cputime[fun2]+=tm2-tm1;
return i; /* 戻り値は 状態数 i, 状態番号は 0,...,i-1 *//*******************************************************//* prepare the mean-field hamiltonian (mfhc) */void MakeMfhc(double lambda_n, double lambda_p, int imax,
splabel label[], matrix mfhc){
int ii,jj;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
for(ii=0;ii<imax;ii++){for(jj=0;jj<imax;jj++) mfhc[ii][jj]=0.0;
mfhc[ii][ii]=label[ii].e;
if(label[ii].t==NEUTRON){ mfhc[ii][ii]-=lambda_n;}else{ mfhc[ii][ii]-=lambda_p;}}}times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun3]+=1; cputime[fun3]+=tm2-tm1;/*******************************************************//* Give initial values to the pairing tensor (prts) */void InitPrts(double lambda_n, double lambda_p, int imax,
splabel label[], matrix prts){
int ii,jj;
double tm1,tm2;
times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
for(ii=0;ii<imax;ii++){prts[ii][ii]=0;
for(jj=0;jj<ii;jj++){prts[ii][jj]=0.5+0.1*(rand()/((double) RAND_MAX));prts[jj][ii]=-prts[ii][jj]; /* prts is an antisymmetric matrix */
68

}

}

}

times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun3]+=1; cputime[fun3]+=tm2-tm1;/*******************************************************//* 平均場 (Hartree-Fock) エネルギーと対相関エネルギーの計算 */void CalcEner(int imax, splabel label[], matrix dnsm, matrix prts,
matrix prpt, double *emf_ptr, double *epair_ptr){
int ii,jj;
double emf,epair;for(emf=0.0, epair=0.0, ii=0;ii<imax;ii++){emf+=(label[ii].e)*dnsm[ii][ii];for(jj=0;jj<ii;jj++){epair+=prpt[ii][jj]*prts[ii][jj];}}*emf_ptr=emf; *epair_ptr=epair;}/*******************************************************//* Calculate the pairing potential (prpt) using the pairing tensor(prts) and the interaction. Depending on the control parametervalue channel_t=0,1,2,respectively, only T=0, only T=1 ,and both T=0 and T=1 interactions are taken into account.*/
void CalcPairPot(int imax, splabel label[], matrix prts, matrix prpt,
int channel_t){
int ii,jj,kk,ll;
int n1,n2,n3,n4, l1,l2,l3,l4, m1,m2,m3,m4, s1,s2,s3,s4, t1,t2,t3,t4;
double v,v0,v1,vd,sum;double tm1,tm2;
if(fabs(V100)<1.0e-14 && fabs(V010)<1.0e-14){for(ii=0;ii<imax;ii++) for(jj=0;jj<imax;jj++) prpt[ii][jj]=0.0;
return;}times (& tbuff); tm1=(double) tbuff.tms_utime/(double) CLK_TCK;
for(ii=0;ii<imax;ii++){prpt[ii][ii]=0; /* pairing potential is an antisymmetric matrix */n1=label[ii].n;l1=label[ii].l;m1=label[ii].m;s1=label[ii].s;t1=label[ii].t;for(jj=0;jj<ii;jj++){n2=label[jj].n; if(n1 != n2) continue;
l2=label[jj].l; if(l1 != l2) continue;
m2=label[jj].m; if(m1 !=-m2) continue;
s2=label[jj].s;t2=label[jj].t;
sum=0.0;for(kk=0;kk<imax;kk++){n3=label[kk].n;l3=label[kk].l;m3=label[kk].m;s3=label[kk].s;t3=label[kk].t;69

vd=(((l1-m1+l3-m3)%2) ? -1:1)/(sqrt((2.0*l1+1.0)*(2.0*l3+1.0)));for(ll=0;ll<kk;ll++){n4=label[ll].n; if(n3 != n4) continue;
l4=label[ll].l; if(l3 != l4) continue;
m4=label[ll].m; if(m3 !=-m4) continue;
s4=label[ll].s;t4=label[ll].t;v0=V010*(Proj_singlet2[t1][t2][t3][t4]*Proj_triplet2[s1][s2][s3][s4]-Proj_singlet2[t1][t2][t4][t3]*Proj_triplet2[s1][s2][s4][s3]);
v1=V100*(Proj_triplet2[t1][t2][t3][t4]*Proj_singlet2[s1][s2][s3][s4]-Proj_triplet2[t1][t2][t4][t3]*Proj_singlet2[s1][s2][s4][s3]);if(channel_t==2) v=v0+v1;
else if (channel_t==1) v=v1;
else v=v0;
sum+=v*vd*prts[kk][ll];}}prpt[ii][jj]=sum*0.25;
prpt[jj][ii]=-prpt[ii][jj];}}}times (& tbuff); tm2=(double) tbuff.tms_utime/(double) CLK_TCK;
calls[fun4]+=1; cputime[fun4]+=tm2-tm1;/*******************************************************/

70