平成30年度
学士学位論文

時間依存Ginzburg-Landau方程式を用いた様々な条件下における超伝導体の量子化磁束線の運動解析

濱田 雄成(学生番号：15232074)

九州工業大学 情報工学部
電子情報工学科

指導教員：小田部 荘司 教授
平成31年2月18日

目次

第1章

序論 .................................................................................................................. - 7 -1.1

はじめに ................................................................................................................ - 7 -1.2

磁束ピンニング ..................................................................................................... - 8 -1.3(Time-Dependent) Ginzburg-Landau 方程式 ...................................................... - 11 -1.3.1

Ginzburg-Landau (G-L) 方程式 .................................................................. - 11 -1.3.2

Time-Dependent G-L (TDGL) 方程式 ........................................................ - 18 -1.3.3

TDGL 方程式の簡易化 .............................................................................. - 19 -1.4

近接効果 ............................................................................................................ - 24 -1.5

ピーク効果 .......................................................................................................... - 25 -1.6

超伝導体の異方性 ............................................................................................. - 26 -1.6.1

有効質量・有効導電率モデル ..................................................................... - 26 -1.6.2

Lawrence-Doniach モデル ........................................................................... - 27 -1.7

超伝導体の温度依存性...................................................................................... - 28 -1.8

本研究の目的 ..................................................................................................... - 29 -第2章

実装手法......................................................................................................... - 30 -2.1

式の離散化......................................................................................................... - 30 -2.1.1

微分の離散化.............................................................................................. - 30 -2.1.2

TDGL 方程式 (オーダーパラメータ) の離散化 ......................................... - 32 -2.1.3

TDGL 方程式(スカラーポテンシャル)の離散化 .......................................... - 40 -2.2𝐽cを得る方法 ...................................................................................................... - 46 -2.2.1

𝐸と常伝導電流 ............................................................................................ - 47 -1-

2.2.22.3

二分探索法の利用 ...................................................................................... - 48 -主要なアルゴリズム............................................................................................. - 52 -2.3.1

ガウス=ザイデル法 ...................................................................................... - 52 -2.3.2

オイラー法 ................................................................................................... - 53 -2.3.3

二分探索法 ................................................................................................. - 53 -2.4

GPGPU ............................................................................................................... - 55 -第3章

計算条件......................................................................................................... - 57 -3.1

SC ナノワイヤについての仮定 ............................................................................. - 57 -3.1.1

電流密度と外部磁束密度の条件 ................................................................ - 57 -3.1.2

オーダーパラメータとスカラーポテンシャルの初期条件 .............................. - 57 -3.1.3

境界条件 ..................................................................................................... - 57 -3.2

シミュレーションモデル ........................................................................................ - 59 -3.2.1

柱状ピンの形状 ........................................................................................... - 60 -3.2.2

基本条件 ..................................................................................................... - 60 -3.2.3

すべてのピン同士が平行な配置で磁束線に対して角度をとる場合 ............ - 61 -3.2.4

スプレーピンの場合 ..................................................................................... - 64 -3.2.5

星状ピン ...................................................................................................... - 67 -第4章

計算結果および考察....................................................................................... - 68 -4.1

基本条件 ............................................................................................................ - 68 -4.2

すべてのピン同士が平行な配置で磁束線に対して角度をとる場合 .................... - 69 -4.2.1

超伝導体に対するピンの角度を変更する場合 ............................................ - 69 -4.2.2

ピンと磁束線の鎖交体積からの考察 ........................................................... - 72 -4.2.3

球状ピンによる考察 ..................................................................................... - 74 -4.2.4

超伝導体に対する外部磁界の角度を変更する場合 ................................... - 77 -2-

4.3

スプレーピンの場合 ............................................................................................ - 81 -4.4

星状ピン ............................................................................................................. - 83 -第5章

結言 ................................................................................................................ - 86 -第6章

付録 ................................................................................................................ - 87 -6.1

GPGPU を用いたプログラミング .......................................................................... - 87 -6.1.1

for 文の削減 ................................................................................................ - 87 -6.1.2

時間計算量の削減 ...................................................................................... - 88 -6.1.3

条件分岐の削減.......................................................................................... - 89 -6.1.4

カーネルやカーネル引数の削減 .................................................................. - 89 -6.1.5

メモリコピー回数の削減 ............................................................................... - 89 -6.2

鎖交体積を求めるために用いた方法 .................................................................. - 90 -6.2.1

柱状ピンと磁束線が 1 本ずつの場合 .......................................................... - 90 -6.2.2

柱状ピンが 2 本ある場合............................................................................. - 93 -参考文献 ........................................................................................................................... - 94 謝辞 .................................................................................................................................. - 96 -

-3-

図目次
Fig. 1.1: 超伝導体に対して電流密度と磁束密度を垂直に与えることで発生する Lorentz
力𝑭L のようす......................................................................................................... - 9 Fig. 1.2: ピンを含む超伝導体内で𝑭L の反対方向に𝑭p が生じて磁束線の動きが妨げられているようす......................................................................................................... - 10 Fig. 1.3: 各種超伝導体におけるピーク効果 .............................................................. - 25 Fig. 2.1: シミュレーションによって得られた𝐸-𝐽特性(図中青線)と，その𝐸から磁束線運動に由来する成分だけを抽出した𝐸-𝐽特性(図中赤線)をあらわすグラフ .................... - 48 Fig. 2.2: 改変を加えた二分探索法のフローチャート .................................................. - 49 Fig. 2.3: 𝐽を増加させつつとった𝐸-𝐽特性(図中青線)と，𝐽を減少させつつとった𝐸-𝐽特性(図中赤線)をあらわすグラフ .................................................................................... - 50 Fig. 2.4: シミュレーション中で𝐽を増加させるとき(a)と減少させるとき(b)の挙動の違い - 51 Fig. 2.5: 二分探索法のフローチャートの一例 ............................................................ - 54 Fig. 3.1: 計算するシミュレーションモデルの立方体空間 ............................................ - 59 Fig. 3.2: 星状ピンの外観........................................................................................... - 67 Fig. 4.1: 基本条件における異方性の強さごとの𝐽c -𝐵特性 .......................................... - 68 Fig. 4.2: すべてのピン同士が平行な配置のまま磁束線に対して様々な角度をとった場合の𝐽c -𝜃p 特性 ........................................................................................................ - 69 Fig. 4.3: すべてのピン同士が平行な配置のまま磁束線に対して0°あるいは90°をとった場
合の𝐵 = 0.4における𝐽c -𝛾𝑧 特性 ........................................................................... - 70 Fig. 4.4: 横磁界に縦のピンが入っている場合の磁束線の変形 ................................. - 71 Fig. 4.5: ピンと磁束線の角度ごとの鎖交体積............................................................ - 72 Fig. 4.6: 磁束線がピンに沿うようす ............................................................................ - 73 -4-

Fig. 4.7: 球状ピンの外観........................................................................................... - 75 Fig. 4.8: 球状ピンによる𝐽c の𝛾𝑧 依存性 ....................................................................... - 76 Fig. 4.9: 基本条件のピン配置に対して，磁束線がさまざまな角度をとった場合の𝐽c -𝜃B 特
性........................................................................................................................ - 77 Fig. 4.10: 基本条件のピン配置に対して磁束線が20°あるいは50°をとった場合の𝐵 = 0.4
における𝐽c -𝛾𝑧 特性 ............................................................................................... - 78 Fig. 4.11: 磁束線によるピンの橋渡し(磁束線が𝑧 = 0の面から侵入し，𝑧 = 10の面から出ていくとき) ........................................................................................................... - 79 Fig. 4.12: 磁束線によるピンの橋渡し(磁束線が𝑦 = 0の面から侵入し，𝑦 = 10の面から出ていくとき) ........................................................................................................... - 80 Fig. 4.13: 楕円を描く磁束線の断面........................................................................... - 81 Fig. 4.14: スプレーピンの𝐽c -𝜃p 特性 ........................................................................... - 82 Fig. 4.15: 星状ピンの𝐽c -𝜃p 特性 ................................................................................. - 83 Fig. 4.16: 星状ピンと面状ピンの𝐽c -𝐵特性 ................................................................. - 84 Fig. 4.17: 星状ピンと比較するための面状ピン .......................................................... - 85 Fig. 6.1: 柱状ピンと磁束線が 1 本ずつの場合の鎖交体積を幾何学的に求めるための模式図 .................................................................................................................... - 90 Fig. 6.2: 𝐿𝑧 が有限である場合に鎖交体積を幾何学的に求めるための模式図 .......... - 91 -

-5-

表目次
Table 3.1: 基本条件におけるピンの中心線の座標 .................................................... - 60 Table 3.2: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合(ピンの角度を変更する場合)の細かな条件の一覧................................................................ - 61 Table 3.3: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合の，いくつかの𝜃p における柱状ピンの外観 ............................................................................. - 62 Table 3.4: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合(外部磁界の角度を変更する場合)の細かな条件の一覧 .................................................... - 63 Table 3.5: スプレーピンの場合の細かな条件の一覧 ................................................. - 64 Table 3.6: いくつかの𝜃p におけるスプレーピンの外観 ................................................. - 65 Table 4.1: 8 つの球状ピンの場合の細かな条件の一覧 .............................................. - 75 -

-6-

第1章
1.1

序論

はじめに

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

-7-

1.2

磁束ピンニング

超伝導体は，それが超伝導状態にあるときに，電流を流した場合や外部磁場をかけた場合，その磁気的な性質，振る舞いの違いにより，第一種超伝導体と第二種超伝導体に区別される．
第一種超伝導体は，電流および外部磁界を与えていない場合，その超伝導体の𝑇c 以下の温度において超伝導状態となり完全反磁性を示す．しかしこれに外部磁界を与えていくと，ある
外部磁界の大きさにおいてその超伝導状態が破壊されてしまう．この磁界は前述した𝐵c である．
一方で，第二種超伝導体は，第一種超伝導体と同じようにある磁界までは完全反磁性を示
すが，その磁界を超えると第一種超伝導体とは異なり，超伝導体内部に一定の磁束（磁束線）
を侵入させ，超伝導状態を維持することができる．磁束線を侵入させた領域は常伝導状態となるが，全体としては超伝導状態を維持している．この状態を混合状態と呼ぶ．これにさらに外部磁界を与えていくと，やがて超伝導状態が破壊される．第二種超伝導体の完全反磁性を示さ
なくなる転移磁界を𝐵c1 ，超伝導状態が破壊される転移磁界を𝐵c2 とする．
現在発見されている超伝導体では，第一種超伝導体の𝐵c と比較すると，第二種超伝導体の𝐵c2 は非常に大きいことが知られている．このために，工学的な応用には第二種超伝導体が用いられていることが一般的である．第二種超伝導体は，前述したとおり混合状態においては超伝導体内部に磁束線が侵入している（この磁束線の磁束密度を𝑩とする）．そのため，超伝導体に流す輸送電流（この電流密度を𝑱とする）により，その磁束線（正確には，その磁束線を留める渦糸電流）にLorentz力𝑭L が与えられる(Fig. 1.1)．

-8-

Fig. 1.1: 超伝導体に対して電流密度と磁束密度を垂直に与えることで発生するLorentz力𝑭L のようす

この𝑭L は𝑭L = 𝑱 × 𝑩(1.1)と表すことができる．また，この𝑭L により磁束線が速度𝒗で運動した場合，Josephsonの式より，誘導起電力
𝑬 =𝑩×𝒗(1.2)

が生じる．この𝑬は，𝑱と同じ向きに生じるので𝑱∙𝑬 > 0

(1.3)

となる．こうした状態が定常的に続くためには，この誘導起電力に見合った損失が発生しなければならない．すなわち，この𝑬は超伝導体に対してOhmicな損失をもたらすこととなり，超伝導体の超伝導状態を破壊する原因となる．しかしながら，実際の第二種超伝導体には磁束の運動を止める（𝒗 = 0）作用があり，第二種超伝導体に含まれる常伝導析出物，空隙，結晶粒界面-9-など，あらゆる欠陥や不均質物質がその作用をする．こうした欠陥などをピンニング・センターと呼び，それらの作用を磁束ピンニングと呼ぶ．磁束ピンニングは，𝑭L がある臨界値を超えるまで磁束線の動きを止めるため，𝑬による損失を生じさせないようにすることができる．単位体積当たりのピンニング・センターが磁束線に及ぼす力をピン力密度𝑭p とすると，超伝導体に𝑬が生じ始める電流密度（これを𝑱c とする）の下では，磁束線に単位体積当たりに𝐹L = 𝐽c 𝐵(1.4)のLorentz力が働いており，これが𝑭p と釣り合っていることから，
𝐽c =

𝐹p
𝐵

(1.5)

の関係がある．このときのようすをFig. 1.2に示す．(1.4)式の𝐽c を臨界電流密度という．第二種超伝導体は，𝑇c ，𝐵c2 ，𝐽c それぞれのパラメータが工学的な応用において重要となっている．

Fig. 1.2: ピンを含む超伝導体内で𝑭L の反対方向に𝑭p が生じて磁束線の動きが妨げられているようす- 10 -

(Time-Dependent) Ginzburg-Landau 方程式

1.3

1.3.1 Ginzburg-Landau (G-L) 方程式
Ginzburg-Landau(以降G-Lと記す)理論は，1950年にV. L. GinzburgとL. D. Landauによってロシアで提唱された，超伝導を説明する現象論である[1]．G-L理論は磁界と超伝導が共存する
場合の相転移を取り扱ったもので，とくに第二種超伝導体の磁気特性をよく記述することが知られている[2, 3]．
このG-L理論では，まず超伝導状態の秩序を表す量として，複素数のオーダーパラメータ
𝛹 = |𝛹|exp(𝑖𝜑)を定義し，𝑛s を超伝導電子密度として以下の式を満たすものとする．
|𝛹|2 ∝ 𝑛s(1.6)

超伝導状態の自由エネルギー𝐸s は𝑛s に依存しているため，(1.6)式より|𝛹|2 の関数となる．ここで，転移点近傍において|𝛹|2 は十分小さいと期待できるため，𝐸s は以下の式にように|𝛹|2 のべき展開ができる．
𝐸s = 𝐸n + 𝛼 |𝛹|2 +

𝛽
|𝛹 |4
2

(1.7)

𝐸nは常伝導状態の自由エネルギーである．また，超伝導-常伝導転移(以降S-N転移と記す)を記述するためには|𝛹|4 の項までの展開で良い．𝛼および𝛽はそれぞれべき展開した際の1次と2次の係数である．𝑇 < 𝑇c では𝛼 < 0，β > 0である．
次に，磁界の存在が𝛹の空間的変化に寄与することを考慮して，(1.7)式に磁界のエネルギ
ー密度と運動エネルギー密度を加算する．
𝐸s = 𝐸n + 𝛼 |𝛹|2 +𝛽1|𝛹|4 + (𝛁 × 𝑨)22
𝜇0

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

(1.8)

𝜇0 は真空中の透磁率，𝐴はベクトルポテンシャル，𝑚 ∗ は超伝導電子の質量，𝑒 ∗ は超伝導電子の電荷量，ℏはプランク定数，iは虚数単位である．

- 11 -

𝐸sを最小とするように，𝛹の共役複素数𝛹 ∗ と𝑨について変分法を適用する．
δ𝐸s
𝜕𝐸s
𝜕𝐸s
]=0
=
− [𝛁 ∙
∗
∗
𝛿𝛹
𝜕𝛹
𝜕𝛁𝛹 ∗

(1.9)

δ𝐸s 𝜕𝐸s𝜕𝐸s]=0=− [𝛁 ∙
𝛿𝑨
𝜕𝑨
𝜕𝛁𝑨(1.10)

まずは(1.9)式を解く．𝐸sを𝛹 ∗で偏微分すると次のようになる．
(𝛁 × 𝑨)2
𝜕𝐸𝑠
𝜕
𝛽
2
4
(
)
|
|
|
|
[𝐸
=
0
+
𝛼
𝛹
+
𝛹
+
𝜕𝛹 ∗ 𝜕𝛹 ∗ 𝑛
2
2𝜇0
+

1
{ℏ2 𝛁𝛹 ∙ 𝛁𝛹 ∗ + 𝑒 ∗2 |𝑨|2 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝑨
2𝑚 ∗∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}](1.11)
𝛽= 0 + 𝛼𝛹 + 𝛹𝛹 ∗ 𝛹
2+= 𝛼𝛹 +

1{0 + 𝑒 ∗2 |𝑨|2 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ (0 − 1 ∙ 𝛁𝛹)}
2𝑚 ∗

𝛽
1
(𝑒 ∗2 |𝑨|2 𝛹 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹)𝛹 |𝛹 |2 +22𝑚 ∗続いて𝐸sを𝛁𝛹 ∗ で偏微分したものの発散は次のようになる．
𝛁∙

𝜕𝐸𝑠
𝜕
𝛽
[𝐸𝑛 (0) + 𝛼|𝛹|2 + |𝛹|4
=𝛁∙
∗
∗
𝜕𝛁𝛹
𝜕𝛁𝛹
2
+

(𝛁 × 𝑨)22𝜇0

1{ℏ2 𝛁𝛹 ∙ 𝛁𝛹 ∗ + 𝑒 ∗2 |𝑨|2 𝛹𝛹 ∗+2𝑚 ∗+ iℏ𝑒 ∗ 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}]- 12 -(1.12)𝛁𝛹 ∗ を持たない項は微分され0となる．
𝛁∙

𝜕𝐸𝑠
𝜕
[0 + 0 + 0 + 0
=𝛁∙
∗
𝜕𝛁𝛹
𝜕𝛁𝛹 ∗
+
=

1
{ℏ2 𝛁𝛹 ∙ 𝛁𝛹 ∗ + 0 + iℏ𝑒 ∗ 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 0)}]
2𝑚 ∗

1
𝜕
(ℏ2 𝛁𝛹 ∙ 𝛁𝛹 ∗ + iℏ𝑒 ∗ 𝑨 ∙ 𝛹𝛁𝛹 ∗ )
𝛁∙
∗
2𝑚
𝜕𝛁𝛹 ∗

1
𝜕𝛁𝛹 ∗
𝜕𝛁𝛹 ∗
2
∗
=
𝛁
∙
(ℏ
𝛁𝛹
∙
+
iℏ𝑒
𝑨
∙
𝛹
)
2𝑚 ∗
𝜕𝛁𝛹 ∗
𝜕𝛁𝛹 ∗
=
=

(1.13)1𝛁 ∙ (ℏ2 𝛁𝛹 + iℏ𝑒 ∗ 𝑨 ∙ 𝛹)
2𝑚 ∗1{ℏ2 𝛁 ∙ 𝛁𝛹 + iℏ𝑒 ∗ 𝛁 ∙ (𝑨 ∙ 𝛹)}
2𝑚 ∗ここで𝛁 ∙ (𝑨 ∙ 𝛹)に対して積の微分公式を用いる．
𝛁∙

𝜕𝐸𝑠
1
[ℏ2 𝛁2 𝛹 + iℏ𝑒 ∗ {(𝛁 ∙ 𝐀)𝛹 + 𝑨 ∙ (𝛁𝛹)}]
=
∗
𝜕𝛁𝛹
2𝑚 ∗

(1.14)

クーロンゲージ𝛁 ∙ 𝑨 = 0より次のようになる．
𝛁∙

𝜕𝐸𝑠
1
(ℏ2 𝛁2 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹)=𝜕𝛁𝛹 ∗ 2𝑚 ∗(1.15)(1.11)式と(1.15)式から(1.9)式は次のようになる．
𝛼𝛹 +

𝛽
1
(𝑒 ∗2 |𝑨|2 𝛹 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹)𝛹 |𝛹 |2 +
2
2𝑚 ∗
−

𝛼𝛹 +

1
(ℏ2 𝛁2 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹) = 0
2𝑚 ∗

𝛽
1
(−ℏ2 𝛁2 𝛹 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹
𝛹 |𝛹 |2 +22𝑚 ∗+ 𝑒 ∗2 |𝑨|2 𝛹) = 0

𝛼𝛹 +

𝛽
1
(−ℏ2 𝛁2 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹 |𝛹 |2 +22𝑚 ∗+ 𝑒 ∗2 |𝑨|2 )𝛹 = 0- 13 -(1.16)

因数分解すると最終的に次のようになる．
𝛼𝛹 +

𝛽
1
(−iℏ𝛁 + 𝑒 ∗ 𝑨)2 𝛹 = 0
𝛹 |𝛹 |2 +22𝑚 ∗(1.17)

次に(1.10)式を解く．まずは𝐸s を𝑨で偏微分する部分を考える．𝑨はベクトルであるので，スカラ
ー𝐸𝑠 を𝑨で偏微分したものも次のように要素ごとに表現することができる．
𝜕𝐸𝑠 𝜕𝐸𝑠
𝜕𝐸𝑠
𝜕𝐸𝑠
=
𝒊𝒙 +
𝒊𝒚 +
𝒊
𝜕𝑨 𝜕𝐴𝑥𝜕𝐴𝑦
𝜕𝐴𝑧 𝒛(1.18)𝐴𝑥 , 𝐴𝑦 , 𝐴𝑧 はそれぞれ𝑨の𝑥, 𝑦, 𝑧成分である．ここで(1.18)式の𝑥成分に注目する．
2∗+ℏ 𝛁𝛹 ∙ 𝛁𝛹(𝛁 × 𝑨)2𝜕𝐸𝑠𝜕
𝛽1∗2 |𝑨|2=[𝐸𝑛 (0) + 𝛼|𝛹|2 + |𝛹|4 +
+
{
+𝑒
𝛹𝛹 ∗
}]
𝜕𝐴𝑥 𝜕𝐴𝑥
2
2𝜇0
2𝑚 ∗
∗
∗
(𝛹𝛁𝛹
+iℏ𝑒 𝑨 ∙− 𝛹 ∗ 𝛁𝛹)=+0
(𝛁 × 𝑨)2𝜕
1∗2 | |2+𝑒𝑨 𝛹𝛹 ∗[0 + 0 + 0 +
{
}]
+
𝜕𝐴𝑥
2𝜇0
2𝑚 ∗ +iℏ𝑒 ∗ 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)
=

1 𝜕
1 𝜕
+𝑒 ∗2 |𝑨|2 𝛹𝛹 ∗(𝛁 × 𝑨)2 +{}∗∗2𝜇0 𝜕𝐴𝑥
2𝑚 𝜕𝐴𝑥 +iℏ𝑒 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)2

𝜕𝐴𝑧 𝜕𝐴𝑦
+(
−
)
𝜕𝑦
𝜕𝑧
=

1 𝜕
2𝜇0 𝜕𝐴𝑥𝜕𝐴𝑥 𝜕𝐴𝑧 2
−
)
𝜕𝑧
𝜕𝑥
2
𝜕𝐴𝑦 𝜕𝐴𝑥
+(
−
)
𝜕𝑥
𝜕𝑦 }
{
+(

+𝑒 ∗2 (𝐴2𝑥 + 𝐴2𝑦 + 𝐴2𝑧 )𝛹𝛹 ∗1 𝜕+𝐴𝑥 𝒊𝒙+∗∗2𝑚 𝜕𝐴𝑥 +iℏ𝑒 (+𝐴𝑦 𝒊𝒚 ) ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)+𝐴𝑧 𝒊𝒛{}=1(0 + 0 + 0)2𝜇0+𝑒 ∗2 (2𝐴𝑥 + 0 + 0)𝛹𝛹 ∗1+1 ∙ 𝒊𝒙+∗ +0 ∙ 𝒊∗∗∗
2𝑚 +iℏ𝑒 (𝒚 ) ∙ (𝛹𝛁𝛹 − 𝛹 𝛁𝛹 )+0 ∙ 𝒊𝒛{}=1+2𝑒 ∗2 𝐴𝑥 𝛹𝛹 ∗{}∗∗2𝑚 +iℏ𝑒 𝒊𝒙 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)- 14 -(1.19)

𝑦成分，𝑧成分についても，(1.19)式と同様に変形する．
𝜕𝐸𝑠
1{2𝑒 ∗2 𝐴𝑦 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒚 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}=𝜕𝐴𝑦 2𝑚 ∗(1.20)

𝜕𝐸𝑠
1{2𝑒 ∗2 𝐴𝑧 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒛 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}=
𝜕𝐴𝑧 2𝑚 ∗(1.21)(1.19)式，(1.20)式，(1.21)式から得られる

𝜕𝐸𝑠,

𝜕𝐸𝑠 𝜕𝐸𝑠,𝜕𝐴𝑥 𝜕𝐴𝑦 𝜕𝐴𝑧

をすべて(1.18)式に代入すると次のようになる．
+{2𝑒 ∗2 𝐴𝑥 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒙 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}𝒊𝒙
𝜕𝐸𝑠
1[+{2𝑒 ∗2 𝐴𝑦 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒚 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}𝒊𝒚 ]=𝜕𝑨 2𝑚 ∗+{2𝑒 ∗2 𝐴𝑧 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒛 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}𝒊𝒛(1.22)

𝜕𝐸𝑠
1{2𝑒 ∗2 𝑨𝛹𝛹 ∗ + iℏ𝑒 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}
=+
𝜕𝑨
2𝑚 ∗

(1.23)

成分をまとめる．

続いて𝐸sを𝛁𝑨で偏微分したものの発散も要素ごとに表現できる．
𝛁∙

𝜕𝐸𝑠
𝜕𝐸𝑠
𝜕𝐸𝑠
𝜕𝐸𝑠
=𝛁∙
𝒊𝒙 + 𝛁 ∙𝒊𝒚 + 𝛁 ∙
𝒊
𝜕𝛁𝑨
𝜕𝛁𝐴𝑥
𝜕𝛁𝐴𝑦
𝜕𝛁𝐴𝑧 𝒛(1.24)

ここで(1.24)式の𝑥成分に注目する．
𝛁∙

𝜕𝐸𝑠
𝜕 𝜕𝐸𝑠
𝜕 𝜕𝐸𝑠
𝜕 𝜕𝐸𝑠
=
∙
+
∙
+
∙
∂𝐴
∂𝐴
𝜕𝛁𝐴𝑥 𝜕𝑥 𝜕 𝑥 𝜕𝑦 𝜕 𝑥 𝜕𝑧 𝜕 ∂𝐴𝑥
∂𝑥
∂𝑦
∂𝑧

(1.25)

𝐸𝑠 のうち微分された𝐴𝑥 を含むのは，ローテーションを持つ項すなわち(1.8)式でいうところの右辺第4項だけであるので，それ以外の項は0となり，(1.25)式は次のようになる．

- 15 -

𝛁∙

𝜕𝐸𝑠
𝜕𝛁𝐴𝑥

𝜕𝐴𝑦 2
𝜕𝐴
+ ( 𝜕𝑦𝑧 − 𝜕𝑧 )𝜕𝐴𝑥 𝜕𝐴𝑧 2−)𝜕𝑧𝜕𝑥𝜕𝐴𝑦 𝜕𝐴𝑥 2+(𝜕𝜕 {𝜕𝑥 − 𝜕𝑦 ) }=∙𝜕𝑥 𝜕 ∂𝐴𝑥
2𝜇0
∂𝑥
+(

2
𝜕𝐴𝑧 𝜕𝐴𝑦+ ( 𝜕𝑦 − 𝜕𝑧 )𝜕𝐴𝜕𝐴 2
+ ( 𝜕𝑧𝑥 − 𝜕𝑥𝑧 )𝜕𝐴𝑦 𝜕𝐴𝑥 2
+
(
−
)
𝜕
𝜕 {
𝜕𝑥
𝜕𝑦 }
+
∙
𝜕𝑦 𝜕 ∂𝐴𝑥
2𝜇0
∂𝑦

𝜕𝐴𝑦 2
𝜕𝐴
+ ( 𝜕𝑦𝑧 − 𝜕𝑧 )𝜕𝐴𝑥 𝜕𝐴𝑧 2−)𝜕𝑧𝜕𝑥𝜕𝐴𝑦 𝜕𝐴𝑥 2+(𝜕𝜕 {𝜕𝑥 − 𝜕𝑦 ) }+ ∙𝜕𝑧 𝜕 ∂𝐴𝑥
2𝜇0∂𝑧+(∂𝐴𝑥 ∂𝐴𝑥 ∂𝐴𝑥∂𝑥, ∂𝑦 , ∂𝑧 による偏微分を行う．

- 16 -(1.26)

𝛁∙

𝜕𝐸𝑠
𝜕𝛁𝐴𝑥

+0
{+0}
𝜕 +0
=
𝜕𝑥 2𝜇0+0+0{}𝜕𝐴𝑦 𝜕𝐴𝑥
(
)
+2
(
−
)
−1
𝜕
𝜕𝑥
𝜕𝑦
+
𝜕𝑦
2𝜇0
+0
𝜕𝐴𝑥 𝜕𝐴𝑧( )}{+2 (𝜕𝑧 − 𝜕𝑥 ) +1
𝜕
+0
+
𝜕𝑧
2𝜇0

(1.27)𝜕 𝜕𝐴𝑦 𝜕𝐴𝑥
𝜕 𝜕𝐴𝑥 𝜕𝐴𝑧
𝜕𝑦 ( 𝜕𝑥 − 𝜕𝑦 ) − 𝜕𝑧 ( 𝜕𝑧 − 𝜕𝑥 )
=−
𝜇0
𝜕
𝜕
(rot𝑨|𝑧 ) − (rot𝑨|𝑦 )
𝜕𝑦
𝜕𝑧
=−
𝜇0
=−

rot(rot𝑨)|𝑥𝜇0

𝑦成分，𝑧成分についても，(1.27)式と同様に変形する．
𝛁∙rot(rot𝑨)|𝑦
𝜕𝐸𝑠
=−
𝜕𝛁𝐴𝑦
𝜇0

(1.28)

𝛁∙𝜕𝐸𝑠rot(rot𝑨)|𝑧
=−
𝜕𝛁𝐴𝑧
𝜇0

(1.29)(1.27)式，(1.28)式，(1.29)式から得られる𝛁 ∙

𝜕𝐸𝑠
𝜕𝛁𝐴𝑥

,𝛁 ∙

𝜕𝐸𝑠
𝜕𝛁𝐴𝑦

,𝛁 ∙

𝜕𝐸𝑠
𝜕𝛁𝐴𝑧

をすべて(1.24)式に代入す

ると次のようになる．
𝛁∙rot (rot𝑨)|𝑦𝜕𝐸𝑠
rot (rot𝑨)|𝑥rot (rot𝑨)|𝑧
=−
𝒊𝒙 −
𝒊𝒚 −
𝒊𝒛
𝜕𝛁𝑨
𝜇0
𝜇0
𝜇0
rot(rot𝑨)=−
𝜇0- 17 -(1.30)

(1.11)式，(1.17)式，(1.23)式，(1.30)式を(1.9)式と(1.10)式に代入すると，次の2式になる．
1(−iℏ𝛁 + 𝑒 ∗ 𝑨)2 𝛹 + α𝛹 + 𝛽|𝛹|2 𝛹 = 02𝑚 ∗(1.31)

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

(1.32)

これら(1.31)式と(1.32)式を，G-L方程式と呼ぶ．

1.3.2 Time-Dependent G-L (TDGL) 方程式
Time-Dependent G-L(以降TDGLと記す)方程式は，G-L方程式に時間依存性を付与したものである． (1.9)及び(1.10)式に対して時間発展する場合を考えると，以下の2式となる．
δ𝐸s
𝜕𝐸s
𝜕𝐸s
𝜕𝛹
] = −γ=− [𝛁 ∙
∗
∗
∗
𝛿𝛹
𝜕𝛹
𝜕∇𝛹
𝜕𝑡

(1.33)

δ𝐸s 𝜕𝐸s
𝜕𝐸s
𝜕𝑨] = −ν
=− [𝛁 ∙
𝛿𝑨
𝜕𝑨
𝜕𝛁𝑨
𝜕𝑡

(1.34)

γとνはそれぞれ𝛹と𝑨の時定数である．さらにゲージ変換を与えるとℏ∂𝛹𝜕→ (ℏ + i𝑒 ∗ 𝑉) 𝛹𝜕𝑡𝜕𝑡

(1.35)

∂𝑨 𝜕𝑨
→
+ 𝛁𝑉
𝜕𝑡
𝜕𝑡

(1.36)となる．ここで，𝑉はスカラーポテンシャルである．(1.35)，(1.36)式を(1.33)，(1.34)式それぞれに代入すると最終的に次のようになる．
𝛾(
𝜈(

𝜕𝛹
1
(−iℏ𝛁 + 𝑒 ∗ 𝑨)2 𝛹 + α𝛹 + 𝛽|𝛹|2 𝛹 = 0+ i𝑒 ∗ 𝑉𝛹) +
𝜕𝑡
2𝑚 ∗

𝜕𝑨
1
iℏ𝑒 ∗ ∗𝑒 ∗𝟐∗)(
|𝛹 |𝟐 𝑨
+ 𝛁𝑉) + 𝛁 × 𝛁 × 𝑨 +
𝛹
𝛁𝛹
−
𝛹𝛁𝛹
−
𝜕𝑡
𝜇0
2𝑚 ∗
𝑚∗
=0

- 18 -(1.37)(1.38)

1.3.3 TDGL 方程式の簡易化
本研究ではTDGL方程式を数値解析するが，(1.37)，(1.38)式をそのまま解くのは困難である
ため，細線近似と規格化による2つの簡易化を行う[4, 5]．
細線近似では，非常に細い超伝導線 (以降SCナノワイヤと記す) に外部磁界𝑩ext を印加した時，SCナノワイヤ全体に𝑩ext が侵入すると仮定する．これにより，𝑨は𝑩ext にのみ依存する変数となる．本研究で𝑩ext は一定とするため，(1.38)式は，左辺第一項の時間偏微分が0となり一定となる．
次に(1.37)，(1.38)に対して規格化を行う．超伝導体のコヒーレンス長𝜉と磁界侵入長𝜆は熱力学的臨界磁界を𝐻c として以下の2式のように表す．
𝜉=ℏ√2𝑚 ∗ |𝛼 |𝜆=𝑒 ∗ 𝜇0 𝐻c√𝑚 ∗ |𝛼|(1.39)(1.40)

そして，(1.41)~(1.45)式に記す規格化を行う．
𝜉𝜵 → 𝜵

(1.41)|𝛼 |𝑡→𝑡
𝛾(1.42)

𝑒 ∗𝛾
𝑉→𝑉
|𝛼 |

(1.43)

𝜆𝑨→𝑨(1.44)𝛽 2( ) 𝛹→𝛹
|𝛼 |(1.45)√2𝜇0 𝐻c1

すると，(1.37)式左辺第1項は次のようになる．

- 19 -

𝛾(𝜕𝛹+ i𝑒 ∗ 𝑉𝛹)𝜕𝑡1|𝛼 | 2
𝜕
→ γ[
( ) 𝛹
𝛾
𝜕 (| | 𝑡) 𝛽
𝛼
1

|𝛼 |
|𝛼 | 2+ i𝑒 ∗ ( ∗ 𝑉) ( ) 𝛹]
𝑒 𝛾
𝛽
1

1

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

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

|𝛼 | 2 𝜕𝛹= |𝛼 | ( ) [+ i𝑉𝛹]
𝛽𝜕𝑡(1.37)式左辺第2項は次のようになる．
1(−𝑖ℏ𝜵 + 𝑒 ∗ 𝑨)2 𝛹
2𝑚 ∗
1

2

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

1

|𝛼 | 2= (−iℏ
𝜵
−
𝑒
𝑨)
(
) 𝛹
𝛽
√2𝑚 ∗ 𝜉√2𝑚 ∗ 𝜆1𝐻∗ √2𝜇0 c(1.47)
21|𝛼 | 2√2𝑚 ∗ |𝛼 |𝐻 √𝑚 ∗ |𝛼|∗ √2𝜇0 c= ( ) (−iℏ
𝜵
−
𝑒
𝑨) 𝛹
𝛽
√2𝑚 ∗ ℏ
√2𝑚 ∗ 𝑒 ∗ 𝜇0 𝐻c1|𝛼 | 2= |𝛼 | ( ) (−i𝜵 − 𝑨)2 𝛹
𝛽(1.37)式左辺第3項は次のようになる．

- 20 -

1|𝛼 | 2α𝛹 → 𝛼 ( ) 𝛹
𝛽(1.48)(1.37)式左辺第4項は次のようになる．
1
2

2

1

|𝛼 |
|𝛼 | 2
𝛽 |𝛹|2 𝛹 → 𝛽 |( ) 𝛹| ( ) 𝛹𝛽𝛽(1.49)12|𝛼 |= |𝛼 | ( ) |𝛹 |2 𝛹𝛽(1.46)~(1.49)式をまとめると，(1.48)式右辺の𝛼は負で，次のようになる．
1|𝛼 | 2 𝜕𝛹𝛼( ) [+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 − 𝛹 + |𝛹|2 𝛹] = 0𝛽𝜕𝑡(1.50)

1|𝛼| 2(1.50)式の両辺を𝛼 ( 𝛽 ) で割ると，最終的に次にようになる．

𝜕𝛹+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 − 𝛹 + |𝛹|2 𝛹 = 0
𝜕𝑡(1.51)(1.38)式も同様に簡易化する．(1.38)式左辺第一項は，𝑨が一定として，次のようになる．
1 |𝛼 |
𝜈𝜵𝑉 → 𝜈 𝜵 ∗ 𝑉
𝜉 𝑒 𝛾
=
=

|𝛼 |
𝜉𝑒 ∗ 𝛾𝜈𝜵𝑉(1.52)√2𝑚 ∗ |𝛼| |𝛼 |∙ ∗ 𝜈𝜵𝑉
ℏ𝑒 𝛾(1.38)式左辺第2項は次のようになる．
11 1
1√2𝜇0 𝐻c
𝜵×𝜵×𝑨→ ∙ 𝜵× 𝜵×
𝑨
𝜇0
𝜇0 𝜉
𝜉
𝜆
=

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

2𝑚 ∗ |𝛼| √𝑚 ∗ |𝛼|=∙ ∗∙ √2𝐻c 𝜵 × 𝜵 × 𝑨ℏ2
𝑒 𝜇0 𝐻c=2√2𝑚 ∗ |𝛼 |√𝑚∗ |𝛼|
𝜵×𝜵×𝑨ℏ2 𝑒 ∗ 𝜇0- 21 -(1.53)

(1.38)式左辺第3項は次のようになる．
iℏe∗ ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )
2𝑚 ∗
1

1

1

1

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

iℏe∗ |𝛼 | 1 ∗∙∙ (𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )2𝑚 ∗ 𝛽 𝜉(1.54)iℏe∗ |𝛼 | √2𝑚 ∗ |𝛼| ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )=∙∙2𝑚 ∗ 𝛽ℏi𝑒 ∗ |𝛼 |√2𝑚 ∗ |𝛼| ∗(𝛹 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )=2𝑚 ∗ 𝛽
(1.38)式左辺第4項は次のようになる．
1

2

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

(1.55)

∗2

=

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

(1.52)~(1.55)式をまとめると，
2√2𝑚 ∗ |𝛼|√𝑚∗ |𝛼|
𝜵×𝜵×𝑨ℏ2 𝑒 ∗ 𝜇0𝑒 ∗ |𝛼 |√2𝑚 ∗ |𝛼 |i√2𝑚 ∗ |𝛼 |
2
∗
∗ )}
(
{|𝛹|
=
𝑨
−
𝛹
𝜵𝛹
−
𝛹𝜵𝛹
−
𝑚∗ 𝛽
2
ℏ
∙

(1.56)

|𝛼 |
𝜈𝜵𝑉
𝑒∗𝛾

となり，(1.56)式の両辺を|𝛼 |√2𝑚 ∗ |𝛼 |で割ると，次のようになる．
2𝑚 ∗
𝜵×𝜵×𝑨
ℏ2 𝑒 ∗ 𝜇0𝑒∗i1= ∗ {|𝛹|2 𝑨 − (𝛹 ∗ 𝜵𝛹 − 𝛹𝜵𝛹 ∗ )} − ∗ 𝜈𝜵𝑉
𝑚 𝛽
2
ℏ𝑒 𝛾
- 22 -(1.57)

ここで，本研究ではTDGL方程式を𝛹と𝑉について解くが，現段階では変数が2つに対して方程式が(1.51)式のみであるため，電流の発散の式
𝛁∙𝑱=0(1.58)

1
𝛁×𝛁×𝑨
𝜇0(1.59)

を第2の方程式として解く．ここで𝑱=である．(1.57)式の両辺に𝜵と内積を取ると
𝑒∗i10 = 𝛁 ∙ [ ∗ {|𝛹|2 𝑨 − (𝛹 ∗ 𝛁𝛹 − 𝛹𝛁𝛹 ∗ )} − ∗ 𝜈𝛁𝑉]
𝑚 𝛽
2
ℏ𝑒 𝛾
𝑒∗ i
1
↔ ∗ { (𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ ) − 𝜵 ∙ (|𝛹|2 𝑨)} = − ∗ 𝜈𝜵2 𝑉
𝑚 𝛽 2
ℏ𝑒 𝛾
i
𝑚∗ 𝛽
2
↔ (𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ ) − 𝜵 ∙ (|𝛹|2 𝑨) = −2 𝜈𝜵 𝑉∗2ℏ𝑒 𝛾(1.60)i↔ 𝜎𝛁2 𝑉 = (𝛹 ∗ 𝛁2 𝛹 − 𝛹𝛁2 𝛹 ∗ ) − 𝛁 ∙ (|𝛹|2 𝑨)
2が得られる．ここで−𝑚∗ 𝛽ℏ𝑒 ∗ 2 𝛾𝜈→𝜎(1.61)とする．𝜎は常伝導領域の導電率である．本研究では(1.51)式と(1.60)式について，数値解析を行う．
一 方 ， SC ナ ノ ワイ ヤ の 対 破 壊 電 流 密 度 𝐽d を考 え る と ， 1 次 元 の G-L 方 程 式 は 𝛹(𝑥) =𝑓exp (i𝑞𝑥)に対して𝑨 = 0として，(1.51)式より
0 = 𝛁2 𝛹 + (1 − |𝛹|2 )𝛹
𝑑2= 2 (𝑓𝑒 i𝑞𝑥 ) + (1 − 𝑓 2 )𝑓𝑒 i𝑞𝑥
𝑑𝑥(1.62)= (−𝑞2 + 1 − 𝑓 2 )𝑓𝑒 i𝑞𝑥
ゆえに
𝑓 = 1 − 𝑞2が成り立つ．このパラメータ𝑓と𝑞は𝐽s から決定される．
- 23 -(1.63)

𝐽s = |𝛹|2

𝜕𝜑𝜕𝑥(1.64)= 𝑓2 𝑞= 𝑞(1 − 𝑞2 )𝑞(1 − 𝑞2 )は𝑞 =

1

2

√3

≈ 0.577のとき最大値3√3 ≈ 0.385を取るため，次の不等式が言える．
𝐽s ≤ 𝐽d =23√3≈ 0.385(1.65)(1.65)式より，規格化された電流密度𝐽をおよそ0.385を超えて印加すると，対破壊電流を超えた電流を流していることになるため，注意する必要がある．

1.4

近接効果

本研究では超伝導領域の内部にピンニング・センターとして様々な形状を持った常伝導析
出物(以降ピンと記す)を導入するが，この超伝導体とピンが接触した系において，ピンの領域に超伝導電子が緩やかに染み出すことが知られている[6]．これによってピンが超伝導性を示す
現象を近接効果 (proximity effect) という．ピンの領域では，(1.51)式に代わる方程式として
𝜕𝛹+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 + 𝜂𝛹 = 0
𝜕𝑡(1.66)

𝜂 = √𝜉/𝜉n(1.67)

を解く[7-10]．ここで，であり，𝜉n はピンの領域のコヒーレンス長である．(1.60)及び(1.61)式より，𝜉n が大きいほど𝜂が小さくなり，示す超伝導性が強いピンとなる．一方，非現実的であるが，ピンの領域に超伝導電子の染み出しが全くない理想的なピンを考慮した場合，ピンの領域において(1.40)式に代わる方程式を，次のように書くことができる．
𝛹=0
今回はこちらの式を用いて計算を行う．

- 24 -(1.68)

1.5

ピーク効果

超伝導体における臨界電流密度𝐽c は通常の場合，外部磁界の増加とともに単調に減少するはずであるが，ある磁界の大きさにおいて𝐽c のピークを示すことがある．この現象をピーク効果(peak effect) という[2]．このピーク効果が観測される超伝導体は，金属系超伝導体や銅酸化物超伝導体がある．ここで， (1.39)と(1.40)式の比で表されるGLパラメータ𝜅 = 𝜆/𝜉

(1.69)

が十分小さい超伝導体は，上部臨界磁界𝐵c2 近傍において𝐽c のピークが現れることが経験的に知られている．そして，𝜅がそこから大きくなるにつれて𝐽c のピークは相対的に低磁界側に移動する．特に𝜅が非常に大きい銅酸化物超伝導体では，かなり低い磁界で𝐽c のピークが現れる．これを特にフィッシュ・テール効果という(Fig. 1.3参照)[2, 11, 12, 13]．

Fig. 1.3: 各種超伝導体におけるピーク効果(a)𝜅が低いNb-50 at%Ta，(b)𝜅が高いTi-22 at%Nb及び(c)Y-Ba-Cu-o高温超伝導体𝐵 ∥ 𝑐軸．
- 25 -

この𝐽c のピークをもたらす要因として，様々な機構が提案されているが，本研究ではその中の1つであるマッチング機構に注目する．

1.6

超伝導体の異方性

多くの超伝導体では𝑎𝑏面内のコヒーレンス長は𝑐軸方向とは異なり，つまり異方的である．このことは，銅酸化物超伝導体において最も顕著である[27]．
超伝導体の異方性を再現するために，次に挙げるような2つのモデルがある．これら2つはアプローチと実装の点で異なるが，シミュレーション時のふるまいとしては同じようになることがわかっている．今回は有効質量・有効導電率モデルを用いるものとする．𝑥, 𝑦, 𝑧軸方向への異方性の強さ(倍率)をそれぞれ𝛾𝑥 , 𝛾𝑦 , 𝛾𝑧 とし，これらを各軸への異方性パラメータと呼ぶ．等方的な場
合，異方性パラメータはすべて1である．

1.6.1 有効質量・有効導電率モデル
異方的な超伝導体内のふるまいを観測するとき，電子質量および導電率が見かけ上，軸ごとに異なっていると捉えることができる．このときの，見かけの質量と導電率をそれぞれ，有効質量と有効導電率と呼ぶ．これらを次のように与えることで異方性を再現することができる[35]．
𝛾𝑥2𝑚⋆ → ( 0
0

0
𝛾𝑦2
0

0
0 ) 𝑚⋆
𝛾𝑧2

1
𝛾𝑥2

0

0

0

1
𝛾𝑦2

0

0

0

1
𝛾𝑧2 )

𝜎→
(

𝜎

(1.70)(1.71)

まずオーダーパラメータに関する式についてだが，(1.51)式の左辺第3項に(1.71)式が反映されることになるので，ここに注目し，異方性パラメータを含ませる．
- 26 -

2
𝜕
− 𝐴𝑥 )
𝜕𝑥
+
𝛾𝑥2
2
𝜕
(−i
−
𝐴
)
𝑦
𝜕𝑦
(−i𝛁 − 𝑨)2 𝛹 →
𝛹
+
2
𝛾𝑦

(−i

(−i
{

+

(1.72)

2𝜕− 𝐴𝑧 )
𝜕𝑧}
𝛾𝑧2

次にスカラーポテンシャルに関する式を考える．(1.60)式が次のようになる．
1
𝛾𝑥2

0

0

0

1
𝛾𝑦2

0

0

0

1
𝛾𝑧2 )

𝜎
(

=

i
2
(

1
𝛾𝑥2

0

0

0

1
𝛾𝑦2

0

0

0

1
𝛾𝑧2 )

1
𝛾𝑥2

0

0

0

1
𝛾𝑦2

0

0

0

1
𝛾𝑧2 )

−𝜵 ∙ |𝛹|2(

𝜵2 𝑉(𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ )(1.73)

𝑨

あとは2.1節で述べるように離散化すれば実装できる形となる．

1.6.2 Lawrence-Doniach モデル
Lawrence-Doniach(以降L-Dと記す)モデルは，1970年にW. E. LawrenceとS. Doniachによって示された[25]．このモデルは，超伝導薄膜が層状に積層し，その層間が弱いジョセフソン結合で影響し合うようすを再現する．すなわち3次元TDGL方程式というよりも，2次元TDGL方程式の積み重ねを行うやりかたである．
超伝導体の𝑎𝑏面がシミュレーションの𝑥-𝑦平面に一致している場合を例にとると，(1.51)式の- 27 -

左辺第3項の𝑧成分を置き換えることで実装できる[36]．具体的にはその項を次のようにする．
2𝜕− 𝐴𝑥 )
𝜕𝑥2
𝜕(−i𝛁 − 𝑨)2 𝛹 = + (−i− 𝐴𝑦 ) 𝛹𝜕𝑦2𝜕
{ + (−i 𝜕𝑧 − 𝐴𝑧 ) }

+ (−i

2
𝜕
+ (−i
− 𝐴𝑥 )
𝜕𝑥
2
𝜕
→
+ (−i
− 𝐴𝑦 )𝜕𝑦−i𝐴ℎ𝑧𝛹 (𝑥, 𝑦, 𝑧 + ℎ )e+ 𝛹(𝑥, 𝑦, 𝑧 − ℎ )ei𝐴𝑧ℎ − 2𝛹 (𝑥, 𝑦, 𝑧)−{}
𝛾𝑧2 ℎ 2(1.74)

ここでℎは超伝導層どうしの層間距離である．ただし空間を離散化する今回のシミュレーション
では，層間距離はシミュレーションにおける空間離散幅と等しいものとする．あとは2.1節で述べるように離散化すれば実装できる形となる．

1.7

超伝導体の温度依存性

超伝導体のふるまいは，その温度に大きく左右される．超伝導体がその性質を保つためには，
𝑇c を下回る温度でなければならないのは，はじめに述べたとおりである．しかし，𝑇c を超えない範囲の温度変化であっても，超伝導体はそのふるまいを変化させる．それを再現するため，(1.37)
式の𝛼に温度依存性を持たせて規格化を行うと次のようになる．
𝜕𝛹+ i𝑉𝛹 + (−i𝜵 − 𝑨)2 𝛹 − (1 − 𝑇)𝛹 + |𝛹|2 𝛹 = 0
𝜕𝑡(1.75)

ここで𝑇は0以上1未満の値をとる温度パラメータである．𝑇 = 0で絶対零度，𝑇 = 1で臨界温度と等しくなる．なお特に断りがない限り，𝑇 = 0の状態でを行う．

- 28 -

1.8

本研究の目的

過去のさまざまな研究から，超伝導体内のピンの条件の違いなどによって，その超伝導体の臨界電流密度𝐽c が変化することが知られている[18–24, 28, 29]．超伝導材料の工業的な応用において，より大きな𝐽c が得られる超伝導材料の開発は極めて重要な課題である．そのため，より大きな𝐽c が得られる超伝導体を調査する研究が求められている．
これまで，TDGLシミュレーションを用いて球状ピン・柱状ピン・面状ピンといった形状のピンについて論じられてきた．その中で柱状ピンは磁束線に沿うような形をしていることから，必要最低限の常伝導体積で十分なピン力を発揮することができるとされている．また，ピーク効果を引
き出しやすいことから，条件によっては面状ピンをも凌ぐピン力を持ちうる．
しかしながら柱状ピンはその形状ゆえ，磁束線に平行である場合でなければ，本来の性能を出せないということも知られている．実際の運用の場面で，材料に対して磁界の向きが完全に一方向でなければならないというのは，大きな制約である．
また，これまでのシミュレーションでは等方的な超伝導体のみを取り扱ってきた．実用される
多くの超伝導体は異方的であるから，等方的な超伝導体のシミュレーションは実用的とは言えなかった．
そこで，本研究では横磁界下における真空中のSCナノワイヤにおいて，柱状ピンの磁束線とのなす角の違いが𝐽c にどのように影響するのかについて，また，超伝導体の異方性が𝐽c にどのように影響するのかについて，簡易化したTDGL方程式を数値解析する．さらには外部磁界の向きに大きく影響されないような柱状ピンを提案することを目的とする．

- 29 -

第2章

実装手法

TDGL方程式をはじめとしたここまでの理論が用意されていたとしても，それを使えば即座に数値解が得られたり，磁束線の動きを見ることができたりするわけではない．それらを実現するには，理論式を数値的に計算するという過程が必要である．しかしTDGL方程式の計算には，次に挙げるような複数の障壁がある．
まず，時間および空間に関して連続的であるTDGL方程式は，連続的なまま数値計算することができない．SCナノワイヤ内のあらゆる地点およびあらゆる時点に対して(1.51)式と(1.60)式を適用することができるが，実数パラメータの集合としてのすべての点を列挙することはできない
ため，不可能である．そこで，TDGL方程式を，時間および空間について離散化する必要があ
る．
次に，それらの計算量は莫大であることが挙げられる．TDGL方程式を離散化し，計算量が無限大でなくなったとはいえ，超伝導体のふるまいを知るだけの精度を得るには，非常に多くの計算の繰り返しが必要となる．したがって，種々の計算を正確かつ高速に行えるような工夫が求
められる．

2.1

式の離散化
2.1.1 微分の離散化(1.51)式と(1.60)式には微分演算が含まれる．微分や積分は式の連続性を前提とした演算であるから，これを離散的に扱えるように変換しなければならない．1 変数関数𝑓(𝑥)の𝑎付近の 1
階微分は，次のような極限である．
𝜕𝑓𝑓(𝑎 + 𝑘) − 𝑓(𝑎)(𝑎) = lim
𝑘→0
𝜕𝑥
𝑘

(2.1)

これを離散化したいという要求から，𝑘を 0 に近づけるのではなく，非零の値𝑘𝑁 とするとよい．そうすると，次のようになる．
- 30 -

𝜕𝑓𝑓 (𝑎 + 𝑘𝑁 ) − 𝑓(𝑎)(𝑎 ) =
𝜕𝑥
𝑘𝑁

(2.2)

また，𝑘𝑁 が十分に小さいのであれば，次のようにも書ける．
𝜕𝑓𝑓 (𝑎) − 𝑓(𝑎 − 𝑘𝑁 )(𝑎 ) =
𝜕𝑥
𝑘𝑁

(2.3)

𝜕𝑓
𝑓 (𝑎 + 𝑘𝑁 ) − 𝑓(𝑎 − 𝑘𝑁 )
(𝑎) =
𝜕𝑥
2𝑘𝑁

(2.4)

次に，2 階微分について考える．2 階微分は，1 階微分を 2 度使って次のように書くことができ
る．
𝜕𝑓
𝜕
𝜕 2𝑓
(𝑎) = 𝜕𝑥 (𝑎)
𝜕𝑥 2
𝜕𝑥

(2.5)

𝜕𝑓
𝜕𝑓
(𝑎 + 𝑘𝑁 ) −
( )
𝜕2𝑓
𝜕𝑥
𝜕𝑥 𝑎
(
)
𝑎
=
𝜕𝑥 2
𝑘𝑁

(2.6)

これに(2.2)式を代入する．

これに(2.3)式を代入する．
𝜕2𝑓
(𝑎) =
𝜕𝑥 2

𝑓(𝑎 + 𝑘𝑁 ) − 𝑓 (𝑎 − 𝑘𝑁 + 𝑘𝑁 ) 𝑓(𝑎 ) − 𝑓(𝑎 − 𝑘𝑁 )
−
𝑘𝑁
𝑘𝑁
𝑘𝑁

𝑓 (𝑎 + 𝑘𝑁 ) − 𝑓(𝑎) − 𝑓(𝑎) + 𝑓(𝑎 − 𝑘𝑁 )
𝑘𝑁
=
𝑘𝑁
=

𝑓(𝑎 + 𝑘𝑁 ) − 2𝑓 (𝑎 ) + 𝑓 (𝑎 − 𝑘𝑁 )𝑘𝑁2

- 31 -(2.7)

2.1.2 TDGL 方程式 (オーダーパラメータ) の離散化
まずは(1.51)式を離散化する．(1.51)式の第1項に(2.2)式を代入し，第3項を展開すると次のようになる．
+

𝛹(𝑡 + 𝛥𝑡) − 𝛹(𝑡)𝛥𝑡+i𝑉(𝑡)𝛹 (𝑡)−𝜵2 𝛹 (𝑡)+2i𝑨 ∙ 𝜵𝛹 (𝑡)(2.8)+|𝑨|2 𝛹(𝑡)−𝛹(𝑡)+|𝛹(𝑡)|2 𝛹(𝑡)=0ここで，𝛥𝑡 は時間の離散幅である．𝛹は複素数であるので，このままでは数値計算できない．そこ
で，式全体を実部と虚部に分解する．(2.9)式，(2.10)式が，それぞれ(2.8)式の実部と虚部に対応する．
+

𝛹R (𝑡 + 𝛥𝑡) − 𝛹R (𝑡)𝛥𝑡−𝑉(𝑡)𝛹I (𝑡)−𝜵2 𝛹R (𝑡)−2𝑨 ∙ 𝜵𝛹I (𝑡)+|𝑨|2 𝛹R (𝑡)−𝛹R (𝑡)+{𝛹R (𝑡)2 + 𝛹I (𝑡)2 }𝛹R (𝑡)=0- 32 -(2.9)

𝛹I (𝑡 + 𝛥𝑡) − 𝛹I (𝑡)𝛥𝑡+𝑉 (𝑡)𝛹R (𝑡)−𝜵2 𝛹I (𝑡)+2𝑨 ∙ 𝜵𝛹R (𝑡)(2.10)+|𝑨|2 𝛹I (𝑡)−𝛹I (𝑡)+{𝛹R (𝑡)2 + 𝛹I (𝑡)2 }𝛹I (𝑡)=0ただし，𝛹R , 𝛹I はそれぞれ𝛹 の実部と虚部であり，次の式を満たすものとする．
𝛹R + 𝛹I = 𝛹(2.11)(2.9)式を𝛹R (𝑡 + 𝛥𝑡)について整理し，(2.10)式を𝛹I (𝑡 + 𝛥𝑡)について整理すると，次の 2 式が得られる．
+𝑉(𝑡)𝛹I (𝑡)+𝛁2 𝛹R (𝑡)+2𝑨 ∙ 𝛁𝛹I (𝑡)𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥𝑡−|𝑨|2 𝛹R (𝑡)+𝛹R (𝑡)2[−{𝛹R (𝑡) + 𝛹I (𝑡)2 }𝛹R (𝑡)](2.12)−𝑉(𝑡)𝛹R (𝑡)+𝛁2 𝛹I (𝑡)−2𝑨 ∙ 𝛁𝛹R (𝑡)𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + 𝛥𝑡−|𝑨|2 𝛹I (𝑡)+𝛹I (𝑡)2[−{𝛹R (𝑡) + 𝛹I (𝑡)2 }𝛹I (𝑡)](2.13)- 33 -

今回はこれを，2 次のシンプレクティック法を参考にしたオイラー法で解く．そのためにさらに次のように分割する．
+𝑉(𝑡)𝛹I (𝑡)+𝛁2 𝛹R (𝑡)11+2𝑨 ∙ 𝛁𝛹I (𝑡)𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥𝑡−|𝑨|2 𝛹R (𝑡)22+𝛹R (𝑡)2[−{𝛹R (𝑡) + 𝛹I (𝑡)2 }𝛹R (𝑡)](2.14)

1−𝑉(𝑡)𝛹R (𝑡 + 𝛥𝑡)2+𝛁2 𝛹I (𝑡)1−2𝑨 ∙ 𝛁𝛹R (𝑡 + 𝛥𝑡)
2
𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + 𝛥𝑡−|𝑨|2 𝛹I (𝑡)+𝛹I (𝑡)2
1− {𝛹R (𝑡 + 𝛥𝑡) + 𝛹I (𝑡)2 } 𝛹I (𝑡)2[](2.15)1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡 + 𝛥𝑡)
2+𝑉(𝑡)𝛹I (𝑡 + 𝛥𝑡)
1+𝛁2 𝛹R (𝑡 + 𝛥𝑡)
2+2𝑨 ∙ 𝛁𝛹I (𝑡 + 𝛥𝑡)
1
1−|𝑨|2 𝛹R (𝑡 + 𝛥𝑡)
+ 𝛥𝑡
2
2
1
+𝛹R (𝑡 + 𝛥𝑡)22
11− {+𝛹R (𝑡 + 2 𝛥𝑡) } 𝛹R (𝑡 + 𝛥𝑡)
2[]+𝛹I (𝑡 + 𝛥𝑡)2(2.16)

これら 3 式を順に解くと，次の時点の𝛹 が判明するというしくみである．これらを，次は空間について離散化する．まずはナブラ演算を展開する．

- 34 -

+𝑉(𝑡)𝛹I (𝑡)𝜕2+ 2 𝛹R (𝑡)𝜕𝑥𝜕2+ 2 𝛹R (𝑡)𝜕𝑦𝜕2+ 2 𝛹R (𝑡)𝜕𝑧11𝜕𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥𝑡
+2𝑨 ∙ 𝛹I (𝑡)
2
2
𝜕𝑥
𝜕
+2𝑨 ∙
𝛹 (𝑡)
𝜕𝑦 I
𝜕
+2𝑨 ∙ 𝛹I (𝑡)𝜕𝑧−|𝑨|2 𝛹R (𝑡)+𝛹R (𝑡)2[−{𝛹R (𝑡) + 𝛹I (𝑡)2 }𝛹R (𝑡)](2.17)

1−𝑉(𝑡)𝛹R (𝑡 + 𝛥𝑡)2
𝜕2+ 2 𝛹I (𝑡)𝜕𝑥𝜕2+ 2 𝛹I (𝑡)𝜕𝑦𝜕2+ 2 𝛹I (𝑡)
𝜕𝑧
𝜕
1
−2𝑨 ∙
𝛹R (𝑡 + 𝛥𝑡)
𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + 𝛥𝑡
𝜕𝑥
2
𝜕
1
−2𝑨 ∙
𝛹R (𝑡 + 𝛥𝑡)𝜕𝑦2𝜕
1−2𝑨 ∙ 𝛹R (𝑡 + 𝛥𝑡)
𝜕𝑧
2
2 ( )
|
|
− 𝑨 𝛹I 𝑡+𝛹I (𝑡)2
1− {𝛹R (𝑡 + 𝛥𝑡) + 𝛹I (𝑡)2 } 𝛹I (𝑡)2[](2.18)- 35 -

𝛹R (𝑡 + 𝛥𝑡)
1= 𝛹R (𝑡 + 𝛥𝑡)
2+𝑉(𝑡)𝛹I (𝑡 + 𝛥𝑡)
𝜕21+ 2 𝛹R (𝑡 + 𝛥𝑡)
𝜕𝑥
2
𝜕2
1
+ 2 𝛹R (𝑡 + 𝛥𝑡)
𝜕𝑦
2
2
𝜕
1
+ 2 𝛹R (𝑡 + 𝛥𝑡)
𝜕𝑧
2
𝜕
+2𝑨 ∙
𝛹 (𝑡 + 𝛥𝑡)
𝜕𝑥 I
1
𝜕
+ 𝛥𝑡
+2𝑨
∙
𝛹 (𝑡 + 𝛥𝑡)2
𝜕𝑦 I𝜕
+2𝑨 ∙ 𝛹I (𝑡 + 𝛥𝑡)𝜕𝑧
1
2−|𝑨| 𝛹R (𝑡 + 𝛥𝑡)2
1+𝛹R (𝑡 + 𝛥𝑡)22
11− {+𝛹R (𝑡 + 2 𝛥𝑡) } 𝛹R (𝑡 + 𝛥𝑡)
2[]+𝛹I (𝑡 + 𝛥𝑡)2
これら 3 つの式の該当箇所に，(2.4)式および(2.7)式を適用する．

- 36 -(2.19)

+𝑉(𝑡)𝛹I (𝑡)+𝛹R (𝑥 + ℎ, 𝑦, 𝑧, 𝑡)
1+ 2 { −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡) }ℎ+𝛹R (𝑥 − ℎ, 𝑦, 𝑧, 𝑡)+𝛹R (𝑥, 𝑦 + ℎ, 𝑧, 𝑡)
1+ 2 { −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡) }ℎ+𝛹R (𝑥, 𝑦 − ℎ, 𝑧, 𝑡)+𝛹R (𝑥, 𝑦, 𝑧 + ℎ, 𝑡)
1+ 2 { −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡) }11ℎ𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥+𝛹R (𝑥, 𝑦, 𝑧 − ℎ, 𝑡)
2
2
2𝑨 +𝛹I (𝑥 + ℎ, 𝑦, 𝑧, 𝑡)}+∙{
2ℎ −𝛹I (𝑥 − ℎ, 𝑦, 𝑧, 𝑡)
2𝑨 +𝛹I (𝑥, 𝑦 + ℎ, 𝑧, 𝑡)}+∙{
2ℎ −𝛹I (𝑥, 𝑦 − ℎ, 𝑧, 𝑡)
2𝑨 +𝛹I (𝑥, 𝑦, 𝑧 + ℎ, 𝑡)}+∙{
2ℎ −𝛹I (𝑥, 𝑦, 𝑧 − ℎ, 𝑡)−|𝑨|2 𝛹R (𝑡)+𝛹R (𝑡)2[ −{𝛹R (𝑡) + 𝛹I (𝑡)2 }𝛹R (𝑡) ]- 37 -(2.20)

1−𝑉(𝑡)𝛹R (𝑡 + 𝛥𝑡) + 𝑉(𝑡)𝛹I (𝑡)2+𝛹I (𝑥 + ℎ, 𝑦, 𝑧, 𝑡)
1+ 2 { −2𝛹I (𝑥, 𝑦, 𝑧, 𝑡) }ℎ+𝛹I (𝑥 − ℎ, 𝑦, 𝑧, 𝑡)+𝛹I (𝑥, 𝑦 + ℎ, 𝑧, 𝑡)
1+ 2 { −2𝛹I (𝑥, 𝑦, 𝑧, 𝑡) }ℎ+𝛹I (𝑥, 𝑦 − ℎ, 𝑧, 𝑡)+𝛹I (𝑥, 𝑦, 𝑧 + ℎ, 𝑡)
1+ 2 { −2𝛹I (𝑥, 𝑦, 𝑧, 𝑡) }ℎ+𝛹I (𝑥, 𝑦, 𝑧 − ℎ, 𝑡)
1
2𝑨 +𝛹R (𝑥 + ℎ, 𝑦, 𝑧, 𝑡 + 2 𝛥𝑡)}𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + 𝛥𝑡 + 2ℎ ∙ {1−𝛹R (𝑥 − ℎ, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2
1
2𝑨 +𝛹R (𝑥, 𝑦 + ℎ, 𝑧, 𝑡 + 2 𝛥𝑡)
}
+
∙{
1
2ℎ
−𝛹R (𝑥, 𝑦 − ℎ, 𝑧, 𝑡 + 𝛥𝑡)
2
1
2𝑨 +𝛹R (𝑥, 𝑦, 𝑧 + ℎ, 𝑡 + 2 𝛥𝑡)
}
+
∙{
1
2ℎ
−𝛹R (𝑥, 𝑦, 𝑧 − ℎ, 𝑡 + 𝛥𝑡)2−|𝑨|2 𝛹I (𝑡)+𝛹I (𝑡)2
1{𝛹− R (𝑡 + 𝛥𝑡) + 𝛹I (𝑡)2 } 𝛹I (𝑡)2[]- 38 -(2.21)

+𝑉(𝑡)𝛹I (𝑡 + 𝛥𝑡)
𝛥𝑡
)
2
1
𝛥𝑡
+ 2 −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + )
ℎ
2
𝛥𝑡
+𝛹
(𝑥
−
ℎ,
𝑦,
𝑧,
𝑡
+
)
R
{
2 }
𝛥𝑡
+𝛹R (𝑥, 𝑦 + ℎ, 𝑧, 𝑡 + )21𝛥𝑡+ 2 −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + )
ℎ
2
𝛥𝑡
+𝛹
(𝑥,
𝑦
−
ℎ,
𝑧,
𝑡
+
)
R
{
2 }
𝛥𝑡
+𝛹R (𝑥, 𝑦, 𝑧 + ℎ, 𝑡 + )
2
1
𝛥𝑡
1
1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡 + 𝛥𝑡) + 𝛥𝑡 + ℎ 2 −2𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 )22𝛥𝑡{+𝛹R (𝑥, 𝑦, 𝑧 − ℎ, 𝑡 + 2 )}2𝑨 +𝛹I (𝑥 + ℎ, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)}+∙{
2ℎ −𝛹I (𝑥 − ℎ, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2𝑨 +𝛹I (𝑥, 𝑦 + ℎ, 𝑧, 𝑡 + 𝛥𝑡)}+∙{
2ℎ −𝛹I (𝑥, 𝑦 − ℎ, 𝑧, 𝑡 + 𝛥𝑡)
2𝑨 +𝛹I (𝑥, 𝑦, 𝑧 + ℎ, 𝑡 + 𝛥𝑡)}+∙{
2ℎ −𝛹I (𝑥, 𝑦, 𝑧 − ℎ, 𝑡 + 𝛥𝑡)
1−|𝑨|2 𝛹R (𝑡 + 𝛥𝑡)2
1+𝛹R (𝑡 + 𝛥𝑡)22
11− {+𝛹R (𝑡 + 2 𝛥𝑡) } 𝛹R (𝑡 + 𝛥𝑡)
2[]+𝛹I (𝑡 + 𝛥𝑡)2
+𝛹R (𝑥 + ℎ, 𝑦, 𝑧, 𝑡 +ここでℎは空間離散幅である．

- 39 -(2.22)

2.1.3 TDGL 方程式(スカラーポテンシャル)の離散化
次に(1.60)式を離散化する．まず左辺に(2.7)式を適用する．
𝜎𝜵2 𝑉 =+𝜎 +𝑉(𝑥 + ℎ, 𝑦, 𝑧) − 2𝑉(𝑥, 𝑦, 𝑧) + 𝑉(𝑥 − ℎ, 𝑦, 𝑧)ℎ2
𝑉(𝑥, 𝑦 + ℎ, 𝑧) − 2𝑉(𝑥, 𝑦, 𝑧) + 𝑉(𝑥, 𝑦 − ℎ, 𝑧)(2.23)ℎ2
𝑉(𝑥, 𝑦, 𝑧 + ℎ ) − 2𝑉(𝑥, 𝑦, 𝑧) + 𝑉(𝑥, 𝑦, 𝑧 − ℎ )
+
{
}
ℎ2

続けて次のように変形する．
𝑉(𝑥 + ℎ, 𝑦, 𝑧) − 2𝑉 (𝑥, 𝑦, 𝑧) + 𝑉(𝑥 − ℎ, 𝑦, 𝑧)ℎ2()(
𝑉 𝑥, 𝑦 + ℎ, 𝑧 − 2𝑉 𝑥, 𝑦, 𝑧) + 𝑉(𝑥, 𝑦 − ℎ, 𝑧)𝜎 +=ℎ2𝑉(𝑥, 𝑦, 𝑧 + ℎ ) − 2𝑉 (𝑥, 𝑦, 𝑧) + 𝑉(𝑥, 𝑦, 𝑧 − ℎ )
+
{
}
ℎ2
+

3 × 2𝑉 (𝑥, 𝑦, 𝑧)ℎ2
𝑉 (𝑥 + ℎ, 𝑦, 𝑧) + 𝑉(𝑥 − ℎ, 𝑦, 𝑧)−ℎ2−𝜎𝑉 (𝑥, 𝑦 + ℎ, 𝑧) + 𝑉(𝑥, 𝑦 − ℎ, 𝑧)−ℎ2
𝑉 (𝑥, 𝑦, 𝑧 + ℎ ) + 𝑉(𝑥, 𝑦, 𝑧 − ℎ ){−}ℎ2
6𝑉(𝑥, 𝑦, 𝑧)
𝜎 −𝑉(𝑥 + ℎ, 𝑦, 𝑧) − 𝑉(𝑥 − ℎ, 𝑦, 𝑧)=− 2ℎ −𝑉(𝑥, 𝑦 + ℎ, 𝑧) − 𝑉(𝑥, 𝑦 − ℎ, 𝑧){−𝑉(𝑥, 𝑦, 𝑧 + ℎ ) − 𝑉(𝑥, 𝑦, 𝑧 − ℎ )}- 40 -(2.24)

右辺第 1 項では，次のような変形を行う．
i(𝛹 ∗ 𝜵2 𝛹 − 𝛹𝜵2 𝛹 ∗ )2i +(𝛹I + i𝛹R )𝜵2 (𝛹R + i𝛹I )}= {2 −(𝛹R + i𝛹I )𝜵2 (𝛹I + i𝛹R )1 +(i𝛹I − 𝛹R )𝜵2 (i𝛹R − 𝛹I )}= {2 −(i𝛹R − 𝛹I )𝜵2 (i𝛹I − 𝛹R )+i𝛹I 𝜵2 (i𝛹R − 𝛹I )1 −𝛹R 𝜵2 (i𝛹R − 𝛹I )=2 −i𝛹R 𝜵2 (i𝛹I − 𝛹R )2
{ +𝛹I 𝜵 (i𝛹I − 𝛹R ) }+i𝛹I 𝜵2 (i𝛹R )+i𝛹I 𝜵2 (−𝛹I )−𝛹R 𝜵2 (i𝛹R )1 −𝛹R 𝜵2 (−𝛹I )=2 −i𝛹R 𝜵2 (i𝛹I )−i𝛹R 𝜵2 (−𝛹R )+𝛹I 𝜵2 (i𝛹I ){ +𝛹I 𝜵2 (−𝛹R ) }−𝛹I 𝜵2 (𝛹R )−i𝛹I 𝜵2 (𝛹I )−i𝛹R 𝜵2 (𝛹R )1 +𝛹R 𝜵2 (𝛹I )=2 +𝛹R 𝜵2 (𝛹I )+i𝛹R 𝜵2 (𝛹R )+i𝛹I 𝜵2 (𝛹I ){ −𝛹I 𝜵2 (𝛹R ) }1 −2𝛹I 𝜵2 (𝛹R )}= {2 +2𝛹R 𝜵2 (𝛹I )= 𝛹R 𝜵2 (𝛹I ) − 𝛹I 𝜵2 (𝛹R )- 41 -(2.25)

さらに(2.7)式を適用する．
𝛹R 𝜵2 (𝛹I ) − 𝛹I 𝜵2 (𝛹R ) =+𝛹R +𝛹I (𝑥 + ℎ, 𝑦, 𝑧) − 2𝛹I (𝑥, 𝑦, 𝑧) + 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)ℎ2𝛹I (𝑥, 𝑦 + ℎ, 𝑧) − 2𝛹I (𝑥, 𝑦, 𝑧) + 𝛹I (𝑥, 𝑦 − ℎ, 𝑧)ℎ2𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) − 2𝛹I (𝑥, 𝑦, 𝑧) + 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )
+
{
}
ℎ2
+

−𝛹I +

𝛹R (𝑥 + ℎ, 𝑦, 𝑧) − 2𝛹R (𝑥, 𝑦, 𝑧) + 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)ℎ2𝛹R (𝑥, 𝑦 + ℎ, 𝑧) − 2𝛹R (𝑥, 𝑦, 𝑧) + 𝛹R (𝑥, 𝑦 − ℎ, 𝑧)ℎ2𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) − 2𝛹R (𝑥, 𝑦, 𝑧) + 𝛹R (𝑥, 𝑦, 𝑧 − ℎ ){+}ℎ2−3 × 2𝛹I (𝑥, 𝑦, 𝑧)𝛹R 𝛹I (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)= 2ℎ 𝛹I (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹I (𝑥, 𝑦 − ℎ, 𝑧){𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}−3 × 2𝛹R (𝑥, 𝑦, 𝑧)𝛹I 𝛹R (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)− 2ℎ 𝛹R (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹R (𝑥, 𝑦 − ℎ, 𝑧){𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}- 42 -(2.26)

右辺第 2 項では，次のような変形を行う．
−𝜵 ∙ (|𝛹|2 𝑨)
𝐴𝑥
𝐴𝑥
2
𝐴
= −𝛁 ∙ ( 𝑦 ) 𝛹𝑅 − 𝛁 ∙ (𝐴𝑦 ) 𝛹𝐼2
𝐴𝑧
𝐴𝑧
𝜕𝛹R2
𝜕𝛹I2
+𝐴𝑥
+𝐴𝑥
𝜕𝑥
𝜕𝑥
2
𝜕𝛹R
𝜕𝛹I2
= − +𝐴𝑦
− +𝐴𝑦
𝜕𝑦
𝜕𝑦
𝜕𝛹R2
𝜕𝛹I2
( +𝐴𝑧 𝜕𝑧 ) (+𝐴𝑧 𝜕𝑧 )𝜕𝛹R𝜕𝛹I
+2𝛹R 𝐴𝑥+2𝛹I 𝐴𝑥
𝜕𝑥
𝜕𝑥
𝜕𝛹R
𝜕𝛹I
= − +2𝛹R 𝐴𝑦− +2𝛹I 𝐴𝑦
𝜕𝑦
𝜕𝑦
𝜕𝛹R
𝜕𝛹I
( +2𝛹R 𝐴𝑧 𝜕𝑧 ) (+2𝛹I 𝐴𝑧 𝜕𝑧 )(2.27)

これに(2.4)式を適用する．
𝛹R (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)2ℎ𝛹R (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹R (𝑥, 𝑦 − ℎ, 𝑧)− +2𝛹R 𝐴𝑦
2ℎ𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹R (𝑥, 𝑦, 𝑧 − ℎ ){ +2𝛹R 𝐴𝑧}
2ℎ+2𝛹R 𝐴𝑥𝛹I (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)+2𝛹I 𝐴𝑥2ℎ𝛹I (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹I (𝑥, 𝑦 − ℎ, 𝑧)− +2𝛹I 𝐴𝑦
2ℎ𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹I (𝑥, 𝑦, 𝑧 − ℎ ){ +2𝛹I 𝐴𝑧}
2ℎ- 43 -(2.28)

まとめると，次のようになる．
6𝑉 (𝑥, 𝑦, 𝑧)
𝜎 −𝑉(𝑥 + ℎ, 𝑦, 𝑧) − 𝑉(𝑥 − ℎ, 𝑦, 𝑧)− 2=ℎ −𝑉(𝑥, 𝑦 + ℎ, 𝑧) − 𝑉(𝑥, 𝑦 − ℎ, 𝑧){−𝑉(𝑥, 𝑦, 𝑧 + ℎ ) − 𝑉(𝑥, 𝑦, 𝑧 − ℎ )}−6𝛹I (𝑥, 𝑦, 𝑧)𝛹R 𝛹I (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)ℎ 2 𝛹I (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹I (𝑥, 𝑦 − ℎ, 𝑧){𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}−6𝛹R (𝑥, 𝑦, 𝑧)𝛹I 𝛹R (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)− 2ℎ 𝛹R (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹R (𝑥, 𝑦 − ℎ, 𝑧){𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}
𝛹R (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)2ℎ𝛹R (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹R (𝑥, 𝑦 − ℎ, 𝑧)− +2𝛹R 𝐴𝑦
2ℎ𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹R (𝑥, 𝑦, 𝑧 − ℎ ){ +2𝛹R 𝐴𝑧}
2ℎ+2𝛹R 𝐴𝑥𝛹I (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)2ℎ𝛹I (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹I (𝑥, 𝑦 − ℎ, 𝑧)− +2𝛹I 𝐴𝑦
2ℎ𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹I (𝑥, 𝑦, 𝑧 − ℎ ){ +2𝛹I 𝐴𝑧}
2ℎ+2𝛹I 𝐴𝑥- 44 -(2.29)

𝜎

両辺を− 2で割る．
ℎ

6𝑉 (𝑥, 𝑦, 𝑧)−𝑉(𝑥 + ℎ, 𝑦, 𝑧) − 𝑉(𝑥 − ℎ, 𝑦, 𝑧)=−𝑉(𝑥, 𝑦 + ℎ, 𝑧) − 𝑉(𝑥, 𝑦 − ℎ, 𝑧)−𝑉(𝑥, 𝑦, 𝑧 + ℎ ) − 𝑉(𝑥, 𝑦, 𝑧 − ℎ )−6𝛹I (𝑥, 𝑦, 𝑧)𝛹R 𝛹I (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)−𝜎 𝛹I (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹I (𝑥, 𝑦 − ℎ, 𝑧){𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}−6𝛹R (𝑥, 𝑦, 𝑧)𝛹I 𝛹R (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)+𝜎 𝛹R (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹R (𝑥, 𝑦 − ℎ, 𝑧){𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}+𝛹R 𝐴𝑥 {𝛹R (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)}ℎ+ [+𝛹R 𝐴𝑦 {𝛹R (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹R (𝑥, 𝑦 − ℎ, 𝑧)}]𝜎+𝛹R 𝐴𝑧 {𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}+𝛹I 𝐴𝑥 {𝛹I (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)}ℎ+ [+𝛹I 𝐴𝑦 {𝛹I (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹I (𝑥, 𝑦 − ℎ, 𝑧)}]𝜎+𝛹I 𝐴𝑧 {𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}- 45 -(2.30)

さらに変形を続けると，最終的に次のようになる．
+𝑉(𝑥 + ℎ, 𝑦, 𝑧) + 𝑉(𝑥 − ℎ, 𝑦, 𝑧)
1
𝑉(𝑥, 𝑦, 𝑧) = {+𝑉(𝑥, 𝑦 + ℎ, 𝑧) + 𝑉(𝑥, 𝑦 − ℎ, 𝑧)}
6+𝑉(𝑥, 𝑦, 𝑧 + ℎ ) + 𝑉(𝑥, 𝑦, 𝑧 − ℎ )−6𝛹I (𝑥, 𝑦, 𝑧)𝛹R 𝛹I (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)−
6𝜎 𝛹I (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹I (𝑥, 𝑦 − ℎ, 𝑧){𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}−6𝛹R (𝑥, 𝑦, 𝑧)𝛹I 𝛹R (𝑥 + ℎ, 𝑦, 𝑧) + 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)+
6𝜎 𝛹R (𝑥, 𝑦 + ℎ, 𝑧) + 𝛹R (𝑥, 𝑦 − ℎ, 𝑧){𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) + 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}(2.31)+𝛹R 𝐴𝑥 {𝛹R (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹R (𝑥 − ℎ, 𝑦, 𝑧)}ℎ+ [+𝛹R 𝐴𝑦 {𝛹R (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹R (𝑥, 𝑦 − ℎ, 𝑧)}]6𝜎
+𝛹R 𝐴𝑧 {𝛹R (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹R (𝑥, 𝑦, 𝑧 − ℎ )}+𝛹I 𝐴𝑥 {𝛹I (𝑥 + ℎ, 𝑦, 𝑧) − 𝛹I (𝑥 − ℎ, 𝑦, 𝑧)}ℎ+ [+𝛹I 𝐴𝑦 {𝛹I (𝑥, 𝑦 + ℎ, 𝑧) − 𝛹I (𝑥, 𝑦 − ℎ, 𝑧)}]6𝜎
+𝛹I 𝐴𝑧 {𝛹I (𝑥, 𝑦, 𝑧 + ℎ ) − 𝛹I (𝑥, 𝑦, 𝑧 − ℎ )}ただしこれでは時間発展性がないため，左辺のみ𝑉(𝑥, 𝑦, 𝑧, 𝑡)を𝑉 (𝑥, 𝑦, 𝑧, 𝑡 + ∆𝑡)と置き換え
る．

2.2

𝐽cを得る方法

シミュレーションを実行することで，あらゆる条件での𝐽c を求めることが，本研究の目的である．
すなわちこれを適切により早く得られるよう工夫することは，理解を加速させるために非常に重要である．

- 46 -

2.2.1 𝑬と常伝導電流
1章で述べたTDGL方程式の変形と3章で述べる条件を用いて数値解析を実行すれば，ある𝐽についての𝐸を求めることができる．しかしそうして得られた𝐸には，磁束線が動くことによる誘導
起電力のほかに，常伝導電流の抵抗に由来する成分が含まれている．これはTDGL方程式を
クーロンゲージを用いて変形したことによる．ファイゼロゲージを用いればこの成分は含まれな
い．
常伝導電流の抵抗に由来する成分が含まれている状態では，𝐸から得られる𝐽c は不明瞭なものになる．なぜなら，𝐽c は磁束線の運動による誘導起電力の有無によって決定されるからである．
磁束線の運動に由来する成分に注目するため，常伝導電流の抵抗に由来する成分は，𝐸から取りのぞくべきである．
SCナノワイヤ全体の常伝導抵抗がわかれば，オームの法則から，その抵抗に由来する成分の量を確定することができる．しかし，SCナノワイヤ全体の常伝導抵抗は，ピン領域の常伝導
導電率，ピンの配置，境界の条件等の様々な要因によって左右されるため，数値解析の引数として与える値から直接求めることは困難であった．
そこで，実際に小さな𝐽と𝐵でシミュレーションを走らせる手法をとった．このときの𝐽と𝐵は，これまでのシミュレーションで磁束線の運動が確認されていない，0.02と0.1とした．そうして得られた𝐸は常伝導電流の抵抗のみに由来するものであると言えるから，SCナノワイヤ全体の常伝導抵抗が決定できる．あとは通常のシミュレーションで得られる𝐸から，オームの法則に従う常伝導電流の抵抗に由来する成分を引けばよい．実際にこの手法で得た𝐸-𝐽特性をFig. 2.1に示す．

- 47 -

normalized electric field 𝐸

0.1

0.05

0
0

0.1

0.2

normalized current density 𝐽

Fig. 2.1: シミュレーションによって得られた𝐸-𝐽特性(図中青線)と，その𝐸から磁束線運動に由来する成分だけを抽出した𝐸 -𝐽特性(図中赤線)をあらわすグラフ
赤線は，𝐽c を超える範囲でのグラフの凹凸を保持したまま，𝐽c を超えない範囲で𝐸を非常に小さく保っている．その値はおよそ10−4 以下と，磁束線運動に由来する成分と比較して十分に小さい．なおシミュレーショ
ンの条件は，第3章で述べる基本条件とした．

2.2.2 二分探索法の利用
磁束線運動の有無が明瞭であることは，コンピュータプログラム上でも有利である．𝐽c を求めるためにはまず𝐸 -𝐽特性を得ることが一般的であるが，シミュレーション中にはっきりと磁束線運動の有無が判定できる場合には0から𝐽c に至るまですべての𝐽を測らずとも，必要最低限の試行回数で𝐽c を求め，かかる時間を短縮することが可能である．今回はその手段として，二分探索- 48 -

法を利用する．
まず，𝐸 -𝐽特性は右肩上がりになるため，本質的にはソート済みリストであるととらえることができる．そして，特定の𝐸を与える𝐽が𝐽c だとすると，それを探索する問題であるものとして二分探索法を実施すればよい．実際は特定の値ちょうどが返ることは滅多にないため，一定の精度が得られるステップだけで探索を打ち切るものとする．シミュレータに合うよう改変を加えた二分探索法をFig. 2.2に示す．

start

0→𝐽
𝐽d → 𝑤true

𝑤 < 𝑤tfalse
true

𝐸𝐽 ≈ 𝐸t𝐽c is 𝐽.

false
𝑤
→𝑤
2

true

𝐸𝐽 < 𝐸t
false
𝐽−𝑤 →𝐽

𝐽+𝑤 →𝐽

Fig. 2.2: 改変を加えた二分探索法のフローチャート
𝐽, 𝐽d はそれぞれ𝐸の測定を行う電流密度と対破壊電流密度を示す．
𝑤, 𝑤t はそれぞれ，そのステップにおける分解能と目標分解能を示す．
𝐸𝐽 はその𝐽で測定した𝐸を示し，𝐸tは，磁束線運動が発生したと判断するための，𝐸に対する閾値である．
- 49 -

二分探索法を利用する場合は𝐽を増加させていくだけでなく，ときによっては減少させる操作も必要になるが，ここで1つ問題が生じる．𝐽を増加させつつ迎える𝐽c と，減少させつつ迎える𝐽c は異なるのである．ピン止め力が摩擦力に見立てられることがあるように，止まっている磁束線を止まったままにしておくことができたとしても，同じ力の大きさで，動いている磁束線を止めることができるとは限らない．この現象を確認するため，Fig. 2.3を用意した．

normalized electric field 𝐸
Normalized
electric field E

0.08

0.06

0.04

0.02

0
0.1

0.12

0.14

0.16

0.18

0.2

normalized current density 𝐽
Normalized
current density J
Fig. 2.3: 𝐽を増加させつつとった𝐸 -𝐽特性(図中青線)と，𝐽を減少させつつとった𝐸-𝐽特性(図中赤線)をあらわすグラフ
シミュレーション条件はどちらも第3章で述べる基本条件で同一である
が，𝐽を増加させつつとった場合の𝐽c は0.144であるのに対し，𝐽を減少させつつとった場合の𝐽c は0.130である．この差は実に0.014に及んでいる．

- 50 -

Fig. 2.3で見られる𝐽c の差は大きく，二分探索を利用する上での都合以前に，𝐽c の決定方法として不適切であるとすら言える．よって，𝐽を増加させつつとった𝐽c と，𝐽を減少させつつとった𝐽c を混在させてはならない．今回はすべての𝐽c を，𝐽を増加させつつとるものとする．
すべての𝐽c を，𝐽を増加させつつとるために，もしシミュレーション中に𝐽を減少させる必要があったときには，休み時間を挟むこととした．休み時間中に測定は行わず，𝐽を非常に小さく設定する．その値は，これまでのシミュレーションで磁束線の運動が確認されていない，0.02とした．こnormalized current density 𝐽のときの様子をFig. 2.4に示す．

測定 II
測定 I

normalized current density 𝐽(a)time 𝑡測定 III
測定 IV

休み時間
(b)

time 𝑡Fig. 2.4: シミュレーション中で𝐽を増加させるとき(a)と減少させるとき(b)の挙動の違い色のついた四角形の高さはそのときの𝐽を，横座標は時間をそれぞれ示している．また色については，青は磁束が停止中，赤は磁束が運動中，
黒は測定中につき未定であることをそれぞれ示している．

- 51 -

二分探索法を利用したやり方では，直前のステップで磁束線が運動していなければ𝐽を増加させ，逆に運動していれば減少させるが，それをFig. 2.4の(a)と(b)がそれぞれ示している．減少
させるときには休み時間を挟み，一旦確実に磁束線の運動を止める．そうすることで，𝐽を増加させながら計測していくときと同じ状況を作り出している．

主要なアルゴリズム

2.3

2.3.1 ガウス=ザイデル法
連立1次方程式の解法には様々な手段があるが，その中に反復法の名で総称される1群の方法がある．この反復法の1つにガウス＝ザイデル法がある[14, 15, 16]．
このガウス＝ザイデル法は反復法の1つであるヤコビ法を改良したもので，連立1次方程式
𝐴𝒙 = 𝒃を仮定すると，ヤコビ法では，𝑖−1

𝑛𝑥𝑖𝑘+1 = (𝑏𝑖 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘 )
𝑗=1
𝑗=𝑖+1

(2.32)と定義され，ガウス＝ザイデル法では，𝑖−1𝑛𝑎𝑖𝑖 𝑥𝑖𝑘+1 = (𝑏𝑖 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘+1 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘 )
𝑗=1
𝑗=𝑖+1

(2.33)と定義される[17]．変化値|𝑥𝑖𝑘+1 − 𝑥𝑖𝑘 |が任意の変化値以下になるまで繰り返しでこの計算を行う．ガウス＝ザイデル法はヤコビ法と比べて実装した際の使用容量と計算速度の両方の点で優れている．

- 52 -

2.3.2 オイラー法
常微分方程式の数値解法の1つにオイラー法がある[15, 17]．これは導関数の定義式
𝑑𝑥𝑓 (𝑡 + Δ𝑡) − 𝑓(𝑡)= 𝑓(𝑡, 𝑥) = lim
Δ𝑡
𝑑𝑡
Δ𝑡

(2.34)

が元となっている．ここで，Δ𝑡が十分小さいと仮定し𝑓 (𝑡, 𝑥) =𝑓(𝑡 + Δ𝑡) − 𝑓 (𝑡)Δ𝑡(2.35)のように差分商で置き換える．これより，𝑡 + Δ𝑡における変数の値𝑥(𝑡 + Δ𝑡)は𝑥 (𝑡 + Δ𝑡) = 𝑥 (𝑡) + Δ𝑡 ∙ 𝑓(𝑡, 𝑥)(2.36)と表現できる．オイラー法は数学的にシンプルであり，プログラムに実装することも容易であるが，1階段数常微分方程式の数値解法としては，𝑡の発展に伴って誤差が蓄積されるため精度が悪い．この計算誤差を低減するためにいくつかの手法が考えられている[15]．

2.3.3 二分探索法
ソート済みリストに対する探索アルゴリズムの1つとして二分探索法がある[26]．二分探索法は，探索範囲のうちの中央の値を見て，求める値との大小関係を判定し，探索範囲のうちの前半か後半かのどちらに求める値があるかを確定させながら探索範囲を狭めていく方法である．この一連のステップで探索範囲は半分以下になるため，𝑛回のステップを行えば探索範囲は0.5𝑛倍以下になる．ここで，探索範囲がちょうど半分になるとは限らないのは，見た中央の値が求める値でなかった場合に，次の探索範囲候補から除外されるからである．見た中央の値が求める値であれば当然，その時点で探索は終了する．探索範囲に含まれるデータが1つ以下となるとき探索は終了するので，リストにあるデータの数が[2𝑛 , 2𝑛+1 )個ならば，𝑛回のステップで探索を完了できる．逆に言えば，𝑁個のデータに必要なステップの回数は，高々log 2 𝑁である．このことは，特に𝑁が大きい場合に探索が短い時間で完了することを示している．このアルゴリズムのフローチャートをFig. 2.5に示す．
- 53 -

start0 → 𝑖mi𝑁 − 1 → 𝑖max

𝑖max = 𝑖mi

true

false
𝑖mi +

𝑖m x − 𝑖mi2→ 𝑖mi
The key was not found.

𝐷𝑖

=𝐾

true

false
𝐷𝑖

<𝐾

The key is 𝑖mi .
true

false
𝑖mi − 1 → 𝑖m x

𝑖mi + 1 → 𝑖min

Fig. 2.5: 二分探索法のフローチャートの一例
求めるデータの番号を key として探索を進める．𝑖min ，𝑖mid ，𝑖max はそれぞれその時点の探索範囲の最小の番号，中央の番号，最大の番号を示す．𝑁はデータの数を示す．𝐷𝑥 は番号に𝑥を持つデータの値を示す．𝐾は求めるデータの値を示す．なお，データが昇順であることを前提としている．
- 54 -

2.4

GPGPU

GPGPUとはGeneral-purpose computing on graphics processing unitsの略称で，GPUの演算資源を画像処理以外の目的に利用すること，またその技術を指す．定型な大量の演算を並列処理によって高速に行うことができ，プログラマブルシェーダの発展によって柔軟性を獲得した
GPUを他の計算に応用することを目的に様々な開発環境や製品が登場してきている．
GPU はメモリに連続的にアクセスし，かつ条件分岐のない密な計算に用いるのに適している．逆に条件分岐の多数存在する処理，木構造やポインタをたどる処理を含む処理は苦手としている．これは，GPUが複数の演算ユニットをまとめてクラスタとしており，演算ユニットに命令を出すインストラクションユニットは多くのGPUでクラスタごとに1つずつしか設置されておらず，同一クラスタ内のすべてのプロセッサが異なるデータを受け取り，同じ命令を処理するためである．
このようなSMID (Single-Instruction Multiple-Data)型のデータ処理は画像処理などの単純で大量な計算処理を得意とする一方で，条件分岐を含む体系ではオーバーヘッドがかさみ，効率を著しく落としてしまう．CPUはこのようなペナルティを避けるために，プリフェッチ，プリデコード，投機実行，あるいはレジスタリネーミングといった機能を持つが，本来画像処理のみを目的として製造されるGPUには一般的に搭載されていない．

- 55 -

また，単精度浮動小数点演算では高いパフォーマンスを発揮するGPUだが，画像処理において必要とされることの少ない倍精度浮動小数点演算ではその限りではない．これは，倍精度浮動小数点演算を行う際に，単精度用に作られた演算器で代用しなければならないため，対象データの分割や複数回の演算といった余分なプロセスが必要になるからである．倍精度専用演算器を搭載した製品も存在するが，倍精度専用演算器では単精度演算が不可能となるため，倍精度性能と単精度性能はトレードオフの関係にあるといえる．
さらに，GPGPUを用いてプログラムを設計する場合は，極力メモリへのアクセスを連続にする，共有メモリを利用する場合はそれを用いるデータを同一ストリーミングマルチプロセッサ内に格納する，条件分岐をできるだけ削減する，データ構造は基本的に配列のみを使うなど，性能を十分に発揮させるための制約が多く存在する．加えて，デバイスとの通信を行うローレベルのAPIを使う必要があるなどの点から，プログラミングの難易度が高いと言える．
上記のように問題点も少なからず抱えるGPGPUだが，実用ソフトウェアも次第に数を増やしており，開発環境は整いつつある．方程式を莫大な回数繰り返す今回の数値解析では，このGPGPUを用いた．超伝導体内の離散化された空間の1点ごとに1つずつのスレッドを当てはめている．詳細は付録6.1節で述べる．

- 56 -

第3章

計算条件

SCナノワイヤについての仮定

3.1
3.1.1

電流密度と外部磁束密度の条件

SCナノワイヤに与える電流密度𝑱と磁束密度𝑩はそれぞれ以下2式のように定義する．
0
𝑱 = (𝐽ext )00
𝑩 = ( 𝐵ext si 𝜃𝐵 )𝐵ext cos 𝜃𝐵

(3.1)

(3.2)

ここで𝐽ext , 𝐵ext はそれぞれ𝑱, 𝑩の大きさであり，𝜃𝐵 は𝑩の𝑧軸に対する角度である．𝜃𝐵 = 0の状態
では，𝑱に対して𝑩を垂直に与えることになるので，横磁界である．一方で𝜃𝐵 = 90の状態では，𝑱に対して𝑩を平行に与えることになるので，縦磁界である．

3.1.2

オーダーパラメータとスカラーポテンシャルの初期条件

オーダーパラメータ𝛹とスカラーポテンシャル𝑉の初期条件は以下のように定義する．
𝛹(𝑡 = 0) = cos𝜃 + i ∙ si 𝜃𝑉(𝑡 = 0) = −𝐽𝑦 ∙ 𝑦
𝜎

(3.3)
(3.4)

ここで，𝜃は𝛹の位相を意味しており，初期値は0 ≤ 𝜃 < 2πの範囲で乱数によって決定される．

3.1.3

境界条件

真空と超伝導体の境界について，超伝導体から真空へ電流が流出・流入しない，つまりは境界面を横切って電流が流れないという条件を与える．ここで，超伝導電流密度𝐽s の𝑥成分について考慮し，𝑥方向の単位面ベクトル𝒏𝑥 との内積を0と置くことでこの条件を求める．

- 57 -

(1.32)式より，まず次のようになる．
i𝐽s = (𝛹 ∗ 𝜵𝛹 − 𝛹𝜵𝛹 ∗ ) − |𝛹|2 𝑨
2

(3.5)

𝛹 = 𝛹R + i𝛹I ，𝛹 ∗ = 𝛹R − i𝛹I として，𝑥成分において展開すると，次のようになる．
i{ (𝛹 ∗ 𝜵𝛹 − 𝛹𝜵𝛹 ∗ ) − |𝛹|2 𝑨}|
2𝒙𝜕
i +(𝛹R − i𝛹I ) 𝜕𝑥 (𝛹R + i𝛹I )} − (𝛹R2 + 𝛹I2 ) ∙ 𝐴𝑥
= {𝜕2−(𝛹R + i𝛹I ) (𝛹R − i𝛹I )𝜕𝑥𝜕𝜕+𝛹R 𝛹R + i𝛹R
𝛹
𝜕𝑥
𝜕𝑥 I
𝜕
𝜕
i −i𝛹I 𝜕𝑥 𝛹R + 𝛹I 𝜕𝑥 𝛹I=− 𝛹R2 𝐴𝑥 − 𝛹I2 𝐴𝑥𝜕𝜕2−𝛹R 𝛹R + i𝛹R
𝛹
𝜕𝑥
𝜕𝑥 I
𝜕
𝜕
−i𝛹
𝛹
−
𝛹
𝛹
I
R
I
(
𝜕𝑥
𝜕𝑥 I )
= 𝛹I

(3.6)

𝜕
𝜕
𝛹R − 𝛹R𝛹 − 𝛹R2 𝐴𝑥 − 𝛹I2 𝐴𝑥
𝜕𝑥
𝜕𝑥 I= 𝛹I (

𝜕𝜕𝛹R − 𝛹I 𝐴𝑥 ) − 𝛹R ( 𝛹I + 𝛹R 𝐴𝑥 )𝜕𝑥𝜕𝑥(2.6)式が0になるには，以下の条件を満たせば良い．
∂
𝛹 − 𝛹I 𝐴𝑥 = 0
𝜕𝑥 R

(3.7)

𝜕
𝛹 + 𝛹R 𝐴𝑥 = 0
𝜕𝑥 I

(3.8)

(2.7)及び(2.8)式をまとめると，次のようになる．
𝜕𝛹 + i𝐴𝑥 𝛹 = 0
𝜕𝑥

(3.9)

超伝導体と真空のすべての境界は，(2.9)式を元に同様の手順で求められる．まとめると
𝜵𝛹 + i𝑨𝛹 = 0(3.10)となる．
次に，スカラーポテンシャル𝑉についての境界条件を定義する．𝑦方向の境界面において，以下の条件を定義する．
- 58 -

𝜵𝑉 = −𝑱/𝜎

(3.11)(2.11)式は，超伝導体を導体で𝑦軸方向に挟んで電流を流すことを意味している．

3.2

シミュレーションモデル

本研究ではFig. 2.2 に示すシミュレーションモデルを基本とする．座標系は直角直交座標(𝑥, 𝑦, 𝑧)で示している．一辺が規格化した長さにおいて10の立方体に超伝導体が満ちていると定義している．この立方体空間の空間離散幅ℎは0.2で定義している．つまりメッシュ数は(

10 3) = 503 = 125000
0.2(3.12)

である．Fig. 3.1の立方体空間の内側に様々な条件を持ったピンを導入して計算を行う．

(0, 10, 10)(10, 10, 10)(0, 10, 0)(0, 0, 10)(10, 10, 0)𝑦(10, 0, 10)𝑧(0, 0, 0)(10, 0, 0)𝑥Fig. 3.1: 計算するシミュレーションモデルの立方体空間- 59 -

3.2.1

柱状ピンの形状

本研究では，解析するSCナノワイヤに比べて十分な長さをもつ，半径𝑟の円柱の形状をした
ピン(柱状ピン)を主に扱う．この円柱の断面の円の中心の集合からなる直線を中心線と呼ぶ．
ただし数値解析に用いる空間は離散的であるから，柱状ピンは厳密には円柱の形状をしていない．中心線からの距離が𝑟以下である要素がピンを構築するとして，円柱に近い形状を再現している．以降の図において，白い枠内にある，白(または明るい灰色)で着色された領域が，ピンの領域である．

3.2.2

基本条件

すべての異方性パラメータを1とし，𝑩 = (0, 0, 0.4)とし，𝑱を𝑦軸と平行に流し，半径𝑟 = 0.5の柱状ピンを𝑧軸と平行に4本設置した状態を基本条件と呼ぶこととする．このときのピンの中心線
はTable 3.1に示すような点の集合である．
Table 3.1: 基本条件におけるピンの中心線の座標(𝑥, 𝑦, 𝑧)座標
ピンA

(3, 3, ∀𝑧 ∈ ℝ)

ピンB

(3, 7, ∀𝑧 ∈ ℝ)

ピンC

(7, 3, ∀𝑧 ∈ ℝ)

ピンD

(7, 7, ∀𝑧 ∈ ℝ)

基本条件では，柱状ピンが磁束線に平行である上に，先行研究[34]で報告された，マッチン
グ機構によるピーク効果が起こるのに適したピン間隔を設定しているため，𝐽c は高くなる．これが再現するのは，柱状ピンにとって理想的な状況である．なお柱状ピンの半径𝑟 = 0.5は，重イオ
ン照射によって形成される人口ピンの半径に則っている[30–32]．

- 60 -

3.2.3

すべてのピン同士が平行な配置で磁束線に対して角度をとる場合
まずは超伝導体内の柱状ピンの向きと外部磁界の向きが合わなくなる状態を想定して，すべてのピン同士が平行な配置のまま，磁束線に対して様々な角度をとる場合を考える．柱状ピンの中心線は常に𝑦-𝑧平面に含まれているとし，中心線と𝑧軸のなす角を𝜃p とする．中心線が𝑧軸と平行となる状態を𝜃p = 0°とし，𝑦軸と平行となる状態を𝜃p = 90°とする．その他の条件はTable
3.2に示す．なおこの条件では𝜃p = 0°においてピーク効果により，𝐵 = 0.4としたときに𝐽c が極大となることがわかっている．ピンの外観はTable 3.3に示す．

Table 3.2: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合(ピンの角度を変更する場合)の細かな条件の一覧
ピンの本数

4

ピンの半径

0.5

ピンの中心線 (𝜃p = 0°)(3, 3, ∀𝑧 ∈ ℝ), (3, 7, ∀𝑧 ∈ ℝ), (7, 3, ∀𝑧 ∈ ℝ), (7, 7, ∀𝑧 ∈ ℝ)

ピンの中心線 (𝜃p = 90°)(3, ∀𝑦 ∈ ℝ, 7), (3, ∀𝑦 ∈ ℝ, 3), (7, ∀𝑦 ∈ ℝ, 7), (7, ∀𝑦 ∈ ℝ, 3)

外部磁界𝐵の強さ

0.4- 61 -

Table 3.3: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合の，いくつかの𝜃p における柱状ピンの外観
右図は3次元的に見たものであり，左図は真横(𝑥 > 10, 𝑦 = 0, 𝑧 = 0)か
ら見たものである．暗い色をしたピンの方が，明るい色をしたピンよりも𝑥座標が小さい．画面内の上向きが𝑦軸の正の向きであり，画面内の左向きが𝑧軸の正の向きである．
𝜃p = 0°𝜃p = 30°

𝜃p = 60°- 62 -

𝜃p = 90°

さらに，すべてのピン同士が平行な配置のまま磁束線に対してさまざまな角度をとる場合のもう1つのパターンとして，外部磁界の向きのみを変更するやり方も考えられる．外部磁界の向きは常に𝑦-𝑧平面に含まれているとし，その向きと𝑧軸のなす角を𝜃とする．𝑧軸と平行となる状態を
𝜃p = 0°とし，𝑦軸と平行となる状態を𝜃p = 90°とする．その他の条件はTable 3.4に示す．なおこの条件でも𝜃p = 0°のときピーク効果により，𝐵 = 0.4としたときに𝐽c が極大となる．

Table 3.4: すべてのピン同士が平行な配置で磁束線に対して角度をとる場合(外部磁界の角度を変更する場合)の細かな条件の一覧
ピンの本数

4

ピンの半径

0.5

ピンの中心線(3, 3, ∀𝑧 ∈ ℝ), (3, 7, ∀𝑧 ∈ ℝ), (7, 3, ∀𝑧 ∈ ℝ), (7, 7, ∀𝑧 ∈ ℝ)

外部磁界𝐵の強さ

0.4- 63 -

3.2.4

スプレーピンの場合

次に，第2項と同じように磁束線に対してピンが特定の角度で交わりつつも，ピン同士が平行でない場合を考える．柱状ピンの中心線は常に𝑦-𝑧平面に含まれているとし，中心線と𝑧軸のなす角を𝜃p とする．中心線が𝑧軸と平行となる状態を𝜃p = 0°とし，𝑦軸と平行となる状態を𝜃p =±90°とする．偶数本の柱状ピンを2グループに分け，一方の𝜃p がある値𝜃𝑒 をとるとき，もう一方のグループの𝜃p は−𝜃𝑒 をとる．このようなピンの組み合わせを，スプレーピン(splayed pins)と呼ぶこととする．その他の条件はTable 3.5に示す．なおこの条件でも𝜃p = 0°のときピーク効果により，𝐵 = 0.4としたときに𝐽c が極大となる．ピンの外観はTable 3.6に示す．

Table 3.5: スプレーピンの場合の細かな条件の一覧
ピンの本数

4

ピンの半径

0.5

ピンの中心線 (𝜃p = 0°)(3, 3, ∀𝑧 ∈ ℝ), (3, 7, ∀𝑧 ∈ ℝ), (7, 3, ∀𝑧 ∈ ℝ), (7, 7, ∀𝑧 ∈ ℝ)

ピンの中心線 (𝜃p = ±90°)(3, ∀𝑦 ∈ ℝ, 7), (3, ∀𝑦 ∈ ℝ, 3), (7, ∀𝑦 ∈ ℝ, 7), (7, ∀𝑦 ∈ ℝ, 3)

外部磁界𝐵の強さ

0.4- 64 -

Table 3.6: いくつかの𝜃p におけるスプレーピンの外観
右図は3次元的に見たものであり，左図は真横(𝑥 > 10, 𝑦 = 0, 𝑧 = 0)か
ら見たものである．暗い色をしたピンの方が，明るい色をしたピンよりも𝑥座標が小さい．画面内の上向きが𝑦軸の正の向きであり，画面内の左向きが𝑧軸の正の向きである．
𝜃p = 0°𝜃p = ±30°𝜃p = ±45°- 65 -

𝜃p = ±60°

𝜃p = ±90°- 66 -

3.2.5星状ピン

最後に，外部磁界の向きに大きく影響されないようなピンを提案する．提案するピンは，柱状
ピンを放射状に配置するもので，星状ピンと呼ぶこととする．𝑦-𝑧平面上に，60°ずつ中心線の向きが異なる3本の柱状ピンを交わらせる．このとき交点は(𝑦, 𝑧) = (5, 5)上に1つだけあるとする．
また，柱状ピンの半径は0.3とする．これを𝑥 = 3, 7の2平面に1セットずつの合計2セット配置する．
𝑥 = 3の平面上にある方のセットには，𝑧軸に平行な中心線をもつ柱状ピンが含まれるように向きを調整する．一方で𝑥 = 7の平面上にある方のセットには，𝑦軸に平行な中心線をもつ柱状ピンが含まれるように向きを調整する．星状ピンの外観をFig. 3.2に示す．

Fig. 3.2: 星状ピンの外観
右図は3次元的に見たものであり，左図は真横(𝑥 > 10, 𝑦 = 0, 𝑧 = 0)か
ら見たものである．暗い色をしたピンの方が，明るい色をしたピンよりも𝑥座標が小さい．

- 67 -

第4章

計算結果および考察

基本条件

4.1

まずは基本条件で得られた，異方性の強さごとの𝐽c -𝐵特性をFig. 4.1に示す．先行研究[34]で示されたとおり，この条件では𝐵 = 0.4に大きなピークがある．また，異方性の強さによる違いはほとんどない．このことから，𝑧軸異方性があったとしてもピーク効果がはたらきうると確認する

1

normalized critical current density 𝐽c ((3)

2 2

3

0.16

√2 𝐽

𝐽→𝐽)

ことができた．

0.11

0.06

0.1

―――
―――
―――

𝛾𝑧 = 1𝛾𝑧 = 2𝛾𝑧 = 4―――𝛾𝑧 = 8

0.2

0.3

0.4

normalized magnetic field 𝐵 (

0.5
𝐵
𝐵c2

0.6

→ 𝐵)

Fig. 4.1: 基本条件における異方性の強さごとの𝐽c -𝐵特性

- 68 -

すべてのピン同士が平行な配置で磁束線に対して

4.2角度をとる場合
すべてのピン同士が平行な配置で磁束線に対して様々な角度をとった場合の，𝐽c -𝜃p 特性を
示す．

4.2.1

超伝導体に対するピンの角度を変更する場合

超伝導体に対してピンの角度を変更した場合の，𝐽c -𝜃p 特性をFig. 4.2に示す． 𝐵 = 0.4の場

1

normalized critical current density 𝐽c ((3)

3

2 2

√2 𝐽

𝐽→𝐽)

合の0°あるいは90°に注目したグラフをFig. 4.3に示す．

0.2

―――
―――
―――

𝛾𝑧 = 1𝛾𝑧 = 2𝛾𝑧 = 4―――𝛾𝑧 = 80.1

0dashed line

𝐵 = 0.1

solid line

𝐵 = 0.4

30

60

90

angle of pins 𝜃p [°]Fig. 4.2: すべてのピン同士が平行な配置のまま磁束線に対して様々な角度をとった場合の𝐽c -𝜃p 特性- 69 -𝐽→𝐽)● 𝜃p = 0°○ 𝜃p = 90°

√2 𝐽

1
3

3

normalized critical current density 𝐽c (( )

2 2

0.2

0.15

0.1

0.05

0.91

2

3

4

5 6 7 8 910

𝑧-axis anisotropy parameter 𝛾𝑧

Fig. 4.3: すべてのピン同士が平行な配置のまま磁束線に対して0°ある
いは90°をとった場合の𝐵 = 0.4における𝐽c -𝛾𝑧 特性

まずは𝐵 = 0.4から見ていく．磁束線と柱状ピンが平行となる𝜃p = 0°で𝐽c は最大となった．このことは予想された柱状ピンのはたらきに一致する．その後𝜃p = 30°まで𝐽c は急激に落ちこみ，
それ以降ではほぼ平坦が続いている．ただし𝜃p = [60, 90]°の範囲では，異方性が小さい方が𝐽c が高い．この理由としては，𝑧軸についての異方性が強い場合には，磁束線の𝑧軸方向へのつながりが弱いということが挙げられる．これは異方性そのものの考え方から説明できることだが，ある軸への異方性が強いということは，その軸への超伝導電子対の動きが制限されているということである．結果として磁束線も、その軸方向へのつながりが弱くなるのである．ピンの角度が大きく，磁束との交点が数点しかないとき，このつながりの弱さゆえに，ピンとの交点から遠い箇所にはたらくLorentz力を支えられない．そうして磁束線の一部が流れてしまうと，なし崩し的に残り- 70 -

の部分も流れてしまう．このときの様子をFig. 4.4に示す．𝜃p が大きい場合に異方性の強さの影
響を強く受けることは，Fig. 4.3からも見て取れる．

𝑧

𝑧

𝑥

𝑥

Fig. 4.4: 横磁界に縦のピンが入っている場合の磁束線の変形
左 は 𝛾𝑧 = 1 で 右 は 𝛾𝑧 = 4 の 場 合 の ， 同 じ 条 件 で 同 じ だ け の 時 間(300sec)が経ったあとのようすである．どちらもシミュレーション領域を真上(𝑥 = 0, 𝑦 > 10, 𝑧 = 0) から見たものである．着彩された箇所が磁束線の侵入領域を示すが，𝛾𝑧 = 4の場合のほうが大きく湾曲している．

つぎに𝐵 = 0.1について見ていく．こちらは逆に，角度が大きいほど𝐽c が高い．これは𝐵 = 0.1で作用するLorentz力は十分に小さいため，ピンが磁束線に沿っておらず数点でしか交差していないとしても，それだけのピン力で釣り合うことができることからきていると考えられる．むしろ角度が大きければSCナノワイヤの上から下までピンが存在することになるので，磁束線を取り逃す可能性が低くなり，結果として𝐽c は高くなる．また，𝜃p = [30, 70]°の範囲で𝐵 = 0.4の場合に見られたように，異方性が小さい方が𝐽c が高い．この傾向の理由も，𝐵 = 0.4の場合と同じだと考えられる．

- 71 -

4.2.2ピンと磁束線の鎖交体積からの考察

ピン力は，ピンの領域と磁束線の領域が重なっている箇所の体積(鎖交体積)に比例する．そnormalized critical current density 𝐽c ((3)

2 2

1

10

0.1

0
0

5

30

60

volume of overlapping area [𝜉 3 ]

3

0.2

√2 𝐽

𝐽→𝐽)

こで，鎖交体積と𝐽c の大きさの関係を確かめるため，Fig. 4.5を示す．

090angle of pins 𝜃p [°]Fig. 4.5: ピンと磁束線の角度ごとの鎖交体積
プロットはFig. 4.2に示した結果のうち，𝐵 = 0.4の場合のものである．赤いプロットは𝛾𝑧 = 1の場合で青いプロットは𝛾𝑧 = 8の場合である．破線はピンと磁束線が1本ずつある場合の鎖交体積の計算結果である．実線
はピンが複数本ある場合の鎖交体積の合計の計算結果である．プロッ
トは左軸を，線は右軸を参照する．両軸間の比は，値が基本条件(𝜃p =0°)で一致するように調整を施してある．

- 72 -

Fig. 4.5は，Fig. 4.2に示した結果に，その𝜃p における磁束線とピンの鎖交体積を重ねたものである．鎖交体積を求める過程は付録6.2節を参照されたい．ここでピンが複数本ある場合というのは，3.2.3項で定義したように，特定の𝑦-𝑧平面に2本のピンが存在する場合を言う．つまりFig.
4.5の破線よりも実線の方がシミュレーションの状況に近い．ピンが2本ある場合，一定角度を超えると，磁束線は2本目のピンとも鎖交を始める．そのため鎖交体積の合計は，ピンが1本である
場合のように単調な減少とならない．なお，鎖交体積を求める過程で，磁束線とピンは変形することがないとしている．Fig. 4.5では左右両軸間の比を，値が基本条件(𝜃p = 0°)で一致するよ
うに調整してあるが，𝜃p = 0°以外にも40°付近で，ピンが2本の場合の鎖交体積を示す実線と，𝐽c を示すプロットが接近している．逆にそれ以外の角度では𝐽c を示すプロットの方が高い．この理由について考察する．
まず𝜃p = [15, 25]°の領域に注目する．この領域では，𝐽c を示すプロットの方が鎖交体積を示す実線よりも高いというよりも，𝐽c の降下のおとずれが遅れているように見える．この理由として，
磁束線がFig. 4.6のようにピンに沿って曲がることが挙げられる． この現象を，磁束線のピンへの追従と呼ぶこととする．

𝑦𝑧Fig. 4.6: 磁束線がピンに沿うようす
シミュレーション領域を真横(𝑥 > 10, 𝑦 = 0, 𝑧 = 0)から見たものである．着彩された箇所が磁束線の侵入領域を示す．磁束線に包まれているため直接は見えないが， 2本のピンが存在する．𝛾𝑧 = 1, 𝐽 = 0.128,
𝐵 = 0.4, 𝜃p = 20°．
- 73 -

磁束線は付近の常伝導領域に吸い込まれるような動きを見せるため，ピンの角度と外部磁界の角度の差を打ち消すように，磁束線は外部磁界の角度とは異なる角度で超伝導体に侵入することができる．これには限界があるが，そうした追従が起こった結果，𝐽c の降下のおとずれが遅れたと考えられる．
つぎに𝜃p = [45, 90]°の領域に注目する．ここでは，𝐽c は平坦であるのに対し，鎖交体積の計算結果は下がり続けているという違いがある．この理由として，十分に小さなLorentz力によって動かされる磁束線は，小さなピンによっても止められるということが挙げられる．該当領域において𝐽c は0.1程度であるので，𝑱と𝑩の積によって表されるLorentz力は小さい．このことから45°を超えるような極端な角度のずれでは，𝐽c は特定の最低ラインを切ることはない．

4.2.3 球状ピンによる考察
球の形状をもったピンを球状ピンと呼ぶ．磁束線と柱状ピンが垂直に交わっている場合はその小さな鎖交体積からピン力が生まれる．そうだとすると，柱状ピンを構成する領域のうち，磁束線と接触していない領域を超伝導領域に置き換えたとしても，ピン力は同じになるはずである．
Fig. 4.7，Table 4.1に示すような8つの球状ピンを用意する．𝐵 = 0.4で得られた結果をFig.
4.8に示す．

- 74 -

Fig. 4.7: 球状ピンの外観

Table 4.1: 8つの球状ピンの場合の細かな条件の一覧(𝑥, 𝑦, 𝑧)座標
ピンA

(3, 3, 3)

ピンB(3, 3, 7)

ピンC(3, 7, 3)

ピンD(3, 7, 7)

ピンE(7, 3, 3)

ピンF(7, 3, 7)

ピンG(7, 7, 3)

ピンH(7, 7, 7)

半径0.5- 75 -

𝐽→𝐽)
1

√2 𝐽

● sphere pins
○ cylinder pins parallel to 𝑦-axis

2 2normalized critical current density 𝐽c ((3)

3

0.15

0.1

0.05
0.91

2

3

4

5 6 7 8 910

𝑧-axis anisotropy parameter 𝛾𝑧

Fig. 4.8: 球状ピンによる𝐽c の𝛾𝑧 依存性
黒いプロットはTable 4.1の球状ピンでの結果，白いプロットはTable 3.3の
𝜃p = 90°での結果を示す．点線はそれぞれのフィッティング結果である．

Fig. 4.8でフィッティングを行っているが，𝛾𝑧 = 8における球状ピンでの結果だけは，フィッティ
ングに含めていない．フィッティング結果の傾きは，2本ともほぼ同じで，𝛾𝑧 が倍になるにつき𝐽c は0.08ずつ減少している．このことは，柱状ピンでも球状ピンでも，異方性の強さから受ける影響がほぼ同じであることを表している．すなわち4.1.1で述べたような，磁束線のつながりの弱さと異方性の強さの関係の予想が，ある程度信頼できると言える．また，フィッティング結果の切片を比較すると，球状ピンよりも柱状ピンの方が，0.06ほど小さい．これは，ピンが存在しない隙間が開いていることで，ピン力の強さ以前に，磁束線がピンと接触しないことがあるためだと考えられ
る．ピンの隙間が埋まるほど𝐽c が高くなることは，弱いLorentz力がはたらく状態においては柱状- 76 -

ピンの𝑧軸に対する角度が90°に近いほど𝐽c が高くなる理由とも一致する．
最後に，球状ピンの結果が𝛾𝑧 = 4と𝛾𝑧 = 8で変化しなかった理由だが，これはマッチング機
構によるピーク効果から来るものだと考えられる．𝑦軸に平行に配置されたピンは，磁束線を固定する𝑦座標を特定しない．つまり，磁束線は𝑦軸に平行に配置されたピンと同一の𝑥座標にとどまるが，どの𝑦座標にとどまるかまではわからない．一方で球状ピンは，磁束線を固定する𝑥座標と𝑦座標をどちらも特定する．こうすることで磁束線同士の間隔が一定に保たれ，先行研究[34]で報告されたようなマッチング機構によるピーク効果がはたらいたと考えられる．

4.2.4 超伝導体に対する外部磁界の角度を変更する場合
超伝導体に対して外部磁界の角度を変更した場合の，𝐽c -𝜃B 特性をFig. 4.9に示す．𝐵 = 0.4

1
2 2

normalized critical current density 𝐽c ((3)

3

0.2

√2 𝐽

𝐽→𝐽)

の場合の20°あるいは50°に注目したグラフをFig. 4.10に示す．

0.1

0

0

―――
―――

𝛾𝑧 = 1𝛾𝑧 = 2――――――𝛾𝑧 = 4𝛾𝑧 = 830

60

90

angle of magnetic field 𝜃[°]Fig. 4.9: 基本条件のピン配置に対して，磁束線がさまざまな角度をとった場合の𝐽c -𝜃B 特性- 77 -

𝐽→𝐽)
1

√2 𝐽

● 𝜃B = 20°○ 𝜃B = 50°

normalized critical current density 𝐽c ((3)

2 2

3

0.2

0.15

0.1

0.05

0.91

2

3

4

5 6 7 8 910

𝑧-axis anisotropy parameter 𝛾𝑧

Fig. 4.10: 基本条件のピン配置に対して磁束線が20°あるいは50°をとった場合の𝐵 = 0.4における𝐽c -𝛾𝑧 特性

𝜃B = [10, 30)°の範囲では異方性が強いほど𝐽c が高く，𝜃B = [30, 75)°の範囲では異方性が弱いほど𝐽c が高く，𝜃B = [75, 85]°の範囲では異方性が強いほど𝐽c が高い．
まず𝜃𝐵 = [10, 30)°の範囲で異方性が強いほど𝐽c が高くなる理由を考える．この範囲における異方性の強さごとの磁束線のふるまいを見比べると，同一の𝜃B であっても，異方性が弱い場合には磁束線がピンに追従できず，反対に異方性が強い場合には，Fig. 4.6に見られるように磁束線がピンに追従できている状況があるとわかる．すなわち磁束線のピンへの追従は，異方性が強いほどより強く作用するということが，磁束線運動の観察から確認できたということになる．これは先述した，𝑧軸異方性が強い場合の磁束線の𝑧軸方向へのつながりの弱さが影響したと考えられる．
- 78 -しかしそれだけでは，外部磁界ではなくピンの角度を変更した場合に𝜃p = [10, 30)°の範囲
で異方性が強いほど𝐽c が高くなるわけではないことに矛盾が生じる．その理由は，境界にある．
𝜃p > 20°においてピンは𝑦 = 0, 𝑦 = 10の境界に触れる．3.2.3項のように外部磁界が𝑧軸に平行に印加されているモデルの場合，プログラム上，𝑦 = 0, 𝑦 = 10の境界から磁束線は出入りできない．つまり磁束線のピンへの追従の強弱にかかわらず，境界で磁束線はピンから脱落する．これにより，ピンの角度を変更した場合には異方性との相関があらわれなかった．
つづいて𝜃B = [30, 75)°の範囲で異方性が弱いほど𝐽c が高くなる理由を考える．ピンの配置(Table 3.3 (𝜃p = 0°))から，4本のうち2本のピンは同一の𝑦-𝑧平面に含まれる．𝜃B が一定以上ある場合，この2本のピンの間をななめに橋渡しするかのように磁束線が大きく曲がる現象が，Fig.
4.11やFig. 4.12のように見られた．このときFig. 4.4で見たような状況に近くなり，結果として同じよ
うに異方性が弱いほど𝐽c が高くなったと考えられる．

𝑦𝑧Fig. 4.11: 磁束線によるピンの橋渡し(磁束線が𝑧 = 0の面から侵入し，
𝑧 = 10の面から出ていくとき)
シミュレーション領域を真横 (𝑥 > 10, 𝑦 = 0, 𝑧 = 0) から見たものである．
着彩された箇所が磁束線の侵入領域を示す．磁束線がピンの間を橋渡しするかのように，ゆるやかにS字を描いている．磁束線に包まれているため見えないが，𝑧軸に平行に2本のピンが存在する．𝛾𝑧 = 1, 𝐽 =0.1, 𝐵 = 0.1, 𝜃B = 45°．
- 79 -

𝑦𝑧Fig. 4.12: 磁束線によるピンの橋渡し(磁束線が𝑦 = 0の面から侵入し，
𝑦 = 10の面から出ていくとき)
シミュレーション領域を真横 (𝑥 > 10, 𝑦 = 0, 𝑧 = 0) から見たものである．
着彩された箇所が磁束線の侵入領域を示す．磁束線がピンの間を橋渡しするかのようにS字を描いている．磁束線に包まれているため見えな
い が ， 𝑧 軸 に 平 行 に 2 本 の ピ ン が 存 在 す る ． 𝛾𝑧 = 4, 𝐽 = 0.1, 𝐵 =0.1, 𝜃B = 45°．

最後に，𝜃B = [75, 85]°の範囲で異方性が強いほど𝐽c が高い理由を考える．この範囲では磁束線の断面は，ほぼ𝑧 -𝑥平面に平行となる．ここで𝑧軸についての異方性が強い場合には，Fig.
4.13に示すように断面は正円ではなく𝑧軸方向につぶれたように楕円を描く．その結果，磁束線の体積そのものが小さくなり，それにともなってLorentz力も弱まったと考えらえれる．なお，𝑧軸の異方性が強く𝑧軸方向への磁束線のつながりが弱くなったとしても，𝑦軸方向への磁束線のつながりが弱くなるわけではない．よってこの場合は，磁束線のつながりの弱さによる𝐽c の低下はほとんどない．

- 80 -

𝑧𝑥
Fig. 4.13: 楕円を描く磁束線の断面シミュレーション領域を真上 (𝑥 = 0, 𝑦 > 10, 𝑧 = 0) から見たものである．
着彩された箇所が磁束線の侵入領域を示す．磁束線の断面が，等方的な場合に見られるような円ではなく，𝑧軸方向に短軸をもつ楕円を描
いている．

4.3

スプレーピンの場合

スプレーピンの𝐽c -𝜃p 特性をFig. 4.14に示す．こちらの場合では，異方性が弱いほど𝜃 =[30, 60]°の範囲で，高い𝐽c が見られた．

- 81 -

𝐽→𝐽)
1

√2 𝐽

2 2

normalized critical current density 𝐽c ((3)

3

0.2

0.1
―――
―――
―――
―――

0

0

𝛾𝑧 = 1𝛾𝑧 = 2𝛾𝑧 = 4𝛾𝑧 = 830

60

90

angle of pins 𝜃p [°]Fig. 4.14: スプレーピンの𝐽c -𝜃p 特性

𝜃p = [30, 60]°の範囲で，ピンはシミュレーションの領域内で互いに交差している．この状態
で高い𝐽c が見られたということは，磁束線が複数の角度の異なるピンをまたぐようにして止められたということを意味している．そのときのようすはFig. 4.4と同様である．
このことから，ピンが磁束線と平行でなくても，ピン同士が平行でなくても，場合によってはピン止めに貢献できるということが言える．ただし，異方性の強さすなわち磁束線のつながりの弱さ次第でその作用は弱まってしまう．
しかしこの結果は，異方性が強く磁束線の柔軟性が高ければピンへの追従が強くあらわれるという，ここまでの考察に反している．その理由は，磁束線同士の斥力にある．スプレーピンに複数本の磁束線が追従するとき，ピンとピンの交点で磁束線も接近してしまう．この状況は，ピン同士が常に平行な場合にはありえなかった．磁束線はお互いに反発しあうので，ピン同士が交- 82 -

わってしまう以上，磁束線の追従は不完全なものとなる．

星状ピン

4.4

1

normalized critical current density 𝐽c ((3)

2 2

3

0.2

√2 𝐽

𝐽→𝐽)

星状ピンの𝐽c -𝜃p 特性はFig. 4.15に示すとおりである．𝜃p によって𝐽c は大きく変化しなかった．

0.1

0

0

―――
―――

𝛾𝑧 = 1𝛾𝑧 = 2――――――𝛾𝑧 = 4𝛾𝑧 = 810

20

30

angle of pins 𝜃p [°]Fig. 4.15: 星状ピンの𝐽c -𝜃p 特性
𝐵 = 0.4．対称性より，すべての角度は0° ≤ 𝜃p ≤ 30°のうちのいずれかと等価である．

- 83 -

次に，幾何学上，角度による形状の変化がまったくない，面状ピンと比較する形で，𝐽c -𝐵特性をFig. 4.16に示す．なお面状ピンは，星状ピンと同じく𝑥 = 3, 7の2平面に1枚ずつ配置し，厚み1
2 2

normalized critical current density 𝐽c ((3)

3

0.2

√2 𝐽

𝐽→𝐽)

は0.3とした．

0.1

0

―――
―――
―――
―――

𝛾𝑧 = 1𝛾𝑧 = 2𝛾𝑧 = 4𝛾𝑧 = 8―――plane pin

0.1

0.2

0.3

0.4

0.5
𝐵

normalized magnetic field 𝐵 (𝐵 → 𝐵)
c2

Fig. 4.16: 星状ピンと面状ピンの𝐽c -𝐵特性
𝐵 = 0.4．

- 84 -

0.6

Fig. 4.17: 星状ピンと比較するための面状ピン

星状ピンの𝐽c は，おおむね面状ピンのそれよりも0.02ほど低い水準であるが， 𝐵 = 0.35,
𝛾𝑧 = 1 においては明らかに面状ピンよりも高い．この理由としては，マッチング機構によるピーク
効果が考えられる．星状ピンは角度の異なる複数の柱状ピンから構成されているが，そのうちの特定の柱状ピンが磁束線を捉えたとき，磁束線の𝑥座標だけでなく𝑦座標も固定する．これにより磁束線同士の距離も固定され，ピーク効果があらわれる．この一連の作用は純粋な柱状ピンでも強く起こるが，面状ピンでは起こりえない．さらに，ピーク以外の箇所の𝐽c も低いわけではなく，平行な柱状ピンと比較すると𝐵の変化に対して安定しているといえる．
一部，磁束線がピンを跨ぎながら留められることを前提とした配置であるから，異方性が強い
場合には面状ピンよりも優れた箇所が見当たらなかった．しかし，より大きく数値解析を実行し，星状ピンを繰り返しパターンで配置できれば，平行となる柱状ピンが増えることから，違った結果が得られる可能性はある．

- 85 -

第5章

結言

本研究では，3次元TDGL方程式を数値的に解くことで超伝導体内の磁束線の動きを可視化し，超伝導体の条件による𝐽c への影響を調査した．柱状ピンは磁束線と平行であるほどより力を発揮するだろうというポイントから，柱状ピンと磁束線のなす角に注目して数値解析を行った．
しかし，実際には平行でなくても，柱状ピンはピン止めに貢献しうることがわかった．そういった調査から，柱状ピンの良い特性を残しつつも，磁束線の向きに左右されにくい星状ピンを提案した．
また，今回は異方性を持つ超伝導体の数値解析を行った．これによりこれまで再現できなかった性質を持つ超伝導体の再現が可能となった．そこから得られた結果からはとくに，異方性の強さが引き起こす磁束線のつながりの弱さについて理解をすることができた．磁束線のつながりが弱ければ磁束線は柔軟性を増し，ピンに容易に沿うことができるようになる一方で，磁束線の一部にかかるピン力を全体に伝えられないようにもなる．つまり異方性の強さは，磁束線とピンの配置の関係によって，有利にはたらく場合もあり，不利にはたらく場合もある．
よりシミュレータを高速化させるなどして，大きな超伝導体を扱うことができれば，さらに複雑な状況も再現できると考える．

- 86 -

第6章
6.1

付録

GPGPUを用いたプログラミング

今回の数値計算の高速化は，GPGPUプログラミングによるところが大きい．GPGPUを利用することで得られる恩恵のうちでもっとも大きいものは，膨大な数のスレッドを同時に働かせる並
列処理能力である．1本のスレッドの計算速度も十分に速いが，これを束ねることによる圧倒的な計算量や，独特なアーキテクチャで問題を解決する．なお，GPGPUはCPUを使った通常のプログラムコードの中に，GPUで動く関数を挿入することで記述されるが，この関数のことをカーネルと呼ぶ．本節では，GPGPUを使ったプログラミングにおける高速化の工夫を簡単に紹介する．

6.1.1 for文の削減
コンピュータプログラムが時間をとるもっとも大きな原因はfor文である．N回for文を繰り返す処理があるとすると，GPGPUを使えばスレッドが不足しない限りN本のスレッドで1度にそれを終了させることができる．これはGPGPUを使うにあたっての基礎的な発想である．1次元的なfor文
だけでなく，高次元的なfor文で記述された処理(例えば2次元配列を使った処理)も，うまく1次
元的なfor文に展開させることで，1度に終了させることができる．
また，カーネル内でのfor文も少なく抑えたほうがよい．カーネル内でfor文が増えるとコンパイ
ル時の最適化が不十分になったり，場合によってはスレッド間でfor文を繰り返す回数が一致しないために不具合が発生したりするからである．これらを避けるために，繰り返す処理単位でのスレッドの分割や，ループアンローリングを行うとよい．

- 87 -

6.1.2 時間計算量の削減
すべてのfor文を削除することができれば，あらゆるアルゴリズムが要する時間計算量は𝑂(1)に帰着する．しかしそれができないのは，次のようなときに並列計算ができないからである．1つ
めは，1つの変数に計算結果を集約する場合である．2つめは，計算結果をほかの計算で利用する場合である．
2つめに関しては自明であるが，1つめに関しては，たとえば数列の総和を求める場合である．
数列の要素1つに対して1本のスレッドを割り当てたとしても，複数のスレッドが答えを格納する変数に同時にアクセスすることになるため，処理順が不定となるうえ，結局は1つずつ足すことと変わらない．この問題をより早く解決するには，並列処理を利用し，アルゴリズム的に時間計算量を削減すればよい．例に挙げたような数列の総和であれば，次のようなアルゴリズムにより，同一変数への多重アクセスを避けながら時間計算量を𝑂(𝑛)から𝑂 (log 𝑛)に削減できる．数列の要
素1つに対して1本のスレッドを割り当てるとし，簡単化のために数列には2の自然数乗個の要素が含まれているとする．

イ) 変数𝑤を定義し1を代入する．
ロ) 数列の2 × 𝑤の倍数番目の要素を担当する2 × 𝑤の倍数号スレッドは，みずからより𝑤だ
け加算した号数のスレッドが担当する要素の値をみずからが担当する要素の値に足す．
ハ) 𝑤を2倍し，それが数列の要素数未満ならばロ)に戻る．
ニ) 開始時点で初項であった要素に，総和が格納されている．

- 88 -

6.1.3 条件分岐の削減
前述したとおりGPUは条件分岐を含む処理を苦手とする．条件判定自体に時間がかかるという問題以外にも，分岐した結果として処理終了のタイミングがずれると，早いスレッドは遅いス
レッドに合わせてしまうという問題もある．これに関してはブロックに関する知識が必要になるが，ここでは詳しくは述べない．
条件判定を減らすためには，if文の代わりに変数を用いる方法がある．たとえば偶数号スレッ
ドには担当する数値をインクリメントしてほしいが，奇数号スレッドには担当する数値をデクリメン
トしてほしい場合が考えられる．この場合は，整数型変数i をスレッドの号数，担当する数値をとして， + 1 − 2 × (i %2) →とすることでif文を避けられる．変数の型などによっては，むしろ素直にif文を使ったほうが早い場合もあるということに留意されたい．

6.1.4 カーネルやカーネル引数の削減
カーネルを増やすことで高速化するということを述べたが，1本のカーネルが要する時間が同じならばカーネルが少ないほうが早い．また，カーネルへ渡す引数は，一様ではないもののデー
タ総量が少ないほど早い．よって複数の変数を1つの変数に圧縮して渡し，カーネル内で解凍すると処理時間が改善されることがある．例えば引数として渡すときだけ，年月日を3つではなく1つの変数にまとめるなど．

6.1.5 メモリコピー回数の削減
CPUが直接アクセスするメモリとGPUが直接アクセスするメモリは一部を除き異なる．GPUで計算した結果を取り出すためにメモリのコピーが必要になるが，これにはかなりの長い時間を要する．今回のシミュレータはメモリコピーの回数を減らすために，計算条件を変えることなく，また，
計算結果を取り出すことなく1000ステップ以上にわたって計算を継続できるように設計された．
- 89 -

6.2

鎖交体積を求めるために用いた方法

4.2.2項では，柱状ピンと磁束線の鎖交体積からピン力に関する考察を行った．鎖交体積を求めるために用いた方法をここに記す．

6.2.1 柱状ピンと磁束線が 1 本ずつの場合
この場合は，鎖交体積を幾何学的に求めることが容易である． Fig. 6.1に示すとおり，超伝導体が𝑥, 𝑦, 𝑧座標にそれぞれ𝐿𝑥 , 𝐿𝑦 , 𝐿𝑧 の長さをもつ直方体であるとする．

𝑦

𝑟

𝐿𝑦
𝜃p

flux line

𝑧

𝜉

0

𝐿𝑧

Fig. 6.1: 柱状ピンと磁束線が1本ずつの場合の鎖交体積を幾何学的に求めるための模式図
図中桃色の部分は柱状ピンと磁束線の鎖交部分である．

そして柱状ピンと磁束線が，どちらも超伝導体と比べて十分に長い正四角柱であると仮定する．
磁束線は2𝜉 × 2𝜉の断面積を持ち，ピンは2𝑟 × 2𝑟の断面積を持つ．さらに，磁束線は常に𝑧軸- 90 -

に平行であるが，ピンは𝑧軸と𝜃p だけ向きが異なるとする．𝜃p = 0°のとき𝑧軸と平行となり，𝜃p =
90°のとき𝑦軸と平行となる．0° ≤ 𝜃p ≤ 90°とする．ピンの回転の中心は超伝導体の中央にあり，ここに磁束線の中心線も通っているとする．
まずは簡単なアプローチを行う．𝐿𝑧 = ∞とする．磁束線とピンが交わっている箇所は，𝑦-𝑧平
面上では，0° < 𝜃p < 90°のとき平行四辺形を描く．この平行四辺形の面積は4𝑟𝜉/si 𝜃P となる．𝑥軸も考慮すると，𝑟 < 𝜉のとき鎖交体積は8𝑟 2 𝜉/si 𝜃Pとなる．
つぎに𝐿𝑧 が有限である場合を考える．Fig. 6.2に示すように𝑙a , 𝑙b , 𝑙c を定義する．

𝑦
𝑙b
𝑙c

𝑙a

𝐿𝑦
+𝜉
2

𝜃p

𝐿𝑦
2

𝑧

𝐿𝑧
2

𝐿𝑧

Fig. 6.2: 𝐿𝑧 が有限である場合に鎖交体積を幾何学的に求めるための模式図
Fig. 6.1に示したうちの右上の部分に注目している．図中黄色の部分は，鎖交部分のうち，超伝導体からはみ出た部分である．

まず𝑙a について次の2式が言える．
𝑙b = 𝑙a cos𝜃P(6.1)

𝜉 = 𝑙a si 𝜃P(6.2)- 91 -

これらから，次の式が成り立つ．
𝑙b
(6.3)

= 𝜉cos𝜃P /si 𝜃P
= 𝜉/t

𝜃P

一方で
(6.4)

𝑙c = 𝑟/si 𝜃Pである．𝐿𝑧 が有限であれば鎖交部分が超伝導体からはみ出ることがあるが，その条件は次のとおりである．
𝑙b + 𝑙c > 𝐿𝑧 /2(6.5)

𝜉cos𝜃P /si 𝜃P + 𝑟/si 𝜃P > 𝐿𝑧 /2
鎖交部分のうち超伝導体からはみ出た部分の面積を𝑆T とする．
1 𝜉 cos 𝜃P + 𝑟 𝐿z 2
𝑆T = (
− ) t
2
si 𝜃P
2

(6.6)

𝜃P

このとき鎖交体積は2𝑟 {= 2𝑟 {=

4𝑟𝜉
− 2 × 𝑆T }
si 𝜃P

4𝑟𝜉
1 𝜉 cos 𝜃P + 𝑟 𝐿z 2
−2× (
− ) t
si 𝜃P
2
si 𝜃P
2

8𝑟 2 𝜉
𝜉 cos 𝜃P + 𝑟 𝐿z 2
− 2𝑟 (
− ) t
si 𝜃P
si 𝜃P
2

𝜃P }

(6.7)

𝜃P

これよりもさらに𝜃p が小さくなると
𝑙b − 𝑙c > 𝐿𝑧 /2(6.8)

𝜉cos𝜃P /si 𝜃P − 𝑟/si 𝜃P > 𝐿𝑧 /2
この条件下では鎖交体積は4𝑟 2 𝐿𝑧 /cos𝜃P になる．以上から𝐿𝑧 が有限である場合の鎖交体積は
つぎのようになる．

- 92 -

8𝑟 2 𝜉
si 𝜃P
2
8𝑟 𝜉
𝜉 cos 𝜃P + 𝑟 𝐿z 2
− 2𝑟 (
− ) t
si 𝜃P
si 𝜃P
2
2
4𝑟 𝐿𝑧
{
cos𝜃P

𝜉 cos 𝜃P + 𝑟 𝐿z
<
si 𝜃P
2
𝜉 cos 𝜃P − 𝑟 𝐿z 𝜉 cos 𝜃P + 𝑟
<
<
si 𝜃P
2
si 𝜃P
𝐿z 𝜉 cos 𝜃P − 𝑟
<
2
si 𝜃P
𝑟<

𝜃P

(6.9)

6.2.2 柱状ピンが 2 本ある場合
ある𝑦-𝑧平面に，柱状ピンの中心線が2本含まれる場合を考える．この場合は，鎖交体積を幾何学的に求めることが容易ではない．そこで，モンテカルロ法を用いて近似的に求めた．3次元空間を定義し，シミュレーションと同様に磁束線やピンを配置する．乱数を用いて空間内の座標を参照し，その座標に磁束線とピンが両方とも存在していれば，その座標は磁束線とピンの鎖交部分に含まれる．乱数で座標を参照した総数と，そのうちで鎖交部分であった数の割合か
ら，空間全体の体積をもとにして鎖交体積を求める．
ただし柱状ピンが1本の場合と異なり，磁束線の𝑦座標によって鎖交体積が変わる．今回は，磁束線の𝑦座標ごとに鎖交体積をとり，最大となった値を採用することとした．これは，磁束線はできるだけ常伝導領域を通ろうとする性質があるため，最もピンと鎖交できる𝑦座標に磁束線が来るのが自然であるだろうとしたことによる．

- 93 -

参考文献[1] V. L. Ginzburg, L. D. Landau, and Zh. Eksp (1950) Teor. Fiz. 20, 1064
[2] T. Matsushita (2014) Flux Pinning in Superconductor, second ed., Springer[3] 松下照男 (2014) 超伝導応用の基礎, 産業図書株式会社, 東京[4] D. Y. Vodolazov (2013) Phys. Rev. B 88, 014525[5] J. F. Blackburn et. al (2000) Philosophical Magazine 80, 1455
[6] E. S. Otabe and T. Matsushita (1993) Cryogenics 33, 531[7] T. Matsushita, E. S. Otabe, and T. Matsuno (1990) Adv. Cry. Eng. Mater. 36, 263[8] T. Matsushita (1983) J. Appl. Phys. 54, 281[9] T. Matsushita (1982) J. Phys. Soc. Jpn. 51, 2755
[10] D. A. Jacobson (1965) Phys. Rev. 138, 1066
[11] K. E. Osborne (1965) Phil. Mag. 23, 1113
[12] Yu. F. Bychkov, V. G. Vereshchagin, V. R. Karasik, and G. B. Kurganov (1969) Sov. Phys.
JETP 29, 276[13] H. Küpfer, I. Apfelstedt, R. Flükiger, C. Keller, R. Meier-Hirmer, B. Runtsch, A. Turowski,
U. Wiech, and T. Wolf (1969) Cryogenics 29, 268[14] 高橋亮一, 棚町芳弘 (1991) 計算力学とCAOシリーズ差分法, 培風館, 東京[15] 田中敏幸 (2006) 数値計算法基礎, コロナ社, 東京[16] 戸川隼人 (1987) 数値解析とシミュレーション, 共立出版, 東京[17] 新濃清志, 船田哲男 (1991) 数値解析の基礎-理論とPAD・PASCAL・C, 培風館, 東京[18] A. Schmid (1966) Phys. Condens. Mat. 5, 302[19] M. Tinkham (1996) Introduction to Superconductivity, second ed., McGraw-Hill[20] e.g. R. S. Thompson and C.-R. Hu (1971) Phys. Rev. Lett. 27, 1352
[21] N. Kopnin (2001) Theory of Nonequilibrium Superconductivity, Clarendon Press, Oxford- 94 -

[22] R. Kato, Y. Enomoto, and S. Maekawa (1991) Phys. Rev. B, 44, 6916
[23] 松本 圭司, 松本 要 (2011) 平成23年度日本金属学会九州支部学術講演会, B5[24] 一野祐亮, 伊藤慎太郎, 吉田隆 (2015) 第75回応用物理学会秋季学術講演会, 17pA21-15[25] W. E. Lawrence and S. Doniach (1970) in Proceedings of the Twelth Conference on Low
Temperature Physics, Kyoto[26] Knuth Donald (1997) Sorting and Searching, The Art of Computer Programming, 3rd ed.,Addison-Wesley, 409[27] J. G. Bednorz and K. A. Müller (1986) Z. Physik, B 64(1), 189[28] R. Kato, Y. Enomoto, and S. Maekawa (1993) Phys. Rev. B, 47, 2016[29] M. Machida and H. Kaburaki (1993) Phys. Rev. Lett. 71, 3206
[30] Tamegai and Nakajima (2009) arXiv, 0906.0444[31] Tamegai and Nakajima (2010) Physica, C 470, 1103
[32] Tamegai and Nakajima (2009) Phys. Rev. B 80, 012510[33] 米塚 里奈, 濱田 雄成, 上地 和典, 谷村 賢太, 吉原 敬貴, 小田部 荘司, 馬渡 康徳,
松野 哲也 (2018) 平成30年度秋季低温工学・超電導学会, 1P-p11[34] T. Yoshihara, K. Tanimura, E. S. Otabe, M. Kiuchi, Y. Mawatari, and T. Matsuno (2018)
the 65th JSAP Spring Meeting, 17p, B401, 8[35] I. A. Sadovskyy, A. E. Koshelev, C. L. Phillips, D. A. Karpeyev, and A. Glatz (2015)
Journal of Computational Physics, 294, 639[36] I. A. Sadovskyy, A. E. Koshelev, A. Glatz, V. Ortalan, M. W. Rupich, and M. Leroux(2016) Phys. Rev. Appl. 5, 014011- 95 -

謝辞
本研究は，熱心なご指導を続けてくださった小田部荘司教授をはじめとした，多くの方々からのご助力をいただき，実現しました．
小田部先生からは研究に関する知識やアイデアを次々とご教示いただきました．またそれだ
けでなく，就職活動に熱心にご協力くださったり，学生あるいは社会人としてのあるべき姿を説いてくださったり，バングラデシュの大学生との交流の機会を与えてくださったりと，その多岐にわたるご支援は枚挙に暇がありません．心から感謝申し上げます．
次に，松下照男名誉教授からは，この研究の土台となる，現代物理学・超伝導の基礎的な知識をご教授いただきました．心から感謝申し上げます．
また，共同研究者の有明工業高等専門学校の松野哲也先生と産業技術総合研究所の馬渡康徳先生に心より感謝申し上げます．お2人からは，共同研究ゆえの的確さを持ちつつ切り口の異なるアドバイスを多くいただきました．そういったご指摘によって，本研究の中で幅広いアプローチをとることができました．
最後に，同研究室の米塚里奈さんと上地和典とは，同じTDGLを扱うグループとして協力させていただき，日々の研究において何度も助けていただきました．また，直接的に研究に関すること以外にも，日常で交わすコミュニケーション等は，私にとって大きな支えとなりました．厚く御礼申し上げ，感謝いたします．

- 96 -