2019 年度計算物理学２第 3 回レポート課題
原子核の基底状態の中には核子間の対相関によって超伝導状態になっているものが多くある。原子核内の中性子の超伝導状態をギャップ方程式を解いて求めてみよう。ギャップ方程式は∆=Ων ∆G∑√
2 ν(eν − λ)2 + ∆2(1)で与えられる。ここで ∆ は対ギャップ、G は超伝導を引き起こす対相互作用の強度、λ は化学ポテンシャル、ν は原子核内の核子が占有する一粒子軌道をラベルする添字で、一粒子軌道のエネルギーは eν でそれぞれの一粒子軌道は 2Ων 重に縮退している。∆ = 0 はギャップ方程式の自明解であるが、系が超伝導状態にあるときは ∆ ̸= 0 の解が存在しうる。
原子核内の中性子のそれぞれの一粒子軌道の占有率は超伝導を記述する BCS 理論で1[Uν (λ, ∆)] =22(1+ √)eν − λ1[Vν (λ, ∆)] =22,(eν − λ)2 + ∆2()eν − λ1− √(eν − λ)2 + ∆2,(2)と与えられる。ここで [Vν (λ, ∆)]2 が軌道の占有確率で、[Uν (λ, ∆)]2 は軌道の非占有の確率に対応しており、各一粒子軌道で [Uν (λ, ∆)]2 + [Vν (λ, ∆)]2 = 1 である (和が 1 で確率と解釈できる)。これらの量は化学ポテンシャル λ と対ギャップ ∆ の関数として与えられ、∆ は∆≡G∑Ων Uν (λ, ∆)Vν (λ, ∆)(3)νで定義され、これに式 (2) を代入するとギャップ方程式 (式 (1)) が導出できることが確かめられる。
一粒子軌道に入っている中性子の数の総和 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 個まで入れることができる。対相互作用 G がない場合や非常に弱い場合は中性子はフェルミ粒子であるため最低エネルギー状態 (ν = 1) の 4 つを占有する。これは∆ = 0 の自明解に対応する。(∆ = 0 とし、e1 < λ < e2 と選ぶことで V1 = 1, V2 = V3 = 0 とできる)
問 1. まずは対ギャップの値が固定値 ∆(0) = 1 MeV で与えられている場合を考える (∆(0) は式 (1) や (3)を満たさないがこれらの式はとりあえず問題 1. と 2. では無視する。) このとき、この３つの一粒子軌道にある中性子数 N は式 (4) と (2) より

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) ) も合わせて書け。

1

問 2. 前問より N (λ, ∆(0) ) は λ に関して単調増加関数であることがわかった。λ の値はパラメータとして N (λ, ∆(0) ) の値が計算したい原子核の中性子数に対応するように決定する。∆(0) = 1 MeV のときにN (λ, ∆(0) ) = 4 となる λ の値を二分法で求めよ。
問 3. 前問で ∆(0) = 1 MeV に対して N (λ, ∆(0) ) が原子核 60 Ni に対応するように λ を決定することができた。しかし本来は ∆ はギャップ方程式を満たすように決める必要がある。∆ は式 (3) より Uν と Vν から計算できるが、Uν や Vν は ∆ に依存する。このようにギャップ方程式は ∆ に関して非線形な方程式となっている(このような求めたい解 ∆ が自分自身に依存する方程式を自己無撞着方程式という)。
自己無撞着方程式は反復法を用いて数値解を求めることができる。以下の反復法のアルゴリズムを用いて、自己無撞着方程式 (2,3,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) ) を式 (2) を用いて計算する。
∑ (k)(k)e) 式 (3) より ∆(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 ループに名前をつけてもよい。(Fortran 文法のまとめ参照。EXIT 文で DO ループを抜けるときも名前を指定できる。)• 問 4. で G を変えて ∆ が小さくなっていくと N (λ, ∆) の λ の関数としての振る舞いが階段関数に近く2なっていくため、もし N (λ, ∆) = 4 の解を Newton 法で求めていると解けなくなる。二分法であれば非常に小さな ∆ に対しても安定して解が求められる。

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

3