計算物理学２第４回レポート課題
原子核の基底状態の中には核子間の対相関によって超伝導状態になっているものが多くある。原子核内の中性子の超伝導状態を考える。原子核内の中性子は 2Ων 重に縮退しているエネルギーが eν の一粒子軌道を占有しているが (ν は一粒子軌道をラベルする添字)、その占有率は超伝導を記述する BCS 理論で1[Uν (λ, ∆)] =22(1+ √)eν − λ

1[Vν (λ, ∆)] =22,(eν − λ)2 + ∆2()eν − λ1− √(eν − λ)2 + ∆2,(1)と与えられる。ここで [Vν (λ, ∆)]2 が軌道の占有確率で、[Uν (λ, ∆)]2 は軌道の非占有の確率に対応しており、各一粒子軌道で [Uν (λ, ∆)]2 + [Vν (λ, ∆)]2 = 1 である (和が 1 なので確率と解釈できる)。これらの量は化学ポテンシャル λ と対ギャップ ∆ の関数として与えられ、∆ は対相互作用の強度 G を用いて∆≡G∑Ων Uν (λ, ∆)Vν (λ, ∆)(2)νで定義され、これに式 (1) を代入すると以下のギャップ方程式を満たすことが確かめられる。

∆=Ων ∆G∑√
2 ν(eν − λ)2 + ∆2(3)∆ = 0 はギャップ方程式の自明解であるが、系が超伝導状態にあるときは ∆ ̸= 0 の解が存在しうる。
一粒子軌道に入っている中性子の数の総和 N は占有確率と縮退度の積の和で与えられる。

N (λ, ∆) = 2∑

Ων [Vν (λ, ∆)]2(4)ν

さて 60 Ni 原子核の中性子の超伝導性を調べてみる。この原子核は 32 個の中性子で構成されているが、そのうち 28 個は閉殻構造をとって安定化しており超伝導性に寄与しない。残りの 4 個が超伝導性に寄与する。３つの一粒子軌道 ν = 1, 2, 3 を考え、4 つの中性子がこの軌道に存在するとする。(上の式での ν に関する和は3∑∑νに置き換える) 一粒子軌道のエネルギーは e1 = 0.0 MeV, e2 = 2.3 MeV, e3 = 2.8 MeV で与えられる

ν=1とする。また、軌道の縮退度 (2Ων ) は Ω1 = 2, Ω2 = 3, Ω3 = 1 とする。つまり、3 つの一粒子軌道に中性子をそれぞれ下から 4 つ、6 つ、2 つ、計 12 個まで入れることができる。相互作用がない場合は中性子はフェルミ粒子であるため最低エネルギー状態 (ν = 1)4 つを占有する。これは ∆ = 0 の自明解に対応する。

1. まずは対ギャップの値が固定値 ∆(0) = 1 MeV で与えられている場合を考える (∆ は式 (2) や (3) を満たす必要があるがこれらの式はとりあえず問題 1. と 2. では無視する。) このとき、この３つの一粒子軌道にある中性子数 N は式 (4) と (1) よりN (λ, ∆(0))=23∑Ων [Vν (λ, ∆(0)2)] =ν=13∑
ν=1(Ωνeν − λ1− √(eν − λ)2 + (∆(0) )2)(5)で与えられる。N (λ, ∆(0) ) は ∆(0) が固定値であれば λ のみの関数である。N (λ, ∆(0) ) を λ の関数としてグラフをかけ。N が 0 から 12 程度の領域が含まれるように λ の範囲を調整すること。また ∆(0) = 0.01 MeVのとき N (λ, ∆(0) ) はどのようなグラフになるか。

2. 前問より N (λ, ∆(0) ) は λ に関して単調増加関数であることがわかった。λ の値はパラメータとして N (λ, ∆(0) ) の値が計算したい原子核の中性子数に対応するように決定する。∆(0) = 1 MeV のときにN (λ, ∆(0) ) = 4 となる λ の値を二分法で求めよ。
1

3. 前問で ∆(0) = 1 MeV に対して N (λ, ∆(0) ) が原子核 60 Ni に対応するように λ を決定することができた。しかし本来は ∆ はギャップ方程式を満たすように決める必要がある。∆ は式 (2) より Uν と Vν から計算できるが、Uν や Vν は ∆ に依存する。このようにギャップ方程式は ∆ に関して非線形な方程式となっている(このような求めたい解 ∆ が自分自身に依存する方程式を自己無撞着方程式という)。
自己無撞着方程式は反復法を用いて数値解を求めることができる。以下の反復法のアルゴリズムを用いて、自己無撞着方程式 (1,2,4) を解き、∆ の値を求めよ。ただし G = 0.6 MeV とする。

a) 対ギャップ ∆ = ∆(0) にゼロでない適当な初期値を与える (前問で使った ∆(0) = 1 MeV でよい)b) 以下の c) から g) を DO ループの中に入れて反復計算を行う。
c) (k 回目の反復計算) 初期値、あるいは前回の反復で求まった ∆ = ∆(k−1) に対して N (λ, ∆(k−1) ) = 4
を二分法で解き、この解 λ = λ(k) を決定する。
(k)(k)d) ∆ = ∆(k−1) , λ = λ(k) のときの Uν (λ(k) , ∆(k−1) ), Vν (λ(k) , ∆(k−1) ) を式 (1) を用いて計算する。
∑ (k)(k)e) 式 (2) より ∆(k) = G ν Uν (λ(k) , ∆(k−1) ), Vν (λ(k) , ∆(k−1) ) を計算する。この式によって前回の対ギャップの値 ∆(k−1) を用いて対ギャップの新しい値 ∆(k) が計算される。

f) ∆ が収束しているか判定する。|∆(k) − ∆(k−1) | が十分に小さければこの DO ループから抜ける。
g) まだ ∆ が収束していない場合、プログラムで ∆ を保存している変数の値を前の値 (∆(k−1) ) から新しいもの (∆(k) ) に更新し、c) に戻って k + 1 回目の反復計算を行う。

h) DO ループはここまで。DO ループの外側で収束解の ∆ を出力
4. 発展問題 (解きたい方のみ) ∆ ̸= 0 となる超伝導解が存在するためには超伝導を引き起こす引力の対相互作用 G > 0 がある程度大きい必要がある。G の値を変え、それぞれの G の値に対して 3. のようにギャップ方程式を解くことで ∆ を G の関数としてプロットすると ∆ ̸= 0 となる解が存在する G の最小の値 Gcrit がある。図より Gcrit のおよその値を示せ。

1. のグラフ、2. の λ の値と 3. を計算するプログラムと ∆ の値、解いた方は 4 のグラフも提出してください。
■ヒント• 3 つの一粒子軌道があるので、eν , Ων , Uν , Vν の対応する量は要素数が 3 つの一次元配列に格納するとよい。

• N (λ, ∆) − 4 を計算する関数副プログラムを作るとこれがゼロになる λ が 2. での解となる。
• 2. を解くときは 1. で書いたグラフから解 λ の位置をおおよそ見積もっておく。N (λ, ∆) は λ に関して単調増加であるので、二分法で解く時に解を含む区間 [λL , λR ] はある程度広く取るとよい (∆ の値を変えても同じ区間から計算をスタートできるため) たとえば N が 0∼1 付近の λ の値を λL ,N = 11 付近の値を λR として初期値に採用しておくとよい。

• 3. は反復法で自己無撞着方程式を解くため DO ループを用いるが、c) の二分法自体も反復法で解くため、二重反復となる。混乱を避けるために DO ループに名前をつけてもよい。(Forran 文法のまとめ 4ページ一番下参照。EXIT 文で DO ループを抜けるときも名前を指定できる。)• 4. で G を変えて ∆ が小さくなっていくと N (λ, ∆) の λ の関数としての振る舞いが階段関数に近くなっていくため、もし N (λ, ∆) = 4 の解を Newton 法で求めていると解けなくなる。二分法であれば非常に小さな ∆ に対しても安定して解が求められる。

• 4. で G があまりに Gcrit に近い値ではギャップ方程式の解が求まりにくいのでその付近は無理に求めなくてよい。

2