令和元年度
卒業論文

アフィン変換数値積分法を用いた超伝導体内の量子化磁束線の可視化

田中宏樹(学籍番号：16232053)

九州工業大学

情報工学部

電子情報工学科
指導教員：小田部荘司

令和 2 年 2 月 20 日

教授

目次
第1章

序論 .......................................................................................................................1

1.1

はじめに ....................................................................................................................1

1.2

Time Dependent Ginzburg-Landau 方程式 .............................................................1

1.2.1

GL 方程式 ..........................................................................................................1

1.2.2

TDGL 方程式 ....................................................................................................2

1.2.3

TDGL 方程式の簡易化 ......................................................................................3

1.3

アフィン変換数値積分法(Affine Integrator(AFI)) ...................................................7

1.3.1

AFI 法のイメージ ..............................................................................................7

1.3.2

1 次元空間において AFI 法の数値積分スキームと時間発展方程式の導出.......8

1.3.3

2 次元空間において AFI 法の時間発展方程式の導出 .....................................13

1.4

研究目的 ..................................................................................................................14

第2章

実装方法と計算条件 ............................................................................................162.1

実装方法 ..................................................................................................................162.2

初期条件 ..................................................................................................................16

2.3

シミュレーション領域内での更新式 ......................................................................17

2.3.1

オーダパラメータの更新式 .............................................................................172.3.2

リンク変数の更新式 ........................................................................................172.4

境界条件 ..................................................................................................................18

2.4.1

オーダパラメータの境界条件 ..........................................................................18

2.4.2

リンク変数の境界条件.....................................................................................18

2.5

シミュレーション領域での更新順序 ......................................................................19

第3章

実行結果および考察 ............................................................................................20

3.1

実行結果 ..................................................................................................................20

3.2

考察および今後の課題 ............................................................................................21

第4章

結言 .....................................................................................................................23

謝辞 ........................................................................................................................................24
参考文献 ................................................................................................................................25i

図目次
図 1-1 空間に関する離散化のイメージ ....................................................................................... 8
図 1-2 時間に関する離散化のイメージ ....................................................................................... 8
図 1-3 1 次元空間での周期的境界条件下の 座標とリンク変数のインデックスの定義 ............. 9
図 1-4 2 次元空間での周期的境界条件下の 座標とリンク変数のダブルインデックスでの定義............................................................................................................................................. 13
図 2-1 シミュレーション領域の オーダパラメータとリンク変数の配置の概念図 .................. 16
図 3-1 作成したプログラムの実行結果(𝐽𝑦 = 0の場合) ............................................................ 20
図 3-2 作成したプログラムの実行結果(𝐽𝑦 = 1の場合) ............................................................ 21
図 3-3 𝐽𝑦 = 0で長時間プログラムを実行した結果 ................................................................... 22ii

第1章 序論
1.1

はじめに

超伝導現象は、超伝導体という物質の電気抵抗が温度の低下とともに消滅する現象であり、1911 年にオランダの物理学者 H. K. Onnes によってはじめて水銀で発見された。その後多くの物質について超伝導現象が確認されている。超伝導体はその電気抵抗が消滅するという性質から様々な工学的な応用に期待されてきた。しかし、その性質はわずかな磁界や温度によって失われてしまう。その転移温度を臨界温度𝑇𝑐 、転移磁界を臨界磁界𝐵𝑐 と呼ぶ。
超伝導現象については、現在までにその発現機構や性質に関する研究が進められており、1933 年にドイツの物理学者 W. Meißner と R. Ochsenfeld によって完全反磁性(MeißnerOchsenfeld 効果)が発見され、さらに 1957 年に J．Bardeen、 L. N. Cooper、 J. R. Schriefferらによって BCS 理論が提唱され、超伝導現象の発現機構における基本的な理解が与えられた。 このとき BCS 理論では、超伝導体の𝑇𝑐 は 30 K を超えない(これは BCS 理論の壁と呼ばれる)と予想されていた。しかし、1986 年にドイツの物理学者 J. G. Bednorz とスイスの物理学者 K. A. Müller らによって𝑇𝑐 が 35 K となる La-Ba-Cu-O 系の超伝導体が発見された。この発見以降、世界各国で𝑇𝑐 の高い超伝導体の探索が行われ、1987 年に窒素の沸点 (77K)より高い𝑇𝑐 をもつ超伝導体が発見された。このことで入手が困難なヘリウムではなく窒素を冷媒として用いることができるようになったため、冷却コストの低減により様々な機器への応用が期待されている。しかしこれらの高温超伝導体にも実用化に向けて様々な課題が残されているため現在でも研究が続けられている。

1.2

Time Dependent Ginzburg-Landau 方程式

こ こ で は 、 超 伝 導 体 の 数 値 解 析 に 用 い る 時 間 依 存 Ginzburg-Landau 方 程 式 (Time
Dependent Ginzburg-Landau 方程式)について説明する。

1.2.1 GL 方程式
Ginzburg-Landau(G-L)理論は、1950 年に V. L. Ginzburg と L. D. Landau によって提唱された超伝導現象を説明する現象論である。G-L 理論は磁界と超伝導が共存する場合の相転移を取り扱ったもので、特に第 2 種超伝導体の磁気特性を良く記述することが知られ1ている。以下では、G-L 理論を成り立たせる仮定とそれによって求められる方程式を説明する。
G-L 理論では、まず超伝導状態の秩序の程度を表す量として複素数であるオーダパラメ
ータ𝛹 = |𝛹|exp(i𝜑)を定義する。そして𝛹は、|𝛹|2 ∝ 𝑛s(1.1)の関係を満たすと仮定する。ここで、𝑛s は超伝導電子密度である。超伝導状態の自由エネル
ギー𝐸s は𝑛s に依存しているから(1.1)式より、|𝛹|2 の関数である。よって、𝐸s は|𝛹|2のべき展
開で表すことが出来る。

𝐸s = 𝐸n + 𝛼|𝛹|2 +

𝛽
|𝛹|4
2

(1.2)

ここで、𝐸n は常伝導状態の自由エネルギー、𝛼および𝛽はそれぞれべき展開した際の 1 次と2 次の係数であり、𝑇 < 𝑇c では𝛼 < 0, 𝛽 > 0である。2 次の項までの展開であるのは|𝛹|2が転移点近傍(𝑇 ≅ 𝑇c )において十分に小さいためである。
次に、磁界の存在が𝛹の空間的変化に寄与することを考慮し、(1.2)式に磁界のエネルギー
密度と運動エネルギー密度を加算する。

𝐸s = 𝐸n + 𝛼|𝛹|2 +𝛽1
1|𝛹|4 +(𝜵 × 𝑨)2 +|(−iℏ𝜵 + 𝑒 ∗ 𝑨)𝛹|2
2
2𝜇0
2𝑚∗

(1.3)

ここで、𝜇0 は真空中の透磁率、𝑨はベクトルポテンシャル、𝑚∗は超伝導電子の質量、ℏはプ
ランク定数を2𝜋で除算したもの、𝑒 ∗は超伝導電子の電荷量、i は虚数単位である。
𝛹と𝑨は(1.3)式の𝐸s を最小とするように決定される。よって𝛹の共役複素数𝛹 ∗と𝑨について変分法を適用する。

δ𝐸s
𝜕𝐸s
∂𝐸s
=
−
∙
[𝜵
]=0
𝛿𝛹 ∗ ∂𝛹 ∗
𝜕𝜵𝛹 ∗

(1.4)

δ𝐸s 𝜕𝐸s∂𝐸s=− [𝜵 ∙
]=0
𝛿𝑨
∂𝑨
𝜕𝜵𝑨

(1.5)

(1.4)式と(1.5)式をそれぞれ解くと、以下の 2 式が得られる。

1(−iℏ𝜵 + 𝑒 ∗ 𝑨)2 𝛹 + 𝛼𝛹 + 𝛽|𝛹|2 𝛹 = 0
2𝑚∗

(1.6)

2

1
iℏ𝑒 ∗ ∗𝑒∗∗)(𝛹|𝛹|2 𝑨
𝜵×𝜵×𝑨=
𝜵𝛹
−
𝛹𝜵𝛹
−
𝜇0
2𝑚∗
𝑚∗

(1.7)

ここで、ゲージ変換には Coulomb ゲージ𝜵 × 𝑨 = 0を用いた。また、条件として超伝導体表面を横切って電流が流れないことを仮定した。この(1.6)、(1.7)式を G-L 方程式と呼ぶ。

1.2.2 TDGL 方程式
G-L 方程式は時間依存性を持たないため、時定数を導入して時間依存性を持った方程式に発展させる。これを Time-Dependent G-L(TDGL)方程式と呼ぶ。
2

𝛹と𝑨の時定数をそれぞれ𝛾, 𝜈とおくと、(1.4)式および(1.5)式は以下のように時間発展する形に書き換えられる。

δ𝐸s
𝜕𝐸s
∂𝐸s
𝜕𝛹
=
− [𝜵 ∙] = −𝛾
∗
∗
∗
𝛿𝛹
∂𝛹
𝜕𝜵𝛹
𝜕𝑡

(1.8)

δ𝐸s 𝜕𝐸s∂𝐸s
𝜕𝑨=− [𝜵 ∙] = −𝜈
𝛿𝑨
∂𝑨
𝜕𝜵𝑨
𝜕𝑡

(1.9)

これに以下のゲージ変換を与える。

𝜕𝛹 𝜕𝛹
→+ i𝑒 ∗ 𝑉𝛹𝜕𝑡𝜕𝑡(1.10)

𝜕𝑨 𝜕𝑨
→
+ 𝜵𝑉
𝜕𝑡
𝜕𝑡

(1.11)

ここで、V はスカラーポテンシャルである。(1.10)式および(1.11)式をそれぞれ(1.8)式、(1.9)
式に代入して解くと、

𝜕𝛹
1(−iℏ𝜵 + 𝑒 ∗ 𝑨)2 𝛹 + 𝛼𝛹 + 𝛽|𝛹|2 𝛹 = 0𝛾(+ i𝑒 ∗ 𝑉𝛹) +𝜕𝑡2𝑚∗(1.12)2𝜕𝑨
1iℏ𝑒 ∗ ∗𝑒∗∗)(𝛹|𝛹|2 𝑨 = 0
𝜈 ( + 𝜵𝑽) + 𝜵 × 𝜵 × 𝑨 +
𝜵𝛹
−
𝛹𝜵𝛹
−
𝜕𝑡
𝜇0
2𝑚∗
𝑚∗

(1.13)

が得られる。この(1.12)式、(1.13)式を解くことによって超伝導体の数値解析を行う。

1.2.3 TDGL 方程式の簡易化(1.12)式、(1.13)式をそのままコンピュータ上で解くことは困難であるため、細線近似と規格化という 2 つの簡易化を行う。
細線近似では、A が外部磁界 B にのみ依存すると仮定する。本研究の数値解析では、B は時間に対して一定に与えるから A も時間に対して一定となる。よって(1.13)式左辺第 1 項の時間偏微分が 0 となる。
次に(1.12)式、(1.13)式に対して規格化を行う。熱力学的臨界磁界を𝐻c として、超伝導体のコヒーレンス長𝜉と磁界侵入長𝜆を以下のように表す。

𝜉=ℏ√2𝑚∗ |𝛼|𝜆=𝑒 ∗ 𝜇0 𝐻c√𝑚∗ |𝛼|(1.14)(1.15)

また、以下の変換を定義する。

𝜉𝜵 → 𝜵

(1.16)|𝛼|𝑡→𝑡
𝛾(1.17)

3

𝑒 ∗𝛾
𝑉→𝑉
|𝛼|
𝜆
√2𝜇0 𝐻c(1.18)

𝑨→𝑨(1.19)

1

𝛽 2(1.20)( ) 𝛹→𝛹
|𝛼|(1.16)-(1.20)式の式変換を(1.14)式、(1.15)式に与えることによって規格化を行う。
まず(1.12)式の規格化を行う。(1.12)式の左辺第 1 項は、
1

1

|𝛼| 2
|𝛼|
|𝛼| 2
𝜕𝛹
𝜕
𝛾(
+ i𝑒 ∗ 𝑉𝛹) → 𝛾 [( ) 𝛹 + i𝑒 ∗ ( ∗ 𝑉) ( ) 𝛹]
𝛾
𝜕𝑡
𝑒 𝛾
𝛽
𝜕 ( 𝑡) 𝛽|𝛼|11|𝛼| 2 𝜕𝛹|𝛼| |𝛼| 2= 𝛾 [( )+i( ) 𝑉𝛹]𝛾𝛽𝛾 𝛽
𝜕 ( 𝑡)|𝛼|12 𝜕𝛹

(1.21)

1
2

|𝛼|
|𝛼|
1
1
= 𝛾 [ |𝛼| ( )+ |𝛼| ( ) i𝑉𝛹]
𝛾
𝛽
𝜕𝑡 𝛾
𝛽
1

|𝛼| 2 𝜕𝛹= |𝛼| ( ) [+ i𝑉𝛹]𝛽𝜕𝑡となる。(1.12)式の左辺第 2 項は、

1(−iℏ𝜵 + 𝑒 ∗ 𝑨)2 𝛹
2𝑚∗
1

2

|𝛼| 2
1
𝜵
√2𝜇0 𝐻c→𝑨) ( ) 𝛹(−iℏ − 𝑒 ∗
∗
2𝑚
𝜉
𝜆
𝛽
2

1

|𝛼| 2
= (−iℏ𝜵−𝑒𝑨) ( ) 𝛹𝛽√2𝑚∗ 𝜉
√2𝑚∗ 𝜆1𝐻∗ √2𝜇0 c(1.22)
21|𝛼| 2√2𝑚∗ |𝛼|√2𝜇0 𝐻c √𝑚∗ |𝛼|= ( ) (−iℏ𝜵 − 𝑒∗𝑨) 𝛹
𝛽√2𝑚∗ ℏ√2𝑚∗ 𝑒 ∗ 𝜇0 𝐻c
1|𝛼| 2= |𝛼| ( ) (−i𝜵 − 𝑨)2 𝛹𝛽となる。(1.12)式の左辺第 3 項は、
4

1

|𝛼| 2
𝛼𝛹 → 𝛼 ( ) 𝛹
𝛽(1.23)となる。(1.12)式の左辺第 4 項は、
1

2

1

|𝛼| 2
|𝛼| 2
𝛽|𝛹|2 𝛹 → 𝛽 |( ) 𝛹| ( ) 𝛹𝛽𝛽(1.24)= |𝛼| (12

|𝛼|
) |𝛹|2 𝛹
𝛽(1.21)-(1.24)式をまとめる。(1.23)式右辺の𝛼は負であるから、
1|𝛼| 2 𝜕𝛹𝛼( ) [+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 − 𝛹 + |𝛹|2 𝛹] = 0𝛽𝜕𝑡(1.25)となる。さらに、(1.25)式の両辺を𝛼(|𝛼|/𝛽)1/2で除算すると、

𝜕𝛹+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 − 𝛹 + |𝛹|2 𝛹 = 0
𝜕𝑡(1.26)となる。(1.26)式は超伝導領域のオーダパラメータ𝛹を記述する。続いて(1.13)式の規格化を行う。(1.13)式の左辺第 1 項は、

1 |𝛼|𝜈𝜵𝑉 → 𝜈 𝜵 ∗ 𝑉
𝜉 𝑒 𝛾
=

|𝛼|
𝜈𝜵𝑉
𝜉𝑒 ∗ 𝛾=√2𝑚∗ |𝛼| |𝛼|∙ ∗ 𝜈𝜵𝑉
ℏ𝑒 𝛾(1.27)となる。(1.13)式の左辺第 2 項は、

1
1 1
1√2𝜇0 𝐻c
𝜵×𝜵×𝑨→ ∙ 𝜵× 𝜵×
𝑨
𝜇0
𝜇0 𝜉
𝜉
𝜆
=

√2𝐻c
𝜵×𝜵×𝑨
𝜉2𝜆

2𝑚∗ |𝛼| √𝑚∗ |𝛼|=∙ ∗∙ √2𝐻c 𝜵 × 𝜵 × 𝑨ℏ2
𝑒 𝜇0 𝐻c=2√2𝑚∗ |𝛼|√𝑚∗ |𝛼|
𝜵×𝜵×𝑨ℏ2 𝑒 ∗ 𝜇0となる。(1.13)式の左辺第 3 項は、

5(1.28)

iℏe∗ ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )
2𝑚∗
1

1

1

1

|𝛼| 2
|𝛼| 2 1|𝛼| 2iℏe∗ |𝛼| 2 ∗ 1→{( ) 𝛹 𝜵 ( ) 𝛹 − ( ) 𝛹 𝜵 ( ) 𝛹 ∗ }
∗
2𝑚
𝛽
𝜉
𝛽
𝛽
𝜉
𝛽
=

iℏe∗ |𝛼| 1 ∗∙∙ (𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )2𝑚∗ 𝛽 𝜉=iℏe∗ |𝛼| √2𝑚∗ |𝛼| ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )
∙
∙
2𝑚∗ 𝛽
ℏ

=

i𝑒 ∗ |𝛼|√2𝑚∗ |𝛼| ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )2𝑚∗ 𝛽(1.29)となる。(1.13)式の左辺第 4 項は、
1

2

𝑒 ∗2
𝑒 ∗ 2 |𝛼| 2√2𝜇0 𝐻c
2
|𝛹|
𝑨
→
𝛹|
𝑨
|(
)
𝑚∗
𝑚∗
𝛽
𝜆
𝑒 ∗ 2 |𝛼| √2𝜇0 𝐻c|𝛹|2 𝑨
= ∗∙
∙
𝑚 𝛽
𝜆

(1.30)∗2=|𝛼|√𝑚∗ |𝛼|𝑒|𝛹|2 𝑨
∙
∙
𝐻
∙
√2𝜇
0
c
𝑚∗ 𝛽
𝑒 ∗ 𝜇0 𝐻c=𝑒 ∗ |𝛼|√2𝑚∗ |𝛼||𝛹|2 𝑨𝑚∗ 𝛽となる。(1.27)-(1.30)式をまとめると、

2√2𝑚∗ |𝛼|√𝑚∗ |𝛼|
𝜵×𝜵×𝑨ℏ2 𝑒 ∗ 𝜇0√2𝑚∗ |𝛼| |𝛼|𝑒 ∗ |𝛼|√2𝑚∗ |𝛼|
i
2
∗
∗ )}
(𝛹
=
𝑨
−
𝜵𝛹
−
𝛹𝜵𝛹
−
∙ ∗ 𝜈𝜵𝑉
{|𝛹|
𝑚∗ 𝛽
2
ℏ
𝑒 𝛾

(1.31)となる。(1.31)式の両辺を|𝛼|√2𝑚∗ |𝛼|で除算すると、

2𝑚∗
𝑒∗
i
1
𝜵
×
𝜵
×
𝑨
=
{|𝛹|2 𝑨 − (𝛹 ∗ 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )} − ∗ 𝜈𝜵𝑉2∗∗ℏ 𝑒 𝜇0𝑚 𝛽
2ℏ𝑒 𝛾(1.32)

が得られる。さらに、(1.32)式を簡単化するために以下の関係を用いる。

𝑱=

1
𝜵×𝜵×𝑨
𝜇0

(1.33)

𝜵∙𝑱=0(1.34)

ここで、J は超伝導体に流れる電流であり、(1.34)式は電流の発散を表す式である。(1.32)式の両辺で𝜵との内積を取ると、(1.33)式、(1.34)式から、

6

0=𝜵∙[↔

𝑒∗i1{|𝛹|2 𝑨 − (𝛹 ∗ 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )} − ∗ 𝜈𝜵𝑉]
∗
𝑚 𝛽
2
ℏ𝑒 𝛾

i
𝑚∗ 𝛽
(𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ ) − 𝜵 ∙ (|𝛹|2 𝑨) = −𝜈𝜵2 𝑉2ℏ𝑒 ∗ 2 𝛾(1.35)i↔ 𝜎𝜵2 𝑉 = (𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ ) − 𝜵 ∙ (|𝛹|2 𝑨)
2が得られる。ここで、−𝑚∗ 𝛽ℏ𝑒 ∗ 2 𝛾𝜈→𝜎(1.36)とする。(1.35)式はスカラーポテンシャル V を記述する。本研究の数値解析では、(1.29)式
および(1.35)式を解くことによって超伝導体の解析を行う。

1.3

アフィン変換数値積分法(Affine Integrator(AFI))

アフィン変換数値積分法(以下 AFI 法と記述する)とは、最近提案された陽的数値積分法
である。拡散方程式のような放物型偏微分方程式や、ゲージ場存在下の Schrödinger 方程式や Time-dependent Ginzburg-Landau(TDGL) 方程式の数値積分のために考案された。
数値的安定性に関して、AFI 法はラプラシアン項に関して無条件で安定であり、方程式を構成する項全体に関して高い数値的安定性を有する。エネルギー保存に関して、AFI 法
は線形保存系(Schrödinger 方程式や TDGL 方程式)において全エネルギーが厳密に保存される性質を持つ。
AFI 法は空間に関して離散化するための格子をチェッカーボード状に分解することから導かれるアフィン変換対で構成される。数学的にシンプルな手順であり、数値的安定性やエネルギー保存に関する理論解析が容易である。[1]1.3.1 AFI 法のイメージ

図 1-1 の左側に空間に関する離散化のためのチェッカーボード格子を示す。基礎方程式
はこのチェッカーボード格子によって、まずは空間に関して離散化される。右側にチェッ
カーボード格子に基づいて構成された高次元常微分方程式の概形を示す。
黒丸格子点の集合を Q 集合、白丸格子点の集合を P 集合と呼ぶことにする。Q 集合お
よび P 集合上で標本化された従属変数の時間発展は図 1-1 の右側に示された方程式のブロック行列が支配する。
図 1-2 に常微分方程式を時間に関して離散化するようすを示す。このように、AFI 法は
7

Q 集合上で定義された従属変数の時間発展と P 集合上で定義された従属変数の時間発展を交互に繰り返していくものである。

図 1-1 空間に関する離散化のイメージ

図 1-2 時間に関する離散化のイメージ

1.3.2 1 次元空間において AFI 法の数値積分スキームと時間発展方程式の導出

1 次元空間での AFI 法の数値積分スキームを導出する。
まずは空間に関する離散化を行う。本研究で用いるオーダパラメータの基礎方程式は(1.37)式で表される。
𝛾𝜕𝜓= (𝜵 − i𝑨)2 𝜓 − 𝛼𝜓 − 𝛽|𝜓|2 𝜓𝜕𝑡(1.37)

1 次元空間では(1.38)式が(1.37)式で適用されており、ベクトルポテンシャルは固定されているとする。
𝜵=𝜕, 𝑨 = 𝐴𝑥
𝜕𝑥

(1.38)

ベクトルポテンシャルはリンク変数𝑤𝑖,𝑗 によって次のように実装される。
𝑤𝑖𝑗 = exp(i𝜃𝑖𝑗 ) ,𝜃𝑖𝑗 = −ℎ𝐴𝑥 (8𝑥𝑖 + 𝑥𝑗),
2𝑖, 𝑗 = 1,2,3,4(1.39)

ここで𝑥𝑖 、𝑥𝑗 はそれぞれ格子点𝑖、𝑗の𝑥座標で、ℎは格子点間隔である。格子点𝑖と格子点𝑗が接続していない場合は𝑤𝑖,𝑗 = 0と約束する。簡単のため、周期的境界条件を採用する。空間に関する離散化のために導入した格子およびリンク変数の状況を図 1-3 に示す。

図 1-3 1 次元空間での周期的境界条件下の座標とリンク変数のインデックスの定義(1.37)式を空間に関して離散化する。𝜓𝑖 (𝑡) = 𝜓(𝑥𝑖 , 𝑡), 𝛼𝑖 = 𝛼(𝑥𝑖 ), 𝛽𝑖 = 𝛽(𝑥𝑖 )とおくと(1.40)式で表される連立微分方程式が得られる。
𝑑𝜓11= 2 ( 𝑤12 𝜓2 + 𝑤
̅41 𝜓4 − 2𝜓1 ) − 𝛼1 𝜓1 − 𝛽1 |𝜓1 |2 𝜓1
𝑑𝑡
ℎ
𝑑𝜓3
1
𝛾
= 2 ( 𝑤34 𝜓4 + 𝑤
̅23 𝜓2 − 2𝜓3 ) − 𝛼3 𝜓3 − 𝛽3 |𝜓3 |2 𝜓3
𝑑𝑡
ℎ
𝑑𝜓2
1
𝛾
= 2 ( 𝑤23 𝜓3 + 𝑤
̅12 𝜓1 − 2𝜓2 ) − 𝛼2 𝜓2 − 𝛽2 |𝜓2 |2 𝜓2
𝑑𝑡
ℎ
𝑑𝜓4
1
(
̅34 𝜓3 − 2𝜓4 ) − 𝛼4 𝜓4 − 𝛽4 |𝜓4 |2 𝜓4{ 𝛾 𝑑𝑡 = ℎ2 𝑤41 𝜓4 + 𝑤
𝛾(1.40)

ここで、𝑤
̅𝑖𝑗 は𝑤𝑖𝑗 の複素共役を表す。
奇数の格子点の組(図 1-3 の黒丸に対応)と偶数の格子点の組(図 1-3 の白丸に対応)の順に並べているが、図 1-1 の右側のイメージに従って表記している。この連立微分方程式は以下の(1.41)式で表せる。
−2 − 𝑈1 ℎ2
𝜓1
𝑑 𝜓
1
0
𝛾 ( 3) = 2
𝑑𝑡 𝜓2
ℎ
𝑤
̅12
𝜓4
𝑤
̅14
(

0
−2 − 𝑈3 ℎ2
𝑤
̅32
𝑤
̅34

𝑤12
𝑤32
−2 − 𝑈2ℎ2
0

𝑤14
𝜓1
𝑤34
𝜓3
( )
𝜓2
0
2
−2 − 𝑈4ℎ ) 𝜓4

(1.41)·
2
𝑈𝑖 = 𝛼𝑖 + 𝛽𝑖 |𝜓̃𝑖 | , 𝑖 = 1,2,3,4(1.42)

𝜓̃𝑖 は𝜓𝑖 の推定値である。理由は後述する。
ここで、
𝜓𝜓𝒒 = ( 1) , 𝒑 = ( 2)𝜓3𝜓4とする。次に、
9(1.43)

𝜎𝑖 = 2 + 𝑈𝑖 ℎ2, 𝑖 = 1,2,3,4(1.44)として、−𝜎𝐷(−𝜎𝑄 ) = ( 100−𝜎) , 𝐷(−𝜎𝑃 ) = ( 2
−𝜎3
0

𝑤12
0
) , 𝑊 = (𝑤
−𝜎4
32

𝑤14
𝑤34 )

(1.45)とすると、(1.41)式は(1.46)式で表せる。
𝛾d 𝒒
1 𝐷(−𝜎𝑄 )𝒒𝑊) (𝒑)(𝒑) = 2 (
𝑇
̅
dt
ℎ
𝑊
𝐷(−𝜎𝑃 )

(1.46)

時間推進演算子 𝐴̂を用いて(1.46)式を表すと、(1.47)式で表される。また、 𝐴̂は(1.48)で定
義される。
d 𝒒
𝒒
( ) = 𝐴̂ (𝒑)
dt 𝒑

𝛾

𝐴̂ =

1
h2

∑

(−𝜎𝑖 𝜓𝑖 +𝑖∈{1,2,3,4}∑(1.47)

𝑤𝑖𝑗 𝜓𝑗 )𝑗={1,2,3,4}𝜕𝜕𝜓𝑖

𝑄 = {1,3}, 𝑃 = {2,4}として時間推進演算子𝐴̂を次のように分解する。
𝐴̂ = 𝐴̂𝑄 + 𝐴̂P(1.48)(1.49)

ただし
𝐴̂𝑄 =

1
𝜕
∑ (−𝜎𝑖 𝜓𝑖 + ∑ 𝑤𝑖𝑗 𝜓𝑗 )h2
𝜕𝜓𝑖

(1.50)1𝜕∑ (−𝜎𝑖 𝜓𝑖 + ∑ 𝑤𝑖𝑗 𝜓𝑗 )2h𝜕𝜓𝑖(1.51)

𝑖∈𝑄

𝐴̂𝑃 =

𝑗∈𝑃

𝑖∈𝑃

𝑗∈𝑄

である。また、𝑤𝑖𝑗 = 𝑤
̅𝑗𝑖 であることに注意する。
次に、(1.47)式を時間に関して離散化することを考える。
時間刻み幅をτ、(𝒒, 𝒑) = (𝒒(t), 𝒑(t))、(𝒒′ , 𝒑′ ) = (𝒒(t + τ), 𝒑(t + τ))として時間発展方程式の(1.52)式が得られる。
𝜏
𝜏
𝒒
𝒒
𝒒′
( ′ ) = exp ( 𝐴̂) (𝒑) = exp ( (𝐴̂𝑄 + 𝐴̂P )) (𝒑)𝒑𝛾𝛾(1.52)

指数関数演算子は Lie-Trotter-Suzuki 分解によって(1.53)式のように近似できる。[4] LieTrotter-Suzuki 分解とは、exp 𝑥(𝐴 + 𝐵)の形のものを(1.54)式のように積に分解する方法の 1𝜏2つである。また、𝑂 ( 𝑛 )は誤差オーダー部分である。

10

𝜏
𝜏
exp ( 𝐴̂𝑄 ) exp ( 𝐴̂𝑃 ) + 𝑂(𝜏 2 )𝛾𝛾
𝜏exp ( (𝐴̂𝑄 + 𝐴̂P )) =𝜏𝜏𝜏𝛾
exp ( 𝐴̂𝑄 ) exp ( 𝐴̂𝑃 ) exp ( 𝐴̂𝑄 ) + 𝑂(𝜏 3 )
2𝛾
𝛾
2𝛾
…
{

(1.53)𝑥𝑥 𝑛
𝜏2
exp 𝑥(𝐴 + 𝐵) = (exp 𝐴 exp 𝐵) + 𝑂 ( )𝑛𝑛𝑛(1.54)

𝜏𝜏
従って、指数関数演算子exp (𝛾 𝐴̂𝑄 )とexp (𝛾 𝐴̂𝑃 )の行列表現が得られたら数値積分スキーム(AFI)を構成できる。
実験として、ここで演算子𝐴̂𝑄 を𝜓𝑖 (𝑖 ∈ 𝑄)に 1 回作用させてみる。
𝐴̂𝑄 𝜓𝑖 =1(−𝜎𝑖 𝜓𝑖 + ∑ 𝑤𝑖𝑗 𝜓𝑗 ) , 𝑖 ∈ 𝑄ℎ2(1.55)𝑗∈𝑃

次に、2 回作用させてみる。
−𝜎𝑖 21
𝐴̂2𝑄 𝜓𝑖 = ( 2 ) (𝜓𝑖 − ∑ 𝑤𝑖𝑗 𝜓𝑗 ) , 𝑖 ∈ 𝑄ℎ𝜎𝑖(1.56)𝑗∈𝑃

従って、𝑚回(𝑚は 1 以上の自然数)作用させると
−𝜎𝑖 𝑚
1
𝐴̂𝑚
𝜓
=
(
) (𝜓𝑖 − ∑ 𝑤𝑖𝑗 𝜓𝑗 ) , 𝑖 ∈ 𝑄
𝑄 𝑖
2
ℎ
𝜎𝑖

(1.56)𝑗∈𝑃となる。よって、マクローリン展開をして
∞

𝑚
𝜏
1 𝜏
exp ( 𝐴̂𝑄 ) 𝜓𝑖 = 𝜓𝑖 + ∑( 𝐴̂𝑄 ) 𝜓𝑖
𝛾
𝑚! 𝛾
𝑚=1
∞

= 𝜓𝑖 + ∑𝑚=1

1 −𝜎𝑖 𝜏 𝑚
1
( 2 ) (𝜓𝑖 − ∑ 𝑤𝑖𝑗 𝜓𝑗 )𝑚! ℎ 𝛾𝜎𝑖𝑗∈𝑃∞

1
1 −𝜎𝑖 𝜏 𝑚1
= 𝜓𝑖 + (𝜓𝑖 − ∑ 𝑤𝑖𝑗 𝜓𝑗 ) ∑() − (𝜓𝑖 − ∑ 𝑤𝑖𝑗 𝜓𝑗 )𝜎𝑖𝑚! ℎ2 𝛾
𝜎𝑖
𝑗∈𝑃

𝑚=1

𝑗∈𝑃

−𝜎 𝜏
1 − exp ( 2𝑖 )
−𝜎𝑖 𝜏
ℎ 𝛾
= exp ( 2 ) 𝜓𝑖 +∑ 𝑤𝑖𝑗 𝜓𝑗
ℎ 𝛾
𝜎𝑖
𝑗∈𝑃

(1.57)= 𝑎𝑖 𝜓𝑖 + 𝑏𝑖 ∑ 𝑤𝑖𝑗 𝜓𝑗 , 𝑖 ∈ 𝑄𝑗∈𝑃となる。ここで、11−𝜎𝑖 𝜏
1 − 𝑎𝑖𝑎𝑖 = exp ( 2 ) , 𝑏𝑖 =ℎ 𝛾
𝜎𝑖(1.58)

とおいた。(1.57)式は 1 次元空間における時間発展方程式である。また、𝐴̂𝑄 𝜓𝑖 = 0, 𝑖 ∈ 𝑃(1.59)であるので、
𝐴̂𝑚𝑄 𝜓𝑖 = 0, 𝑖 ∈ 𝑃(1.60)であることがわかる。演算子𝐴̂𝑃 でも同様な計算を行い、(1.61)式と(1.62)式が得られる。
𝜏𝒒𝒒exp ( 𝐴̂𝑄 ) (𝒑) = (𝐷(𝑎𝑄 ) 𝐷(𝑏𝑄 )𝑊) (𝒑)𝛾0𝐼(1.61)

𝜏𝒒0
exp ( 𝐴̂𝑃 ) (𝒑) = (
̅𝑇
𝐷(𝑏
𝛾
𝑃 )𝑊

(1.62)𝒒𝐼)( )𝐷(𝑎𝑃 ) 𝒑

ここで、𝐼は単位行列で𝑎
𝐷(𝑎𝑄 ) = ( 100𝑏) , 𝐷(𝑏𝑄 ) = ( 1
𝑎3
0

0
𝑎
) , 𝐷(𝑎𝑃 ) = ( 2
0
𝑏3

0
𝑏
) , 𝐷(𝑏𝑃 ) = ( 2
𝑎4
0

0
)
𝑏4

(1.63)とおいた。
𝜏𝜏

(1.61)式と(1.62)式はそれぞれ指数関数演算子exp (𝛾 𝐴̂𝑄 )とexp (𝛾 𝐴̂𝑃 )の行列表現を表しており、数値積分スキーム(AFI)が得られた。
2
次に、(1.42)式の非線形項𝛽𝑖 |𝜓̃𝑖 | で、𝜓𝑖 ではなく推定値 𝜓̃𝑖 が用いられる理由を説明する。推定値は(1.64)式から(1.67)式で計算される。
1𝜓̃1 = (𝑤12 𝜓2 + 𝑤
̅41 𝜓4 )2(1.64)1𝜓̃3 = (𝑤34 𝜓4 + 𝑤
̅23 𝜓2 )2(1.65)1𝜓̃2 = (𝑤23 𝜓3 + 𝑤
̅12 𝜓1 )2(1.66)1𝜓̃4 = (𝑤41 𝜓1 + 𝑤
̅34 𝜓3 )2(1.67)

注意することとして、𝜓𝑖 (𝑖 ∈ 𝑄)の値は𝜓𝑖 (𝑖 ∈ 𝑃)で推定され、𝜓𝑖 (𝑖 ∈ 𝑃)の値は𝜓𝑖 (𝑖 ∈ 𝑄)で推定されることである。このことにより、指数関数演算子の行列表現が容易になる。
𝜏

具体的に言うと、指数関数演算子exp (𝛾 𝐴̂𝑄 )の行列表現を得るために、(1.56)式において
𝐴̂𝑄 をベクトル𝒒の要素に何回も作用させた。その際、𝜎𝑖 (𝑖 ∈ 𝑄)はベクトル𝒒の要素に依存せず、ベクトル𝒑の要素に依存することが重要である。つまり、時間推進演算子𝐴̂𝑄 にとって𝜎𝑖 (𝑖 ∈ 𝑄)は定数となる。このことが行列表現の導出を可能としている。
12

1.3.3 2 次元空間において AFI 法の時間発展方程式の導出

2 次元空間での AFI 法の数値積分スキームを導出する。
基礎方程式は 1 次元空間で用いたものと同じで(1.37)式を用いる。ただし、
𝜕 𝜕
𝜵 = ( , ) , 𝑨 = (𝐴𝑥 , 𝐴𝑦 )𝜕𝑥 𝜕𝑦(1.68)である。
ここで、格子点やリンク変数を指定するためにダブルインデックスを用いる。そうする
ことで格子点やリンク変数を空間的な位置に対応させることができる。図 1-4 にダブルイ𝑦𝑥
ンデックスの定義を示す。また、リンク変数は𝑥成分の𝑤𝑖,𝑗と、𝑦成分の𝑤𝑖,𝑗に分けて定義さされる。境界条件は𝑥方向𝑦方向ともに周期的境界条件とする。

図 1-4 2 次元空間での周期的境界条件下の座標とリンク変数のダブルインデックスでの定義
まず空間に関して離散化すると(1.69)式が得られる。
γ

𝑑𝜓𝑖,𝑗
1
𝑦
𝑦
= 2 (𝑤𝑥𝑖,𝑗 𝜓𝑖+1,𝑗 + 𝑤
̅ 𝑥𝑖−1,𝑗 𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅ 𝑖,𝑗−1 𝜓𝑖,𝑗−1 − 4𝜓𝑖,𝑗 )𝑑𝑡ℎ2−𝛼𝑖,𝑗 𝜓𝑖,𝑗 − 𝛽𝑖,𝑗 |𝜓𝑖,𝑗 | 𝜓𝑖,𝑗

(1.69)

2
𝑈𝑖,𝑗 = 𝛼𝑖𝑗 + 𝛽𝑖,𝑗 |𝜓̃𝑖,𝑗 |(1.70)1 𝑥𝑦𝑦𝑥
̃ = (𝑤𝑖,𝑗𝜓
𝜓𝑖+1,𝑗 + 𝑤
̅𝑖−1,𝑗𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅𝑖,𝑗−1 𝜓𝑖,𝑗−1 )𝑖,𝑗4(1.71)

ここでは、13

とおく。また、
𝜎𝑖,𝑗 = 4 + 𝑈𝑖,𝑗ℎ2(1.72)

とおく。
ここまで格子点とリンク変数をダブルインデックスで定義したが、図 1-4 において格子点を(1,1)、(1,2)、…、(4,3)、(4,4)と定義しているところをシングルインデックス
1、２、…、15、16 と定義したものを考えると、1 次元空間で導出した数値積分スキームの(1.61)式と(1.62)式、時間発展方程式の(1.57)式がそのまま使える。数値積分スキームは
図 1-4 の黒丸に対して(1.61)式、白丸に対して(1.62)式が対応する。
具体例として格子点(1,1)での時間発展方程式を(1.73)式に示す。なお、𝜓1,1が時間発展した後のものを𝜓′1,1とする。(1.57)式を適用するが、使用されている格子点とリンク変数
に関してはダブルインデックスで定義した表記を使用している。ほかの格子点についても同様に(1.57)式を用いて導出可能である。
𝑦𝑦𝜓 ′1,1 = 𝑎1,1 𝜓1,1 + 𝑏1,1 (𝑤𝑥1,1 𝜓2,1 + 𝑤
̅ 𝑥4,1 𝜓4,1 + 𝑤1,1 𝜓1,2 + 𝑤
̅ 1,4 𝜓1,4 )(1.73)

𝑎1,1 = exp(−𝜎1,1 𝜏⁄ℎ2 )(1.74)(1.75)

𝑏1,1 = (1 − 𝑎1,1 )⁄𝜎1,1
2

𝜎1,1 = 4 + (𝛼1,1 + 𝛽1,1 |𝜓̃1,1 | ) ℎ2
𝜓̃1,1 =

(1.76)

1 𝑥𝑦𝑦𝑥(𝑤 𝜓 + 𝑤
̅4,1
𝜓4,1 + 𝑤1,1 𝜓1,2 + 𝑤
̅1,4 𝜓1,4 )
4 1,1 2,1

(1.77)

これを一般化すると、2 次元空間における AFI 法での時間発展方程式の(1.78)式を得ら
れる。
𝑦

𝑦

𝑥
𝑥
𝜓 ′ 𝑖,𝑗 = 𝑎𝑖,𝑗 𝜓𝑖,𝑗 + 𝑏𝑖,𝑗 (𝑤𝑖,𝑗𝜓𝑖+1,𝑗 + 𝑤
̅𝑖−1,𝑗𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅𝑖,𝑗−1 𝜓𝑖,𝑗−1 )(1.78)

𝑎𝑖,𝑗 = exp(−𝜎𝑖,𝑗 𝜏⁄ℎ2 )(1.79)

𝑏𝑖,𝑗 = (1 − 𝑎𝑖,𝑗 )⁄𝜎𝑖,𝑗(1.80)
2

𝜎𝑖,𝑗 = 4 + (𝛼𝑖,𝑗 + 𝛽𝑖,𝑗 |𝜓̃𝑖,𝑗 | ) ℎ2
𝜓̃𝑖,𝑗 =

1.4

1 𝑥
𝑦
𝑦
𝑥
(𝑤 𝜓
+𝑤
̅𝑖−1,𝑗𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅𝑖,𝑗−1 𝜓𝑖,𝑗−1 )4 𝑖,𝑗 𝑖+1,𝑗(1.81)(1.82)

研究目的

G-L 理論は第 2 種超伝導体の磁気特性を良く記述している。G-L 理論を成り立たせる仮定によって求められる G-L 方程式に時間発展性を加えた TDGL 方程式は第 2 種超伝導体内の磁束線の動きを調査するのに用いられる。TDGL 方程式に限らず、常微分方程式を数値解析する際によく用いられるのがオイラー法である。オイラー法は常微分方程式の数値
解法の 1 つで、微分の定義に基づいて数値解析するので数学的に簡単であり、プログラムを記述する上でも簡単である。しかし、数値解析する上では誤差が蓄積されるため精度が悪いことが知られている。
14

対して AFI 法は、オイラー法と同じく数学的に簡単であり、プログラムを記述する上でも簡単であることに加えて、方程式を構成する項全体に対して高い数値的安定性を有している。したがって、TDGL 方程式の数値解析を行う上で AFI 法を用いると高い安定性のもと数値解析を行うことができる。
本研究では、量子化磁束線可視化シミュレーションプログラムを作成する上で通常用いられるオイラー法で実装するところを AFI 法で実装し、超伝導体内の磁束線を可視化すること、および観察することを目的とする。

15

第2章 実装方法と計算条件
2.1

実装方法

実装するにあたって、統合開発環境でありプログラミング言語でもある Processing で実装を行う。Processing は、アーティストよるコンテンツ制作作業のために、詳細な設定を行
う関数を排除しており、非プログラマーでも容易にプログラムを実装できる。
シミュレーション領域は 2 次元平面であり、平面上に超伝導領域が存在するとしてプログ
ラムを実装する

2.2

初期条件

空間離散化格子を以下の図に示す。図 3-1 において、黒丸はオーダパラメータで、四角がリンク変数である。また、陰を付けている部分が超伝導領域で、その周囲が境界領域である。
なお、図 3-1 に振ってある番号は𝑥軸方向、𝑦軸方向それぞれのインデックス番号である。また、𝑥座標、𝑦座標はインデックス番号に格子間隔ℎをかけたものとなる。オーダパラメータの格
𝑦

𝑥
子点座標、
リンク変数の𝑥成分𝑤𝑖,𝑗と𝑦成分の𝑤𝑖,𝑗 の定義は図 1-4 のようにし、図 3-1 に適用する。

また、シミュレーション領域内においては黒丸と白丸の対応も適用する。

図 2-1 シミュレーション領域のオーダパラメータとリンク変数の配置の概念図
本稿では、シミュレーション領域は𝑥軸方向を𝑁𝑥 = 100で、𝑦軸方向を𝑁𝑥 = 100で分割する
16

ことにする。
オーダパラメータの初期条件は式(3.1)、ゲージ場の初期条件は式(3.2)とする。𝑥、𝑦はそ
れぞれ𝑥座標、𝑦座標を示す。
𝜓(𝑥, 𝑦, 0) = cos(𝜋𝑚𝑥 ⁄𝐿) cos(𝜋𝑛𝑦⁄𝐿)(3.1)(𝐴𝑥 , 𝐴𝑦 ) = (0, 𝐵𝑧 𝑥)(3.2)

ここで、𝐵𝑧 は 2 次元平面に垂直に印加された一様磁場である。また、式(3.1)の𝑚、𝑛は 1 以上の自然数であり、今回は𝑚 = 2、𝑛 = 3で実装することにする。

2.3

シミュレーション領域内での更新式

2.3.1 オーダパラメータの更新式
オーダパラメータの従属変数𝜓𝑖,𝑗 は式(3)に従って時間𝜏後の値𝜓′𝑖,𝑗 に更新される。
𝑦

𝑦

𝑥
𝑥
𝜓 ′ 𝑖,𝑗 = 𝑎𝑖,𝑗 𝜓𝑖,𝑗 + 𝑏𝑖,𝑗 (𝑤𝑖,𝑗𝜓𝑖+1,𝑗 + 𝑤
̅𝑖−1,𝑗𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅𝑖,𝑗−1 𝜓𝑖,𝑗−1 )

(3.3)

𝑥
𝑤𝑖,𝑗
= exp(− iℎ𝐴𝑥 (𝒓𝑖,𝑗 + 𝒓𝑖+1,𝑗 )⁄2)

(3.4)

𝑦
𝑤𝑖,𝑗 = exp(− iℎ𝐴𝑦 (𝒓𝑖,𝑗 + 𝒓𝑖,𝑗+1 )⁄2)

(3.5)

𝑎𝑖,𝑗 = exp(−𝜎𝑖,𝑗 𝜏⁄ℎ2 )

(3.6)

ただし

(3.7)

𝑏𝑖,𝑗 = (1 − 𝑎𝑖,𝑗 )⁄𝜎𝑖,𝑗
2

𝜎𝑖,𝑗 = 4 + (𝛼𝑖,𝑗 + 𝛽𝑖,𝑗 |𝜓̃𝑖,𝑗 | ) ℎ𝜓̃𝑖,𝑗 =

2

1 𝑥
𝑦
𝑦
𝑥
(𝑤 𝜓
+𝑤
̅𝑖−1,𝑗𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤
̅𝑖,𝑗−1 𝜓𝑖,𝑗−1 )4 𝑖,𝑗 𝑖+1,𝑗

(3.8)
(3.9)

である。シミュレーション領域のオーダパラメータの更新範囲は𝑖 = 1,2, … , 𝑁𝑥 、𝑗 = 1,2, … , 𝑁𝑦 である。ここで𝒓𝑖,𝑗 は格子点(𝑖, 𝑗)の位置ベクトル、ℎはシミュレーション領域内の格子間隔、𝜏は更新時間間隔である。式(3.4)と式(3.5)はそれぞれ𝑥軸方向、𝑦軸方向のリンク変数である。リンク
𝑦𝑥
変数の存在範囲は 、𝑤𝑖,𝑗が𝑖 = 0,1,2, … , 𝑁𝑥 、𝑗 = 0,1,2, … , 𝑁𝑦 + 1、𝑤𝑖,𝑗 が𝑖 = 0,1,2, … , 𝑁𝑥 + 1、𝑗 =

0,1,2, … , 𝑁𝑦 である。また、式(3.8)のパラメータ𝛼𝑖,𝑗 と𝛽𝑖,𝑗 はそれぞれ𝛼𝑖,𝑗 = −10、𝛽𝑖,𝑗 = 10で実装する。

2.3.2 リンク変数の更新式
ゲージ場の時間変化を表す式を、リンク変数を用いて(3.10)式と(3.11)式に示す。
𝜏𝐴

𝑥
𝑑𝑤𝑖,𝑗
1
𝑧
̅̅̅𝑧 𝑖,𝑗−1 𝐶𝑖,𝑗
= −iIm[𝜓̅𝑖,𝑗 𝑤𝑖𝑥 ,𝑗 𝜓𝑖+1,𝑗 ]𝑤𝑖𝑥 ,𝑗 − 2 (𝐶− 1)𝑤𝑖𝑥 ,𝑗
𝑑𝑡
ℎ
𝑦

𝜏𝐴

𝑑𝑤𝑖,𝑗
𝑑𝑡

𝑦

𝑦

= −iIm[𝜓̅𝑖,𝑗 𝑤𝑖 ,𝑗 𝜓𝑖+1,𝑗 ]𝑤𝑖 ,𝑗 −

ただし、
17

1
𝑦
𝑧
̅̅̅𝑧 𝑖,𝑗 𝐶𝑖−1,𝑗(𝐶− 1)𝑤𝑖 ,𝑗ℎ2(3.10)(3.11)

𝑦𝑧
𝐶𝑖,𝑗= 𝑤𝑖𝑥 ,𝑗 𝑤𝑖+1 ,𝑗 ̅̅̅̅𝑤 𝑥 𝑖,𝑗+1 ̅̅̅̅𝑤 𝑦 𝑖,𝑗(3,12)𝑦𝑥である。ここで、𝑥軸方向のリンク変数の偏角𝜃𝑖,𝑗と𝑦軸方向のリンク変数の偏角𝜃𝑖,𝑗 を用いて
𝑦𝑦𝑥
𝑤𝑖𝑥 ,𝑗 = exp(i𝜃𝑖,𝑗) , 𝑤𝑖+1 ,𝑗 = exp(i𝜃𝑖,𝑗 )
𝑥
𝜔𝑖,𝑗
=

𝑑 𝑥 𝑦
𝑑 𝑦
𝜃 ,𝜔 = 𝜃𝑑𝑡 𝑖,𝑗 𝑖,𝑗 𝑑𝑡 𝑖,𝑗(3.13)(3.14)

を用意する。すると次の関係式(3.15)式と(3.16)式が得られる。
𝑥
𝜔𝑖,𝑗
=−

𝑦

𝜔𝑖,𝑗 = −1
1 𝑧𝑧
Im [𝜓̅𝑖,𝑗 𝑤𝑖𝑥 ,𝑗 𝜓𝑖+1,𝑗 + 2 ̅̅̅
𝐶 𝑖,𝑗−1 𝐶𝑖,𝑗]𝜏𝐴ℎ

(3.15)1
1 𝑧
𝑦𝑧
Im [𝜓̅𝑖,𝑗 𝑤𝑖 ,𝑗 𝜓𝑖,𝑗+1 + 2 ̅̅̅
𝐶 𝑖,𝑗 𝐶𝑖−1,𝑗]𝜏𝐴ℎ(3.16)

これらの関係式より、(3.17)式と(3.18)式で表されるようなリンク変数の偏角に関するオイラー法が導き出せる。𝜏はオーダパラメータと同じ更新時間間隔である。
𝑥𝑥𝑥(𝑡 + 𝜏) = 𝜃𝑖,𝑗(𝑡) + 𝜏𝜔𝑖,𝑗
𝜃𝑖,𝑗
(𝑡)

(3.17)𝑦𝑥𝑥(𝑡) + 𝜏𝜔𝑖,𝑗
𝜃𝑖,𝑗 (𝑡 + 𝜏) = 𝜃𝑖,𝑗(𝑡)

(3.18)

リンク変数はリンク変数の偏角を更新し、その値を用いて更新時間後のリンク変数を求める形となる。また、プログラムで実装する際、安定性を考えて𝜏𝐴 = 10として実装する。ただし、𝑖 = 1,2, … , 𝑁𝑥 、𝑗 = 1,2, … , 𝑁𝑦 である。

2.4

境界条件

2.4.1 オーダパラメータの境界条件

オーダパラメータの境界条件は以下の式(3.19)から式(3.22)で示す。
𝜓0,𝑗 = 𝑤𝑥0 ,𝑗 𝜓1,𝑗

(3.19)

𝜓𝑁𝑥+1,𝑗 = ̅̅̅̅𝑤𝑥𝑁𝑥,𝑗 𝜓𝑁𝑥,𝑗(3.20)𝑦𝜓𝑖,0 = 𝑤𝑖 ,0 𝜓𝑖,1

(3.21)

𝜓𝑖,𝑁𝑦 +1 = ̅̅̅̅𝑤𝑦 𝑖,𝑁𝑦 𝜓𝑖,𝑁𝑦(3.22)

ただし、𝑖 = 1,2, … , 𝑁𝑥 、𝑗 = 1,2, … , 𝑁𝑦 である。

2.4.2 リンク変数の境界条件

超伝導体の𝑦軸方向に電流𝐽𝑦 を印加した場合の更新式を示す。
東西境界においてのリンク変数は以下の(3.23)式から(3.26)式に示す。
𝑥𝑤0𝑥 ,𝑗 = 𝑤1,𝑗18

(3.23)

𝑤𝑁𝑥𝑥 ,𝑗 = 𝑤𝑁𝑥𝑥 −1 ,𝑗𝑦𝑦𝑤0 ,𝑗 = 𝑤0𝑥 ,𝑗 𝑤1 ,𝑗 ̅̅̅̅𝑤 𝑥 0,𝑗+1 exp [iℎ2 (𝐵𝑧 + 𝜇0
𝑦(3.24)
𝐿𝑥
𝐽 )]
2 𝑦

𝑦

𝑤𝑁𝑥+1 ,𝑗 = ̅̅̅̅𝑤 𝑥 𝑁𝑥,𝑗 𝑤𝑁𝑥𝑥 ,𝑗+1 𝑤𝑁𝑥,𝑗 exp [−iℎ2 (𝐵𝑧 − 𝜇0

𝐿𝑥
𝐽 )]
2 𝑦

(3.25)(3.26)

ただし、(3.23)式と(3.24)式では𝑗 = 1,2, … , 𝑁𝑦 、(3.25)式と(3.26)式では𝑗 = 1,2, … , 𝑁𝑦 − 1である。また、(3.25)式と(3.26)式はアンペールの法則より容易に求まり、𝜇0 は真空透磁率、𝐿𝑥はシミュレーション領域の𝑥軸方向の大きさである。
南北境界におけるリンク変数の更新は以下の(3.27)式から(3.30)式に示す。
𝑥𝑥𝑤𝑖,0= 𝑤𝑖,1(3.27)𝑥𝑥𝑤𝑖,𝑁= 𝑤𝑖,𝑁𝑦 +1
𝑦𝑦𝑦𝑤𝑖,0 = 𝑤𝑖,1𝑦𝑦𝑤𝑖,𝑁𝑦 = 𝑤𝑖,𝑁𝑦−1(3.28)(3.29)(3.30)

ただし、(3.23)式と(3.24)式では𝑖 = 1,2, … , 𝑁𝑥 − 1、(3.25)式と(3.26)式では𝑖 = 1,2, … , 𝑁𝑥 である。

2.5

シミュレーション領域での更新順序

オーダパラメータは境界条件、図 1-4 における黒丸更新(0.5 更新時間分)、白丸更新(更新
時間分)、黒丸更新(0.5 更新時間分)の順で計算し、オーダパラメータの位相を色相で、オーダパラメータの絶対値の 2 乗の値を明度で表し、描画を行う。
リンク変数の更新は、境界条件を計算した後に超伝導領域内のリンク変数を計算する。
全体の更新順序は、リンク変数の境界条件処理、超伝導領域内のリンク変数の更新、オー
ダパラメータの境界条件処理、超伝導領域内のオーダパラメータの更新を 1 つのループとして更新を繰り返す。

19

第3章 実行結果および考察
3.1

実行結果

まず、電流を印加しなかった場合( 𝐽𝑦 = 0 )の結果を示す。第 2 章での計算条件を元にProcessing を用いてプログラムを実装した結果、以下の図 3-1 が得られた。ただし、オーダ
パラメータの位相を色相で、オーダパラメータの絶対値の 2 乗を明度で示している。また、Processing で描画する際、下が𝑦軸正の向き、右が𝑥軸正の向きになって描画されるので上が𝑦軸正の向き、右が𝑥軸正の向きになるように描画し直している。これより先の描画図は特に指定が無い限り図 3-1 と同様に描画しているものである。黒丸になっているところが量子化磁束線を示している。量子化磁束線が三角格子を組むように現れている。

−𝜋

0

𝜋

図 3-1 作成したプログラムの実行結果(𝐽𝑦 = 0の場合)
次に、y 軸方向への電流𝐽𝑦 = 1を印加したときの結果を示す。その結果図 3-2 が得られた。左右両方から磁束線が侵入し、内部の磁束線とぶつかると消滅するような動きが右半分の領域で観察された。

20

図 3-2 作成したプログラムの実行結果(𝐽𝑦 = 1の場合)

3.2

考察および今後の課題

𝐽𝑦 = 0でシミュレーションを実行して長時間放置すると図 3-3 のようになり、量子化磁束線が移動をやめてその場にとどまるようになった。これは、超伝導体内に電流が流れておらず、磁場のみが加えられているのでローレンツ力がはたらかず、磁束線がその場にとどまっているものだと考えられる。
次に、𝐽𝑦 = 1でシミュレーションを実行して左からだけでなく、右から磁束線が侵入することを考える。電流を𝑦軸正の向きに、磁場を𝑧軸正の向き、つまり奥から手前側に印加した
場合、磁束線はローレンツ力により𝑥軸正の向きへと一方通行で移動すると考えるのが自然である。これは、印加した磁場と電流の大きさの比率が関係しており、印加した磁場が電流と比較して大きい場合右側から磁束線が侵入する。[1]

21

−𝜋

0

𝜋

図 3-3 𝐽𝑦 = 0で長時間プログラムを実行した結果
今後の課題としては、3 次元に拡張した場合のプログラムの作成や、リンク変数の更新にオイラー法を用いて計算を行ったのでオーダパラメータだけでなくリンク変数の更新に
AFI 法を適用して計算するものを作成することが挙げられる。また、それを 3 次元に拡張して数値的安定性の高さから先行研究[3]で行われたことを精度良く調査することが期待される。

22

第4章 結言
本研究では、超伝導体内での量子化磁束線の可視化を行った。ただし、通常用いられるオイラー法ではなく AFI 法により TDGL 方程式を用いた計算を行い量子化磁束線の可視化を行った。
オイラー法で実装するときと同様に、AFI 法は数学的にもシンプルで、プログラムによる計算が容易であることが本研究で理解できた。また、今後の課題として、先行研究[3]では、
3 次元空間においてオイラー法でシミュレーションプログラムを作成していた。そして縦磁
界下での磁束線の動きについて調査を行っていた。このことから 3 次元空間でも AFI 法によってオーダパラメータとゲージ場の時間変化を計算し、精度良く磁束線の動きを調査することが今後の課題となる。

23

謝辞
日頃の研究指導に私の指導教員である情報工学研究院物理情報工学研究系の小田部荘司
教授に、超伝導の基礎知識に関しての指導に松下照夫名誉教授に、プログラム実装の指導に
有明工業高等専門学校創造工学科の松野哲也教授にお世話になりました。
小田部先生からは研究テーマ決めから、研究する上での助言をしていただきました。他にもバングラデシュやインドネシアの大学生との学術的な交流や文化的な交流の機会を下さったり、研究以外にも貴重な機会にたくさん招いていただき大変感謝しています。
松下先生からは、この研究の土台となる超伝導の知識を熱心に教えてくださいました。心より感謝申し上げます。
松野先生からは、2 月半ばまでも熱心にプログラムや AFI 法の指導をしていただきました。また、松野先生が参加された研究発表会に呼んでいただき、本稿完成の一助となりました。本当にありがとうございました。
また、小田部研究室メンバーの方々には公私ともに助言をいただいたり、直接研究に関係しない日常のコミュニケーションなどが、自分の心の大きな支えとなりました。この場を借りて深くお礼申し上げます。

本研究は JSPS 科研費 19H00771 の助成を受けたものです。

24

参考文献
[1] 松野哲也 有明工業高等専門学校創造工学科：私信による(2019 年)[2] Cases Reas, Ben Fry 著、船田巧

訳、Processing をはじめよう 第 2 版 (Make:

PROJECTS)： 2016 年 9 月 7 日
[3]増田

嘉道

九州工業大学大学院情報工学府

平成 27 年度修士学位論文

先端情報工学専攻 電子情報工学分野

時間依存 Ginzburg-Landau 方程式を用いた縦磁界下の超伝導

体内の磁束線に関する研究： 2016 年 2 月 9 日[4]中田 真秀 論文紹介:Lie-Trotter-鈴木分解について
https://www.slideshare.net/NakataMaho/lietrottersuzuki-159710114(参照 2020 年 2 月 15 日)

25