平成 27 年度
卒業論文
時間依存性 Ginzburg-Landau 方程式を用いた超伝導体内の磁束線評価

久野昌隆(学籍番号：12232066)

九州工業大学

情報工学部

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

平成 28 年 2 月 23 日

教授

目次
第1章

序章 ................................................................................................................................................................. 1

1.1 はじめに ................................................................................................................. 1
1.2 縦磁界効果 .............................................................................................................. 2
1.2.1 磁束カッティングモデル ............................................................................... 3
1.2.2 force-free モデル ........................................................................................... 5
1.3 磁束ピンニング ....................................................................................................... 8
1.4 Time-Dependent Ginzburg-Landau 方程式 ........................................................... 9
1.5 オイラー法 ............................................................................................................ 12
1.6 シンプレクティック法 .......................................................................................... 13
1.7 ガウス・ザイデル法 .............................................................................................. 14
1.8 Finite Difference Time Domain 法....................................................................... 15
1.8.1 差分法........................................................................................................... 15
1.8.2 差分近似の精度............................................................................................. 16
1.8.3 差分近似の注意点 ......................................................................................... 16
1.8.4 FDTD 法 ....................................................................................................... 17
1.9 波動方程式 ............................................................................................................ 17
1.10 拡散方程式 .......................................................................................................... 18
1.11 シュレディンガー方程式 ....................................................................................... 18
1.12 本研究の目的 ........................................................................................................ 19
第2章

計算手法 .................................................................................................................................................... 21

2.1 計算方法 ............................................................................................................... 21
2.1.1 波動方程式 ........................................................................................................................................ 21
2.1.2 拡散方程式 ........................................................................................................................................ 22
2.1.3 シュレディンガー方程式 ............................................................................................................ 22
2.1.4 TDGL 方程式 ..................................................................................................................................... 24
第3章

実行結果 .................................................................................................................................................... 39

3.1 波動方程式 ................................................................................................................................................... 39
3.2 拡散方程式 ................................................................................................................................................... 39
3.3 シュレディンガー方程式....................................................................................................................... 40
3.4 横磁界下における 2 次元 TDGL 方程式 .......................................................................................... 41
3.5 横磁界下における 3 次元 TDGL 方程式 .......................................................................................... 43
3.6 縦磁界下における 3 次元 TDGL 方程式 .......................................................................................... 45
3.7 ピンを添加した横磁界下における TDGL 方程式 ....................................................................... 48
3.8 ピンを添加した縦磁界下における TDGL 方程式 ....................................................................... 49

第4章

考察 .............................................................................................................................................................. 51

4.1 波動方程式 ................................................................................................................................................... 51
4.2 拡散方程式 ................................................................................................................................................... 51
4.3 シュレディンガー方程式....................................................................................................................... 51
4.4 横磁界下における 2 次元 TDGL 方程式 .......................................................................................... 51
4.5 横磁界下における 3 次元 TDGL 方程式 .......................................................................................... 53
4.6 縦磁界下における 3 次元 TDGL 方程式 .......................................................................................... 54
4.7 ピンを添加した場合の 3 次元 TDGL 方程式................................................................................. 55
4.8 パラメータが異なる場合の実行結果............................................................................................... 56
第5章

まとめ ......................................................................................................................................................... 57

謝辞 ................................................................................................................................................................................. 58
参考文献........................................................................................................................................................................ 59
付録 ................................................................................................................................................................................. 60

図目次
図 1.1 超伝導体の電気特性[1]。第一種超伝導体と第二種超伝導体の相違。..................... 2
図 1.2 超伝導体に対し磁界と電流を平行に印加した状態 ................................................. 3
図 1.3 磁束カッティングの模式図[2]。(a)は inter-cutting、(b)は intra-cutting。 ........... 4
図 1.4 磁束構造の歪み[2]。(a)は磁束密度の勾配、(b)は磁束の曲げ、(c)は force-free 電流が流れる時のせん断歪み。矢印方向にトルクが働くと考えられている。 ............... 6
図 1.5 ねじれた磁束線の運動[2]。𝑣1 は force-free トルクによる回転運動成分であり、𝑣2 はこれに誘導された並進運動成分 .............................................................................. 7
図 1.6 超伝導円柱におけるらせん状磁束フロー[2]。(a)は円柱の中心を通る磁束線の動きを示し、(b)は円柱の長さ方向の各部分での磁束線の運動を示す。角度は相対的な磁束線フロー方向の角度を示す。 .............................................................................. 7
図 1.7 オーダーパラメータ𝛹と超伝導体内における局所的な磁束密度𝑏の関係[2] .......... 12
図 1.8 オイラー法の概念図 ............................................................................................. 13
図 2.1 2 次元横磁界の模式図 .......................................................................................... 37
図 2.2 3 次元の模式図..................................................................................................... 37
図 2.3 超伝導体を𝑦軸正方向から見た図 ......................................................................... 38
図 2.4 超伝導体を𝑧軸正方向から見た図 .......................................................................... 38
図 3.1 波動方程式の計算結果。(a)初期状態𝑡 = 0。(b)境界にぶつかり反射したあとの状態
𝑡 = 3 .................................................................................................................... 39
図 3.2 拡散方程式の計算結果。𝑡 = (a)0。(b)10。(c)30 .................................................. 39
図 3.3 シュレディンガー方程式の計算結果。𝑡 = (a)0, (b) 800, (c) 1600, (d) 2000。 ...... 40
図 3.4 横磁界下における 2 次元 TDGL 方程式の計算結果。𝑡 =(a)5, (b)12, (c)65, (d)300。
............................................................................................................................. 42
図 3.5 横磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 =(a)5, (b)12, (c)20, (d)65。
............................................................................................................................. 44
図 3.6 縦磁界下における 3 次元 TDGL 方程式の計算結果。
𝑡 = (a)10, (b) 30, (c)45, (d)60。
............................................................................................................................. 46
図 3.7 縦磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 = (a)150, (b) 250, (c)270,(d)300。 ................................................................................................................ 47
図 3.8

ピンを添加した場合の横磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 =(a)150, (b) 155。 ................................................................................................... 48

図 3.9

ピンを添加した場合の縦磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 =(a)150, (b) 155, (c)160, (d) 165。 .......................................................................... 50

図 4.1 2 次元横磁界におけるエネルギーの時間依存性.................................................... 52

図 4.2 3 次元横磁界におけるエネルギーの時間依存性.................................................... 53
図 4.3 3 次元縦磁界におけるエネルギーの時間依存性.................................................... 55
図 4.3 パラメータが異なる場合の相図 ........................................................................... 56

表目次
表 2.1 2 次元横磁界下における計算条件のパラメータ.................................................... 29
表 2.2 3 次元横磁界下における計算条件のパラメータ.................................................... 32
表 2.3 3 次元縦磁界下における計算条件のパラメータ.................................................... 36

第1章

序論
1.1 はじめに1908 年にオランダの Leiden 大学の Kamerlingh Onnes が世界で初めてヘリウムの液化に成功した。1911 年に液体ヘリウムを用いて水銀の電気抵抗を測定する過程において、4.2 K 以下で電気抵抗が突然ゼロになることを発見した。また、1933 年に Fritz Walther Meissner によって超伝導の完全反磁性(マイスナー効果)が発見された。これらの現象を超伝導現象といい、一定の条件下で超伝導現象を示す物質を超伝導体という。超伝導体は電気抵抗ゼロの性質を持つことから様々な機器への応用が期待された。しかし、当初発見されたものの多くが、ある磁界を超えると超伝導状態を維持することができなくなるということがわかった。超伝導状態を示さなくなる温度、磁界をそれぞれ臨界温度𝑇c、臨界磁界𝐵c と呼ぶ。ある一定の温度や磁界でしか特性を示さないという超伝導状態になる条件はあるが、非常に有用性が高いため、現在も研究が続いている。
1957 年には John Bardeen, Leon Neil Cooper, John Robert Schrieffer らのBCS 理論によって超伝導現象発現の基本的な発現機構が説明された。BCS 理論によると𝑇cは 30 K を超えないと予想されていたが、1986 年に Johannes Georg
Bednorz と Karl Alex Müllar らによって臨界温度が 30 K を超える酸化物超伝導体が発見された。これ以降、世界各国で高温超伝導体の探索が続けられ、その 1 年後には𝑇cが液体窒素の常圧下での温度を超えるYBa2 Cu3 O7−δ (𝛿 ≈ 0.14)が発見された。このような高い𝑇c持つ超伝導体は高温超伝導体と呼ばれる。また、その中でも銅酸化物であるものを銅酸化物超伝導体と呼ぶ。これらの超伝導体
は液体ヘリウムに比べて安価な液体窒素中で超伝導状態になるため、様々な機器への応用の可能性や冷却コストの低減などの点から大きな注目を浴びた。しかし、高温超伝導体は実用化に向けて、加工が難しいことや臨界電流密度を高めることが難しいという課題があり、今日も研究が続けられている。
超伝導材料の性能を表す指標として、電気抵抗ゼロの状態で流せる電流の最大値である臨界電流𝐼c 、その電流密度である臨界電流密度𝐽c が存在する。断面積が同じであれば𝐽c が大きいほど超伝導体の電気容量が大きくなる。したがって、
超伝導体の𝐽c を大きくすることが工学的な応用の際に重要である。
超伝導体は磁界に対する特性の違いから第一種超伝導体と第二種超伝導体に1分類される。第一種超伝導体は超伝導状態になるとマイスナー効果を示すが、外部磁界を加えていくと、一定の磁界を超えるとマイスナー効果が消滅し、同時に超伝導状態も失われる。第二種超伝導体も同様にマイスナー効果を示し、外部磁界を加えると一定の磁界の強さに達したところでマイスナー効果を失う
が、磁束の一部が超伝導体の内部に侵入した状態で超伝導状態を保つことができる。この状態を混合状態といい、磁束が侵入し始める磁界を下部臨界磁界𝐵c1 という。また、混合状態からさらに外部磁界を加えると、一定の磁界を超えて超伝導状態が消滅する。この磁界を上部臨界磁界𝐵c2 という。𝐵c2 は𝐵c1 よりも非常に大きいため工学的な応用への期待があり、超伝導体をエネルギー機器等へ応用することが現実的なものとなった。

図 1.1:超伝導体の電気特性[1]。第一種超伝導体と第二種超伝導体の相違。

1.2 縦磁界効果
縦磁界とは、超伝導体に流している電流に対し、平行に磁界を印加した状態のことを言う。一般的に、超伝導体に流れる電流によって生じる自己磁界を含めて磁界は垂直方向となるが、図 1.2 に示すように、磁界と電流が平行な場合、垂直な場合とは異なる様々な現象が観測される。以下にその時の例を示す。
1. 電流によって磁界と同方向の磁化が正となる。これを常磁性効果と呼ぶ。
2. 外部磁界を増加させると交流電流による損失が減少する。
3. 縦磁界の場合は、
𝑱×𝑩=0
2

(1.1)

のように磁束線に対して Lorentz 力が働かないため臨界電流密度が横磁界の場合に比べ大幅に増加する[3, 4]。
4. 磁束線の運動と電磁現象を結びつける Josephson の関係式( 𝑬 = 𝑩 × 𝒗 )が成立しない[5]。
5. 臨界電流密度を超えた抵抗状態において、らせん状に縦方向に負の電界を有する領域が存在する。
これらを総称して縦磁界効果と呼ぶ。3 に関しては、超伝導体の形状によって臨界電流密度の増加量が大きく変化する。この特徴を用いて、超伝導体ケーブルの設計がなされている[6]。
一方で、4 について、実際にどのような磁束線の運動が起きているかについて
は、磁束線同士が交錯する際に、磁束線が切れるという磁束カッティングモデ
ルや、電流と磁界が平行なため局所的に Lorentz 力が存在しないことを基にした force-free モデルなどが提唱されている｡

𝐽𝑩
図 1.2:超伝導体に対し磁界と電流を平行に印加した状態

1.2.1 磁束カッティングモデル
Cave らの実験[5]によると、縦磁界下で超伝導円柱に直流電流に微小交流電流を重畳して流した結果、誘導電界はほとんど縦方向を向いている一方、磁束の縦成分はほとんど変化せず、方位角成分のみ変化していた。これより、一定角度を保ったまま超伝導体内へ磁束線が進入するという磁束線運動が起きておらず、Josephson の関係式が成り立たないことが明らかとなった。Cave らは磁束3線の縦成分と方位角成分とが交差しあって交流電流による方位角成分のみが出入りするという、磁束カッティングモデルを提案した。
磁束カッティングモデルには大きく分けることができ、異なった角度をなす
磁束線同士が交差するもの(inter-cutting)と、異なる磁束線同士で編成しなおす
もの(intra-cutting)の 2 つに分類できる。Cave らのモデル[3]は inter-cutting の例の 1 つであり、Clem[7]や Brandt[8]の分断と再結合のモデル(intersection andcross-joining)は intra-cutting の 1 つである。図 1.3 の(a)に Cave らが提唱した
直接的なカッティングモデルの、(b)に Clem や Brandt が提唱した分断と再結合のモデルの模式図を示す。図 1.3(a)では、外部磁界𝐵𝑧 による磁束線と輸送電流から生じる自己磁界𝐵𝜙 による磁束線の 2 種類が存在する。この𝐵𝜙 による磁束線が中心に向かって移動していくが、その際に𝐵𝑧 による磁束線を通り抜けていくのである。次に図 1.3(b)についてだが、(1)では手前側と奥に互いに異なった磁束線があり、両者の間で最も近寄った部分で磁束線が切断され、(2)のように継ぎ換えが行われる。これがまっすぐに戻って(3)のようになる。(1)と(3)を比較してみると、横成分のみやり取りを行い、縦成分に関しては不変である。この磁束線の組み換えが連続して行うことによって、横成分のみ進入することが実現
される。

(a)

(b)
(1)

(2)

(3)

図 1.3:磁束カッティングの模式図[2]。(a)は inter-cutting、(b)は intra-cutting。
の例

4

1.2.2 force-free モデル
一般的に電流が流れる磁束構造を考えてみる。磁束密度を𝑩 = 𝐵𝒊𝐵と書く。ここで、𝒊𝐵 は𝑩の単位ベクトルである。すると電流密度は
11
1(𝒊𝐵 × 𝛻 )𝐵 + 𝐵 rot𝒊𝐵
𝑱 = rot𝑩 = −
𝜇0
𝜇0
𝜇0

(1.2)

(1.3)

と表される。この電流は 3 つの成分から、すなわち、磁束密度の大きさ𝐵の勾配による成分、磁束の曲がりによる成分および force-free 電流成分である。第 1の成分は磁気圧に、第 2 の成分は線長力に寄与する。すなわち、これらの成分
は共に Lorentz 力に寄与する。磁気圧に関する成分は(1.3)式の第 1 項にのみ含
まれ、一方、force-free 電流成分は第 2 項にのみ含まれる。そして線張力に関係する成分は両方の項から生じる。
以上のように、一般に電流が流れた時には、磁束構造は歪みをもつ。図 1.4に各種の基本的な歪みを示す。(a)は磁束線の勾配、(b)は磁束の曲げ、(c)はforce-free 電流が流れるときのせん断歪みである。それぞれにおける矢印方向にトルクが働くと考えられる。
ここで、実際の実験について考える。縦磁界下における臨界状態は force-free
電流の駆動トルクとピン力が釣り合っている状態である。図 1.5 のような磁束線構造を考えると、臨界角δ𝜃c を超えた磁束線は、それを緩和させようと𝛿𝜃を減少させる𝑣1 方向に動く。さらに、
∂𝜃
𝜕𝜃
= 𝑣2
(1.4)
𝜕𝑡
𝜕𝑦
という定常状態の条件を満たすためには𝑣1 の回転運動の他に、𝑣2 の並進運動も必要となる。このとき、円柱状超伝導体の場合に期待される磁束線の動きを図1.6 に示す。(a)は円柱の中心を通る磁束線の運動の方向を示したものであり、(b)は円柱の長さ方向の各部分での磁束線の運動の方向を示したものである。こうした磁束線運動は一様な円柱内の直線的な磁束線フローを円柱の軸の周りにねじった形となっている。

5

(a)

(b)

(c)

図 1.4:磁束構造の歪み[2]。(a)は磁束密度の勾配、(b)は磁束の曲げ、(c)はforce-free 電流が流れる時のせん断歪み。矢印方向にトルクが働くと考えられている。

6

図 1.5:ねじれた磁束線の運動[2]。𝑣1 は force-free トルクによる回転運動成分であり、𝑣2 はこれに誘導された並進運動成分
(b)

(a)

図 1.6:超伝導円柱におけるらせん状磁束フロー[2]。(a)は円柱の中心を通る
磁束線の動きを示し、(b)は円柱の長さ方向の各部分での磁束線の運動を示す。角度は相対的な磁束線フロー方向の角度を示す。

7

1.3 磁束ピンニング
超伝導体は臨界温度𝑇cと臨界磁界𝐵c によって、超伝導状態を維持できる環境が制限されている。また、超伝導体を実際の機器で応用する際においては、数ある超伝導状態の特性の中でも、「電気抵抗ゼロ」という特性が最も魅力的であると考えられる。しかし、電気抵抗ゼロのままに無制限に電流を流せるわけでなく、一定以上の電流を流すと電気抵抗が発生し、電気抵抗ゼロとはならなくなる。このときの電気抵抗ゼロの状態で超伝導体に流すことができる最大の電流量を、臨界電流𝐼c と呼び、そのときの電流密度を、臨界電流密度𝐽c と呼ぶ。𝐼c が大きいと、単純にその超伝導体で扱える電流量が増えるので、𝐼c を大きくすることは重要になってくる。同様に、同じサイズの線材でも電流密度𝐽c が大きくなると𝐼c も大きくなるので、𝐽c を大きくすることも重要であり、超伝導体自体のポテ
ンシャルを評価するうえで、𝐽c は重要なパラメータの一つとなる。
この𝐽c の値を決定するものは、基本的に磁束ピンニングという機構である。磁界中において超伝導体に電流を流すと、超伝導体内部に侵入している磁束にLorentz 力が生じる。この Lorentz 力によって磁束線が運動してしまうと、誘導起電力が生じて電気抵抗が生じる。よって、電気抵抗をゼロにするには誘導起電力を生じさせないように、磁束線の運動を止める必要がある。この磁束線の運動を止めようとする作用を、磁束ピンニングと呼ぶ。この磁束ピンニング
は、超伝導体内部の転移、常伝導析出物、空隙、結晶粒界面などといった様々な欠陥により引き起こされる。こうした欠陥をピンニング・センターと呼ぶ。
この作用により、Lorentz 力がある臨界値を超えるまで磁束線の動きを止めることで、電気抵抗ゼロで臨界点（𝐼c ）まで超伝導体に電流を流すことができる。このピンニング・センターが単位体積あたりに及ぼす力を𝐹p とすると、臨界電流密
度𝐽c は

𝐽c =

𝐹p
𝐵

(1.5)

と表される。これは単位体積あたりに及ぼす Lorentz 力にあたる𝐽c 𝐵と、単位体積当たりのピン力𝐹p が釣り合っていることを示す。以上のことから、材料によって先天的に決定される臨界温度𝑇cや臨界磁界𝐵c と違い、臨界電流密度𝐽c は、後天的に上昇させることが可能なパラメータであるため、工学的には𝑇cや𝐵c より、𝐽cの方がより重要なパラメータといえる。つまり、線材加工技術の向上により、𝐹p を大きくすることで、より大きな𝐽c を得られる。

8

1.4 Time-Dependent Ginzburg-Landau 方程式
Ginzburg-Landau 方程式(以下、GL 方程式)とは、1950 年に第二種超伝導体に注目し、熱力学的特性及び磁気的性質をうまく記述する方程式として V. L.
Ginzburg と L. D. Landau によって提唱された[10]。
臨界温度周辺(𝑇 ≈ 𝑇c )では超伝導状態が弱いと考えられるため、磁場の無い場
合の自由エネルギーは𝐸s = 𝐸n + 𝛼 |𝛹 |2 +

𝛽
|𝛹 |4
2

(1.6)

となり、ここでの𝐹s および𝐹n はそれぞれ超伝導状態、常伝導状態における自由エ
ネルギーであり、𝛹 は超伝導状態の秩序状態を示すオーダーパラメータ、
𝛼および𝛽はそれぞれ冪級数展開した際の 1 次と 2 次の係数である。𝑇 < 𝑇c では
𝛼 < 0、 𝛽 > 0である。今回は臨界点近傍であるため、超伝導電子密度となる|𝛹 |2は十分小さい値となると考えられるため、|𝛹 |6 以降の項は無視するものとする。
また、磁場がある場合を考えると、𝛹の空間的な変化を考慮に入れて、量子力学
における運動エネルギーと同様な
𝐸K =1|(−iℏ∇ − 2𝑒𝑨)𝛹 |2
4𝑚

(1.7)

を加える必要がある。ここで𝑚と𝑒 はそれぞれ電子 1 つの質量と電荷量である。
さらに反磁性効果による磁気エネルギーを考えると、
𝐸M =

𝐵2
2𝜇0

(1.8)

を加算する必要があり、ここで𝐻は外部磁界である。
(1.6)～(1.8)式より、自由エネルギーは𝐸 = 𝐸s − 𝐸n =1𝛽𝐵2|(−iℏ∇ − 2𝑒𝑨)𝛹 |2 + 𝛼 |𝛹 |2 + |𝛹 |4 +
4𝑚
2
2𝜇0

(1.9)

となる。ここで𝛹 ∗ と𝑨について変分法を適用すると、
δ𝐸
=0
δ𝛹 ∗

(1.10)

δ𝐸=0
δ𝑨(1.11)

1(−iℏ∇ − 2𝑒𝑨)2 𝛹 + 𝛼𝛹 + 𝛽|𝛹 |2 𝛹 = 04𝑚(1.12)iℏ𝑒 ∗2𝑒 2(𝛹 ∇𝛹 − 𝛹∇𝛹 ∗ ) −|𝛹 |2 𝑨2𝑚𝑚

(1.13)

となり、

𝑱=−

9

を得る。

さらに、時間発展する場合を考える。この時、(1.10)式及び(1.11)式は
δ𝐸
δ𝛹
= −𝛾
∗
δ𝛹
δ𝑡

(1.14)

δ𝐸
δ𝑨
= −𝜈
δ𝑨
δ𝑡

(1.15)となり、γ 及び νはそれぞれ𝛹と𝑨の時定数である。ここでゲージ変換を考えると、ℏ∂𝛹∂→ (ℏ + 2i𝑒𝑉) 𝛹∂𝑡∂𝑡(1.16)∂𝑨 𝜕𝑨(1.17)→+ ∇𝑉∂𝑡𝜕𝑡
であり、𝑉はスカラーポテンシャルである。これを(1.14)式及び(1.15)式に適用す
ると、
γ(

∂𝛹
1
(−iℏ∇ + 2𝑒𝑨)2 𝛹 + 𝛼𝛹 + 𝛽 |𝛹 |2 𝛹 = 0 (1.18)+ 2i𝑒𝑉𝛹) +
∂t
4𝑚

∂𝑨
iℏ𝑒 ∗
2𝑒 2
∗
(𝛹 ∇𝛹 − 𝛹∇𝛹 ) −|𝛹 |2 𝑨 = 0 (1.19)
𝜈 ( + ∇𝑉) + ∇ × ∇ × 𝑨 −
∂𝑡
2𝑚
𝑚
となる。
ただし、(1.18)式及び(1.19)式をそのまま解くのは困難である。そこで、今回は 2 つの簡易化を行うものとする。1 つ目は非常に細い超伝導体を仮定すること
で、超伝導体全体に外部磁界が進入することを仮定する。これにより、𝑨(𝑨の具体的の標式は第 2 章において議論をする。)が外部磁界にのみ依存する変数とな
るため、(1.19)式が一定となる。2 つ目は定数が多いため、規格化を行う。コヒーレンス長𝜉と磁場侵入長𝜆はℏ𝜉=(1.20)√4𝑚|𝛼 |
𝜆=−√2𝑒𝜇0 𝐵c√𝑚|𝛼 |(1.21)と表すことができるため、

−

1
𝑥→𝑥
𝜉

(1.22)|𝛼 |𝑡→𝑡
𝛾(1.23)

2𝑒𝛾
𝑉→𝑉
|𝛼 |

(1.24)

10

𝜆

𝑨→𝑨

(1.25)𝛽 2( ) 𝛹→𝛹
|𝛼 |(1.26)√2𝜇0 𝐵c1のような規格化を行うと、(1.18)式は∂𝛹+ i𝑉𝛹 + (−i∇ − 𝑨)2 𝛹 − 𝛹 + |𝛹 |2 𝛹 = 0∂𝑡という式になる。また、対破壊電流密度𝐽c は(1.27)

32 2 𝐵c𝐽c = ( )3 𝜆(1.28)

を用いると、(1.27)式より規格化された電流密度𝐽は
𝟑

2 𝟐 1
( )
𝐽→𝐽
3 √2 𝐽c(1.29)となる。よって、2/3√3 ≅ 0.385を超えるような電流密度を考えると、対破壊電流を超えた電流を流していることになるため、注意する必要がある。
TDGL 方程式を解く際に、𝛹と𝑉について解くが、方程式が(1.27)式しかない
ため、電流の発散についての式∇∙𝑱=0(1.30)
を第 2 式として解く。電流には超伝導電流𝐽s と常伝導電流𝐽n の 2 つがあり、それぞれ規格化した結果、
𝑱s =1i(𝛹 ∗∇𝛹 − 𝛹∇𝛹 ∗ ) − |𝛹 |2 𝑨
2𝑱n = 𝜎∇𝑉(1.31)(1.32)と表される。ここで、𝛹 ∗ は𝛹の共役複素数、𝜎は常伝導体における導電率である。
これより、(1.30)式は𝑱 = 𝑱s + 𝑱n(1.33)
を用いると、
1(1.34)i(𝛹 ∗ ∇2 𝛹 − 𝛹∇2 𝛹 ∗ ) − ∇ ∙ (|𝛹 |2 𝑨)
2となる。図 1.6 に常伝導核周辺における𝛹と超伝導体内における局所的な磁束密度𝑏の関係を示す。図において左端や右端のように、|𝛹 |2 が高いところでは超伝導状態が強く、その部分では磁束線が進入しにくい。逆に図の中央部分のような、|𝛹 |2 が低いところでは超伝導状態が弱くなって常伝導となっているため、磁𝜎∇2 𝑉 =

束線が侵入しやすい。
11

図 1.7:オーダーパラメータ𝛹と超伝導体内における局所的な磁束密度𝑏の関係[2]1.5 オイラー法
オイラー法は常微分方程式の初期値問題を解く最も基本的な解法である。刻み幅と呼ばれる量ℎを仮に置き、独立変数𝑡を
𝑡𝑛 = 𝑎 + 𝑛ℎ (𝑛 = 0, 1, 2 … )(1.35)のように刻み幅ごとにずらしていったときに、その都度未知関数の値𝑥 𝑖 (𝑡𝑛 )は刻
み幅ℎが十分小さいとき(1.36)𝑥 𝑖 (𝑡𝑛 ) ≈ 𝑥𝑛𝑖というように近似できる。そのとき
𝑓𝑛𝑖 ≡ 𝑓 𝑖 (𝑥𝑛1 , … , 𝑥𝑛𝑚 , 𝑡𝑛 )(1.37)(1.38)と次々定めていくことで𝑥 (𝑡)を求める。𝑚 = 1のときの図を以下の図 1.7 に示す。
なお、この方法は初期条件𝑥0𝑖 が与えられている場合に有効である。
𝑖
𝑥𝑛+1≈ 𝑥 𝑖 + ℎ𝑓𝑛𝑖12

図 1.8:オイラー法の概念図

1.6 シンプレクティック法
シンプレクティック法はハミルトン系を表す常微分方程式を数値積分する方法のうち、計算法そのものが正準変換の性質をもつものである。正準変換とは元の変数を新たなハミルトンの運動方程式を満たす新しい変数に写すという変数変換のことである。平面上の図形を別の平面上に写したとき、面積を保ちつ
つ、形ができるだけ変化しないように精度を高めていく。また、係数パラメー
タが満たす条件式から解法の精度が決まる次数条件式である。
簡単な例として実数値関数𝑞(𝑡)、𝑝(𝑡)に関する
d𝑞
d𝑝(1.39)
= 𝑝,
= −𝑞
d𝑡
d𝑡
のような単振動の方程式を考える。この式は関数𝐻(𝑞, 𝑝) = (1/2 )(𝑝2 + 𝑞2 )をハ
ミルトニアンとするハミルトン系である。この方程式にオイラー法を適用すると、
𝑞= 𝑞n + ℎ𝑝n{ n+1
𝑝n+1 = 𝑝n − ℎ𝑞n(1.40)の差分方程式が得られるため、これをベクトルと行列を用いて表すと、
𝑞n+1
1 ℎ 𝑞n[𝑝 ] = [][ ](1.41)
n+1−ℎ 1 𝑝nのようになる。右辺の行列式は1 + ℎ2 であることから(1.41)の式から(𝑞n , 𝑝n )平面の図形を(𝑞n+1 , 𝑝n+1 )平面上に移動すると、面積が必ず(1 + ℎ2 )倍に拡大される。
また、後退オイラー法を適用すると、

13

𝑞= 𝑞n + ℎ𝑝n+1{ n+1
𝑝n+1 = 𝑝n − ℎ𝑞n+1(1.42)となるので、同様に行列で表すと、
𝑞n+1
1 −ℎ −1 𝑞n(1.43)[𝑝 ] = [] [𝑝 ]n+1
nℎ 1が得られる。(1.43)式によって(𝑞n , 𝑝n )平面の図形を(𝑞n+1 , 𝑝n+1 )平面上に移動すると、面積が必ず(1/(1 + h2 ))倍に縮小される。
シンプレクティック法はオイラー法の(1.40)式の第 2 式の𝑞n だけを𝑞n+1 に置き換えて、𝑞= 𝑞n + ℎ𝑝n{ n+1
𝑝n+1 = 𝑝n − ℎ𝑞n+1(1.44)のようにする。第 1 式を第 2 式に代入すると
𝑞n+1
𝑞n
1
ℎ
[𝑝 ] = [] [𝑝 ]
2n+1
n−ℎ 1 − ℎ(1.45)のように書き直される。この場合は、行列式の値が 1 となることから(1.45)式は面積を不変に保つ写像である。この写像をシンプレクティック写像という。このように、ハミルトン系に適用した際の計算式がシンプレクティック写像となるような数値解析法のことをシンプレクティック数値解法という。

1.7 ガウス・ザイデル法
大きな行列を扱う場合、多くの値を記憶しておく必要があり、大きなメモリが必要となる。その問題の対策としてメモリを節約できる方法がガウス・ザイ
デル法である。
3 元の連立一次方程式
𝑎11
(𝑎21
𝑎31

𝑎12
𝑎22
𝑎32

𝑎13 𝑥1𝑏1𝑎23 ) (𝑥2 ) = (𝑏2 )𝑎33 𝑥3𝑏3(1.46)(𝑘)を解くことを考える。𝑘回目の反復で得られる𝑥1 の値を𝑥1 と表せる。
(𝑘)(𝑘)𝑏1 − 𝑎12 𝑥2 − 𝑎13 𝑥3
(𝑘+1)
𝑥1
=
𝑎11
(𝑘+1)
(𝑘)
𝑏2 − 𝑎21 𝑥1− 𝑎23 𝑥3
(𝑘+1)
𝑥2
=
𝑎22
(𝑘+1)
(𝑘+1)
𝑏3 − 𝑎31 𝑥1− 𝑎32 𝑥2
(𝑘+1)
𝑥3
=
𝑎33

14

(1.47)

(𝑘+1)という反復を繰り返していく。(1.47)式の 2 番目の式で𝑥1が使われており、
(𝑘+1)
新たに𝑥𝑖
を順次求めて使用している。したがって、使用メモリが少なく済む。
(𝑘+1)
即ち、𝑥𝑖
を計算するときに
(𝑘+1)

(𝑘+1)
𝑥𝑖
=

{𝑏𝑖 − (𝑎i1 𝑥1(𝑘+1)+ 𝑎i2 𝑥2

(𝑘+1)

+ ⋯ + 𝑎ii−1 𝑥𝑖−1(𝑘)+ ⋯ + 𝑎in 𝑥n )}𝑎ii(1.48)とする。実際の計算では
(𝑘)

(𝑘+1)
𝑥1
=

(𝑘)

𝑎11
(𝑘+1)

(𝑘+1)
𝑥2
=

𝑏2 − (𝑎21 𝑥1

(𝑘)

(𝑘)

+ 𝑎23 𝑥3 + ⋯ + 𝑎2n 𝑥n )
𝑎22

(𝑘+1)

𝑏3 − (𝑎31 𝑥1(𝑘+1)+ 𝑎32 𝑥2(𝑘)

+ ⋯ + 𝑎3n 𝑥n )(1.49)

𝑎33
…

(𝑘+1)
𝑥3
=

(𝑘)

𝑏1 − (𝑎12 𝑥2 + 𝑎13 𝑥3 + ⋯ + 𝑎1n 𝑥n )

(𝑘+1)

(𝑘+1)
𝑥n
=

𝑏n − (𝑎n1 𝑥1(𝑘+1)+ 𝑎n2 𝑥2

(𝑘+1)

+ ⋯ + 𝑎nn−1 𝑥n−1 )𝑎nnと計算することができる。これを繰り返して連立方程式の解を求める方法がガ
ウス・ザイデル法である。

1.8 Finite Difference Time Domain 法
Finite Difference Time Domain (FDTD 法)は差分法の方式の一つである。まずは差分法との基礎を説明することで関連をつかんでいく。

1.8.1

差分法

差分法とは微分可能な関数𝑓 (𝑥 )を差分式で代用する方法である。ある点𝑥において点𝑥からΔ𝑥はなれた点の関数𝑓の傾きを考えると
𝜕𝑓 (𝑥 )𝑓 (𝑥 + Δ𝑥 ) − 𝑓 (𝑥 )}={(1.50)
𝜕𝑥Δ𝑥Δ𝑥→0のようになる。これを差分式とよぶ。よって(1.50)式の左辺の微分を右辺の差分15

式に近似していることを示している。このときΔ𝑥が小さいほど近似が正確にできるようになる。また、Δ𝑥 → 0の極限で左辺と右辺が等しくなる。差分近似を行う際は安定条件や誤差、境界条件などを注意しておく必要がある。それについては後述する。

1.8.2 差分近似の精度
差分近似が元の微分とどの程度の差があるかを確認するための方法としてテーラー展開を使う。
FDTD 法による式の計算においては時間微分、空間微分を中心差分で置き換えられている。そこで、差分式の変数をテーラー展開する。
𝑓 (𝑥 ± 𝛥𝑥 ) = ±𝜕𝑓
1 𝜕2𝑓 2 1 𝜕3 𝑓 3
Δ𝑥 +
Δ𝑥 ±
Δ𝑥 + ⋯𝜕𝑥
2! 𝜕𝑥 23! 𝜕𝑥 3(1.51)

𝑓(𝑥 + 𝛥𝑥 ) − 𝑓 (𝑥 − Δ𝑥 ) 𝜕𝑓 1 𝜕 3 𝑓 2(1.52)=+
Δ𝑥 + ⋯
2𝛥𝑥
𝜕𝑥 3! 𝜕𝑥 3
上記の式にみられるようにΔ𝑥 2 以下の項は打ち切りによる誤差が入る。この場
合、中心差分における誤差のオーダーは𝑂(Δ𝑥 2 )であり、近似の次数は 1 (𝑂(Δ𝑥 2 )における 2 乗の 2 から 1 を引く)である、という。

1.8.3 差分近似の注意点
微分方程式の解には三通りある。これは時間微分に対する事項である。すな
わち、
1．時間を追って、逐次計算できる解を求めるもの。
2．静的な力が加わって変形しているが、時間的に変化しないもの。
3．一定の周波数で変数が変化しているもの。
の三つである。1 を初期値問題、2 を静的問題、3 を定常問題とよぶ。FDTD 法に関連する問題は初期値問題である。初期値問題では、計算の途中で値が発散したり、異常に振動したりすることがある。これを不安定現象とよび、その原因は
1．定式化に問題がある場合。
2．時間離散間隔と空間離散間隔の比が適当ではない場合。
の二通りある。
解析対象と他の物理系が接する場所を境界と呼び、境界上では、変数の満たすべき条件がありその条件を境界条件という。この境界条件の設定によって解析結果に影響が出る。
16

1.8.4 FDTD 法
FDTD 法では差分式の時間及び空間微分を中心差分で近似する。例として
𝜕𝑓
𝜕𝑢
=𝑐
𝜕𝑡
𝜕𝑥
𝜕𝑢 1 𝜕𝑢=𝜕𝑡 𝜌 𝜕𝑥

(1.53)(1.54)という式を与える。これらの式を近似すると次のような差分式になる。
1
Δ𝑡 𝑛−1
1
1{𝑢 2 (𝑖 + ) − 𝑢𝑛−2 (𝑖 − )}
Δ𝑥
2
2
1
1
1 𝛥𝑡 𝑛𝑢𝑛+1 (𝑖 + 1) = 𝑢𝑛−2 (𝑖 + ) +{𝑓 (𝑖 + 1) − 𝑓 𝑛 (𝑖)}2𝜌 𝛥𝑥𝑓 𝑛 (𝑖 ) = 𝑓 𝑛−1 (𝑖 ) + 𝑐(1.55)(1.56)

ここで𝑛は時間ステップ、𝛥𝑡は離散間隔である。また、𝑖は空間離散点の番号、Δ𝑥は離散間隔である。このように FDTD 法では(1.53), (1.54)式から式の微分を
中心差分で近似し計算する。従って、2 種類の変数を交互に計算していく。

1.9 波動方程式
波を記述する 2 階偏微分方程式のことを波動方程式という。式の形は波の形を表す関数を𝑓(𝑥)としてこの波が一定の速さ𝑣で𝑥方向に波の形を変えずに動いているとすると𝑓(𝑥 − 𝑣𝑡)とおける。すると 1 次元の波動方程式は
∂2
𝜕2
2
𝑓
=
𝑣
𝑓
∂𝑡 2
𝜕𝑥 2
と表せる。2 次元、3 次元も同様に
∂2
𝜕2
𝜕2
2
𝑓
=
𝑣
(
+
)𝑓
∂𝑡 2
𝜕𝑥 2 𝜕𝑦 2
∂2
𝜕2
𝜕2
𝜕2
2
𝑓
=
𝑣
(
+
+
)𝑓
∂𝑡 2
𝜕𝑥 2 𝜕𝑦 2 𝜕𝑧 2

(1.57)(1.58)(1.59)(1.60)と考えることができる。
(1.60)式の右辺はラプラシアン∇2 =△を使うと∂2
𝑓 = 𝑣2 △ 𝑓∂𝑡 2のように表すこともできる。
17(1.61)

1.10 拡散方程式
拡散が生じている物質や物理量の密度の揺らぎを記述する方程式のことを拡散方程式という。
1 次元空間を𝛥𝑥の間隔で離散化し、各格子点上にはある物理量𝐶(𝑥, 𝑡)が割り当
てられているものとする。ここで時刻𝑡から𝑡 + Δ𝑡の間に各格子点上の物理量が隣の格子に等方的に移動することを考えると物理量𝐶(𝑥, 𝑡 + Δ𝑡)期待値は
1
1(1.62)𝐶 (𝑥, 𝑡 + Δ𝑡) = 𝐶 (𝑥 + Δ𝑥, 𝑡) + 𝐶(𝑥 − Δ𝑥, 𝑡)22と表現できる。ここで𝛥𝑥、𝛥𝑡が小さいものとして、(1.62)式の各項をテイラー展
開し、整理すると 1 次元の拡散方程式
∂𝐶
𝜕2𝐶
=𝜅 2
∂𝑡
𝜕𝑥

(1.63)(Δ𝑥 )2(1.64)(> 0)2Δ𝑡が得られる。1 次元格子を正方格子や立方格子に変えて同様に考えると、それぞ𝜅=

れ 2 次元、3 次元の拡散方程式になる。
一般的に拡散方程式は
𝜕𝐶
=𝜅△𝐶
𝜕𝑡

(1.65)の形で表される。

1.11 シュレディンガー方程式
物質と波動をつなぐ式はプランクの関係式からアインシュタインが光量子のエネルギーの式として出した𝐸 = ℎ𝜈
(1.66)と、ド・ブロイの関係式がある。この二式は光や物質で一般に成立し、実験的にも確認されている。
振動数𝜈で波長λをもち、𝑥軸正方向へと伝播する波は𝑥

𝜓λ = 𝑒 2𝜋i(𝜆−𝜈𝑡)(1.67)という式で表すことができる。
(1.67)式は平面波であり、空間の端から端まで同じ振幅で振動している波である。これから作る方程式は線型同次方程式であることを要求する。線型同次であれば、解の重ね合わせができるため、いろんなλに対して𝜓λ を求めると、その18

重ね合わせでさらに多くの解を作ることができる。これを重ね合わせの原理という。この重ね合わせの原理が量子力学でも満たされているならば、複雑な波も簡単な平面波の重ね合わせで表現できる。
(1.67)式においては
𝑝=

ℎ
ℎ 𝜕
𝜕
→ −i
= −iℏ
𝜆
2𝜋 𝜕𝑥
𝜕𝑥

(1.68)

ℎ 𝜕
𝜕
= iℏ
2𝜋 𝜕𝑡
𝜕𝑡
という置き換えをしてもよい。なぜなら(1.69)

𝐸 = ℎ𝜈 → i−iℏ

𝜕 2𝜋i(𝑥−𝜈𝑡)
2𝜋i 2𝜋i(𝑥−𝜈𝑡) ℎ 2𝜋i(𝑥−𝜈𝑡)𝜆𝜆𝜆𝑒= −iℏ ×
𝑒
= 𝑒
𝜕𝑥
𝜆
𝜆

𝑥
𝑥
𝜕 2𝜋i(𝑥−𝜈𝑡)𝜆𝑒= −iℏ × (−2𝜋𝑖𝜈)𝑒 2𝜋i(𝜆−𝜈𝑡) = ℎ𝜈𝑒 2𝜋i(𝜆−𝜈𝑡)𝜕𝑡となるからである。このように

iℏ(1.70)(1.71)(1.72)(演算子) × (関数) = (値) × (関数)となるような関数をその演算子に対する固有関数と呼び、右辺に出てくる値を固有値という。そして(1.72)式のような形の方程式を固有値方程式と呼ぶ。
古典力学においては、エネルギーはハミルトニアン𝐻(𝑝, 𝑥 )として、運動量や∂座標の関数として表された。量子力学におけるエネルギー𝐸 = iℏ ∂tも、同様に運動量や座標と関係付けられるはずである。その関係を波動方程式の形で表したものがシュレディンガー方程式である。
1

非相対論的な 1 次元古典粒子の場合、𝐸 = 𝐻 = 2𝑚 𝑝2 + 𝑉(𝑥)であるから、そのような粒子を表す波は
𝜕ℏ2 𝜕 2iℏ 𝜓 (𝑥, 𝑡) = (−+ 𝑉 (𝑥 )) 𝜓(𝑥, 𝑡)
𝜕𝑡
2𝑚 𝜕𝑥 2

(1.73)のような方程式を満たすであろうと考えることができる。これがシュレディン
ガー方程式である。この𝜓(𝑥, 𝑡)は複素数に値をとる関数で、波動関数と呼ばれる。

1.12 本研究の目的
現在、超伝導と常伝導の相転移を現象的にとらえるためにオーダーパラメー
タについて議論している Ginzburg-Landau 方程式(GL 方程式)は第二種超伝導体の特性を説明するための理論として確立している。それに対して本稿で議論し19ている Time-Dependent Ginzburg-Landau 方程式(TDGL 方程式)は GL 方程式を時間発展させたものである。オーダーパラメータ𝛹の絶対値の二乗と超伝導電子密度𝑛s は比例関係にあるため[2]、𝛹の大小を把握することで磁束線が侵入している様子を観測することができる。したがって、超伝導体に起こる変化を可視化することが可能である。既存の研究としては、電流に対して磁界を垂直に与える横磁界の場合が多く、横磁界中の超伝導体の電圧-電流特性を 2 次元 TDGL 方
程式の計算により求めたものなどがある[16]。
それに対して本稿では電流に磁界を平行に与える縦磁界の場合においてシミュレーションを行う。縦磁界下においては横磁界下に比べて、ピンニングの強さに依存はするが、臨界電流密度が大幅に増加することが分かっている[2]。したがって、臨界電流密度特性が優れている縦磁界下において調べることを目的とする。また、TDGL 方程式の数値解析に必要な手法を理解するために、解析に必要な計算手法を利用したシミュレーションを行う。

20

第2章

計算手法
2.1 計算方法
本研究では Time-Dependent Ginzburg-Landau 方程式(TDGL 方程式)を解くまでの過程を詳細に記述することで、TDGL 方程式を解く環境を構築する。今回の計算では TDGL 方程式を解くために必要な差分化、オイラー法、シンプレックテ
ィック法などの知識を得るために、波動方程式、拡散方程式、シュレディンガー方程式と段階を経て解いていくものとする。

2.1.1 波動方程式
波動方程式の基本的な式
𝜕 2 𝑢(𝑥, 𝑡) 𝜕 2 𝑢(𝑥, 𝑡)
(2.1)
=
𝜕𝑡 2
𝜕𝑥 2
を差分化することを考える。𝑥方向の微小変位をΔ𝑥、時間軸方向の微小変位をΔ𝑡
とおいて𝑢(𝑥, 𝑡)をテーラー展開すると
𝑢(𝑥 + Δ𝑥, 𝑡) = 𝑢(𝑥, 𝑡) +

𝜕𝑢
1 𝜕2𝑢
1 𝜕2𝑢
2
(
)
(Δ𝑥 )3 + ⋯
Δ𝑥 +
Δ𝑥
+
𝜕𝑥
2! 𝜕𝑥 23! 𝜕𝑥 2

𝜕𝑢
1 𝜕2𝑢
1 𝜕2𝑢
2
(
)
(Δ𝑥 )3 + ⋯
Δ𝑥 +
Δ𝑥
−
𝜕𝑥
2! 𝜕𝑥 23! 𝜕𝑥 2となり、これらを足し合わせると𝑢(𝑥 − Δ𝑥, 𝑡) = 𝑢(𝑥, 𝑡) −𝜕2 𝑢
1[𝑢(𝑥 + Δ𝑥, 𝑡) − 2𝑢(𝑥, 𝑡) + 𝑢(𝑥 − Δ𝑥, 𝑡)] − 𝑂(Δ𝑥 2 )=𝜕𝑥 2 Δ𝑥 2

(2.2)
(2.3)

(2.4)

時間軸方向についても同様に考えて𝜕2 𝑢
1= 2 [𝑢(𝑥, 𝑡 + Δ𝑡) − 2𝑢(𝑥, 𝑡) + 𝑢(𝑥, 𝑡 − Δ𝑡)] − 𝑂(Δ𝑡 2)
2
𝜕𝑡
Δ𝑡

(2.5)

以上の式を(2.1)式に代入して整理すると差分化された式
𝑢(𝑥, 𝑡 + Δ𝑡) =Δ𝑡 2{𝑢(𝑥 + Δ𝑥, 𝑡) + 𝑢(𝑥 − Δ𝑥, 𝑡) − 2𝑢(𝑥, 𝑡)} − 𝑢(𝑥, 𝑡 − Δ𝑡) + 2𝑢(𝑥, 𝑡)
Δ𝑥 2
となる。
このときの境界条件は
21

(2.6)

𝜕𝑢
𝜕𝑢
(0, 𝑡) =(𝑥, 𝑡) = 0
𝜕𝑥
𝜕𝑥

(2.7)

𝑢(𝑥, 𝑡) = exp(−𝑥 2 )(2.8)

で、初期条件はである。この条件下で数値計算を行った。

2.1.2 拡散方程式
拡散方程式の基本的な式
𝜕𝑢 (𝑥, 𝑡) 𝜕 2 𝑢(𝑥, 𝑡)
(2.9)
=
𝜕𝑡
𝜕𝑥 2
を差分化することを考える。時刻𝑡に関して前進差分空間𝑥に関して中心差分を行うと𝑢(𝑥, 𝑡 + Δ𝑡) − 𝑢(𝑥, 𝑡) 𝑢(𝑥 + Δ𝑥, 𝑡) − 2𝑢(𝑥, 𝑡) + 𝑢(𝑥 − Δ𝑥, 𝑡)
=
Δ𝑡
Δ𝑥 2
となり、この式を変形すると(2.10)Δ𝑡(𝑢(𝑥 + Δ𝑥, 𝑡) − 2𝑢(𝑥, 𝑡) + 𝑢(𝑥 − Δ𝑥, 𝑡)) (2.11)
Δ𝑥 2となる。右辺は時刻𝑡における値であるのに対し、左辺は時刻𝑡 + Δ𝑡の値、すな
わち 1 ステップ後の値だけで表すことができている。従って、時刻𝑡における各格子点の値がわかっていれば𝑡 + 1における値を計算可能である。
このときの境界条件は
𝑢(𝑥, 𝑡 + Δ𝑡)＝𝑢(𝑥, 𝑡) +𝜕𝑢
𝜕𝑢(0, 𝑡) =(𝑥, 𝑡) = 0𝜕𝑥𝜕𝑥(2.12)

で、初期条件は𝑢(𝑥, 𝑡) = exp(−𝑥 2 )である。この条件下で数値計算を行った。

(2.13)

2.1.3 シュレディンガー方程式
一次のシュレディンガー方程式
𝜕𝜓
𝜕 2𝜓
(2.14)= − 2 + 𝑉′(|𝜓|2) 𝜓𝜕𝑡𝜕𝑥に関して差分化することを考える。ここでの𝜓 = 𝜓(𝑥, 𝑡)は複素従属変数、𝑉(𝑥)はi1

ポテンシャルを表す関数である。今回は𝑉(𝑥 ) = 2 𝜔𝑥 2 である。
22

まず差分法を適用させるために方程式(2.14)の右辺第一項の微分項を差分項に置き換える。空間の微小変位をΔ𝑥、離散変数を𝜓𝑖 = 𝜓(𝑖𝛥𝑥)として、中心差分を行うと、
𝜕𝜓𝑖+1 − 2𝜓𝑖 + 𝜓𝑖−1 𝜕 2 𝜓(2.15)𝜓𝑖 ==+ 𝑜(Δ𝑥 2 )
𝜕𝑡
Δ𝑥 2
𝜕𝑥 2
となる。シンプレクティック解法は任意のハミルトン系に対して構成できる。
一方、ハミルトニアンが𝐻 = 𝐻1 + 𝐻2 のように 2 つの和で書かれていて、𝐻1 , 𝐻2がそれぞれ単独で解が書ける場合、この計に対してシンプレクティック解法を適用することができる。1 次の解法は𝐻1 をハミルトニアンとしてΔ𝑡だけ時間発展させ、次に𝐻2 をハミルトニアンとしてΔ𝑡だけ時間発展させたものの合成として 1i

時間ステップが定義される。これは非可換な作用素𝐴, 𝐵の指数関数に関する関係式
exp[(𝐴 + 𝐵)𝛥𝑡] = exp[𝐴𝛥𝑡] exp[𝐵𝛥𝑡] + 𝑂(𝛥𝑡 2 )(2.16)
が基礎となっている。ここで𝐻1 、𝐻2 とのポアソン括弧を作る微分作用素を(2.17)𝐴 ≡ {・, 𝐻1 }, 𝐵 ≡ {・, 𝐻2 }と定義すると時間ステップΔ𝑡に対する 1 次のシンプレクティック数値解法は
𝑆1 (𝛥𝑡) = exp[𝛥𝑡𝐴]exp[𝛥𝑡𝐵](2.18)
で与えられる。同様にして時間反転対称な 2 次のシンプレクティック数値解法
は、
𝑆2 (𝛥𝑡) = exp [𝛥𝑡𝐴𝛥𝑡𝐴] exp[𝛥𝑡𝐵] exp []22(2.19)

である。今回の数値解析ではこの 2 次のシンプレクティック法を利用している。
(2.19)式を𝜓の実部、虚部をそれぞれ𝑞, 𝑝として適用させると、
𝛥𝑡𝜓𝑝
𝛥𝑡𝐴𝜓𝑝
𝑞(𝑡 + 𝛥𝑡)(
) = exp [] exp[𝛥𝑡𝜓𝑞 ] exp []𝑝(𝑡 + Δ𝑡)22(2.20)となり、実際に解く式は
11
1 𝜕 2 𝜓𝑝 1 2
𝑞 (𝑡 + 𝛥𝑡) = 𝑞 (𝑡) + Δ𝑡 (−+ 𝜔𝑥 𝑝(𝑡))222 𝜕𝑥 22(2.21)

1 𝜕 2 𝜓𝑞 1 2
1(𝑡𝑝(𝑡 + 𝛥𝑡) = 𝑝(𝑡) − Δ𝑡 (−
+
𝜔𝑥
𝑞
+
Δ𝑡))
2 𝜕𝑥 22
2(2.22)

1
1 𝜕 2 𝜓𝑝 1 2
𝑞 (𝑡 + 𝛥𝑡) = 𝑞 (𝑡 + Δ𝑡) + Δ𝑡 (−+ 𝜔𝑥 𝑝(𝑡 + Δ𝑡))2
2 𝜕𝑥 22(2.23)の 3 式を計算する。
23

今回の各パラメータは𝛥𝑥 = 0.0500, 空間メッシュ数𝐿 = 10.0, 𝛥𝑡 = 0.000500,
全 Time step 𝑁 = 200とする。境界条件は
𝜓(0) = 𝜓(𝑥 ) = 0(2.24)
である。初期条件は(𝑥 − 𝜇 )2
1
exp (−)(2.25)2𝜎 2√2𝜋𝜎の正規分布である。

2.1.4 TDGL 方程式
まず、TDGL 方程式におけるオーダーパラメータ𝛹を求めるために(1.27)式について考える。(1.27)式を時間𝑡に関して差分を行うと
𝛹 (𝑡 + 𝛥𝑡) − 𝛹(𝑡)+ i𝑉(𝑡)𝛹(𝑡) − 𝛻 2 𝛹(𝑡) + 2i𝑨･𝛻𝛹(𝑡)𝛥𝑡+ |𝑨|2 𝛹(𝑡) − 𝛹 (𝑡) + |𝛹 (𝑡)|2 𝛹(𝑡) = 0となる。(2.26)式を実部と虚部に分解すると(2.26)

𝛹R (𝑡 + 𝛥𝑡) − 𝛹R(𝑡)− 𝑉 (𝑡)𝛹I (𝑡) − 𝛻 2 𝛹R(𝑡) − 2𝑨･𝛻𝛹I (𝑡)𝛥𝑡+ |𝑨|2 𝛹R(𝑡) − 𝛹R (𝑡) + (𝛹R(𝑡)2 + 𝛹I (𝑡)2 )𝛹R(𝑡) = 0(2.27)

𝛹I(𝑡 + 𝛥𝑡) − 𝛹I (𝑡)+ 𝑉 (𝑡)𝛹I (𝑡) − 𝛻 2 𝛹I(𝑡) + 2𝑨･𝛻𝛹R(𝑡)𝛥𝑡+ |𝑨|2 𝛹I(𝑡) − 𝛹I (𝑡) + (𝛹R (𝑡)2 + 𝛹I(𝑡)2 )𝛹I (𝑡) = 0(2.28)の 2 式で表すことができる。ここではオーダーパラメータ𝛹 = 𝛹R + i𝛹Iである。
また、𝑉はスカラーポテンシャル、Δ𝑡は時間の微小変位を表している。(2.27),(2.28)式をそれぞれ𝛹R(𝑡 + 𝛥𝑡)、𝛹I (𝑡 + 𝛥𝑡)に関して整理すると
𝛹R(𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + Δ𝑡(𝑉(𝑡)𝛹I (𝑡) + 𝛻 2 𝛹R(𝑡) + 2𝑨･𝛻𝛹I(𝑡)− |𝑨|2 𝛹R(𝑡) + 𝛹R (𝑡) − (𝛹R(𝑡)2 + 𝛹I (𝑡)2 )𝛹R(𝑡))(2.29)

𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + Δ𝑡(−𝑉(𝑡)𝛹R (𝑡) + 𝛻 2 𝛹I (𝑡) − 2𝑨･𝛻𝛹R (𝑡)− |𝑨|2 𝛹I (𝑡) + 𝛹I(𝑡) − (𝛹R (𝑡)2 + 𝛹I (𝑡)2 )𝛹I (𝑡))(2.30)となる。今回は 2 次のシンプレクティック法を参考にしたオイラー法で解くため、実際に解く式は
1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡)21+ Δ𝑡(𝑉(𝑡)𝛹I (𝑡) + 𝛻 2 𝛹R (𝑡) + 2𝑨･𝛻𝛹I (𝑡) − |𝑨|2 𝛹R(𝑡)2+ 𝛹R (𝑡) − {(𝛹R(𝑡)2 + 𝛹I (𝑡)2 )𝛹R(𝑡)})24(2.31)

𝛹I(𝑡 + 𝛥𝑡) = 𝛹I (𝑡)1+ Δ𝑡 (−𝑉(𝑡)𝛹R (𝑡 + Δ𝑡) + 𝛻 2 𝛹I (𝑡)21− 2𝑨･𝛻𝛹R (𝑡 + Δ𝑡) − |𝑨|2 𝛹I(𝑡) + 𝛹I (𝑡)2(2.32)

2
1
{(𝛹
(𝑡
−
+ Δ𝑡) + 𝛹I (𝑡)2 ) 𝛹I(𝑡)})R
2

1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡 + Δ𝑡)
2
1+ Δ𝑡 (𝑉 (𝑡)𝛹I (𝑡 + Δ𝑡) + 𝛻 2 𝛹R (𝑡 + Δ𝑡)2
1
1+ 2𝑨･𝛻𝛹I(𝑡 + Δ𝑡) − |𝑨|2 𝛹R (𝑡 + Δ𝑡) + 𝛹R (𝑡 + Δ𝑡)22(2.33)

2
1
1
{(𝛹
(𝑡
−
+ Δ𝑡) + 𝛹I (𝑡 + Δ𝑡)2 ) 𝛹R (𝑡 + Δ𝑡)})R
2
2の 3 式を順番に計算する。
次にスカラーポテンシャル𝑉を求めるために(1.34)式について考える。(1.34)
式の離散化を行うと
𝑛

𝜎∑
𝑘=1

𝑉 (𝑥𝑘 + Δ𝑙, 𝑡) − 2𝑉 (𝑥𝑘 , 𝑡) + 𝑉(𝑥𝑘 − Δ𝑙, 𝑡)Δ𝑙 2(2.34)= (𝛹R𝛻 2 𝛹I − 𝛹I 𝛻 2 𝛹R − 2𝛹R 𝑨･𝛻𝛹R − 2𝛹I𝑨･𝛻𝛹I )となる。ここで𝑛は次元数、Δ𝑙は距離の微小変位である。(2.34)式を𝑉に関して整理すると
𝑛

𝑉 (𝑡) = (∑(𝑉(𝑥𝑘 + Δ𝑙, 𝑡) + 𝑉(𝑥𝑘 − Δ𝑙, 𝑡)) −
𝑘=1

Δ𝑙 2
(𝛹R 𝛻 2 𝛹I
𝜎(2.35)− 𝛹I 𝛻 2 𝛹R − 2𝛹R 𝑨･𝛻𝛹R − 2𝛹I 𝑨･𝛻𝛹I ) ) /2𝑛となる。(2.35)式を時間発展させるために、左辺を𝛥𝑡だけ進めると、

25𝑛

𝑉 (𝑡 + Δ𝑡) = (∑(𝑉(𝑥𝑘 + Δ𝑙, 𝑡) + 𝑉 (𝑥𝑘 − Δ𝑙, 𝑡))
𝑘=1

−

Δ𝑙 2
(𝛹R 𝛻 2𝛹I − 𝛹I 𝛻 2𝛹R − 2𝛹R𝑨･𝛻𝛹R
𝜎(2.36)− 2𝛹I𝑨･𝛻𝛹I ) ) /2𝑛となる。(1.34)式はラプラス方程式であるため、この方程式を解く際にはガウス・ザイデル法を用いて解く。
上記の式(2.31) – (2.33), (2.36)を 2 次元横磁界、3 次元縦磁界の場合でそれぞれ適用することによって数値解析の式を導く。
(1)2 次元横磁界
今回は磁場侵入長より細い超伝導体を仮定したため、オーダーパラメータは外部磁界にのみ依存する変数となる。ベクトルポテンシャル𝑨は外部磁場に影響を受ける変数であるため、電流をy軸に沿って流した場合
𝑨 = 𝐵ext 𝑥𝒊𝑦(2.37)となる。𝐵ext は外部磁界を表している。これを用いるとオーダーパラメータ𝛹は1
𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)2= 𝛹R(𝑥, 𝑦, 𝑡)
1+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑡)𝛹I (𝑥, 𝑦, 𝑡)2𝜕 2 𝛹R(𝑥, 𝑦, 𝑡) 𝜕 2 𝛹R(𝑥, 𝑦, 𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑡)+(+) + 2𝐵ext 𝑥
2
2
𝜕𝑥
𝜕𝑦
𝜕𝑦
− (𝐵ext 𝑥 )2 𝛹R(𝑥, 𝑦, 𝑡) + 𝛹R (𝑥, 𝑦, 𝑡)− (𝛹R (𝑥, 𝑦, 𝑡)2 + 𝛹I (𝑥, 𝑦, 𝑡)2 )𝛹R (𝑥, 𝑦, 𝑡))26(2.38)𝛹I(𝑥, 𝑦, 𝑡 + 𝛥𝑡)= 𝛹I (𝑥, 𝑦, 𝑡)
1+ 𝛥𝑡 (−𝑉(𝑥, 𝑦, 𝑡)𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡)2
𝜕 2 𝛹I (𝑥, 𝑦, 𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑡)
+(
+
)
𝜕𝑥 2
𝜕𝑦 2
1
𝜕𝛹R (𝑥, 𝑦, 𝑡 + 2 Δ𝑡)− 2𝐵ext 𝑥− (𝐵ext 𝑥 )2 𝛹I(𝑥, 𝑦, 𝑡)𝜕𝑦(2.39)+ 𝛹I (𝑥, 𝑦, 𝑡)2
1− (𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡) + 𝛹I (𝑥, 𝑦, 𝑡)2 ) 𝛹I(𝑥, 𝑦, 𝑡))
2

𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑡)𝛹I(𝑥, 𝑦, 𝑡 + Δ𝑡)+(

1
𝜕 2 𝛹R (𝑥, 𝑦, 𝑡 + 2 Δ𝑡)𝜕𝑥 21
𝜕 2 𝛹R (𝑥, 𝑦, 𝑡 + 2 Δ𝑡)
+
)
𝜕𝑦 2

𝜕𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)
1+ 2𝐵ext 𝑥− (𝐵ext 𝑥 )2 𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡)𝜕𝑦2
1+ 𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡)
2
2
1
(𝑥,
− (𝛹R
𝑦, 𝑡 + Δ𝑡) + 𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)2 ) 𝛹R (𝑥, 𝑦, 𝑡
21
+ Δ𝑡))
2となる。また、スカラーポテンシャル𝑉は

27(2.40)

𝑉(𝑥, 𝑦, 𝑡 + Δ𝑡) = (𝑉(𝑥 + Δ𝑙, 𝑦, 𝑡) + 𝑉(𝑥 − Δ𝑙, 𝑦, 𝑡) + 𝑉(𝑥, 𝑦 + Δ𝑙, 𝑡)+ 𝑉(𝑥, 𝑦 − Δ𝑙, 𝑡)−

Δ𝑙 2
𝜕 2 𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)(𝛹R(𝑥, 𝑦, 𝑡 + Δ𝑡) (
𝜎
𝜕𝑥 2

+

𝜕 2 𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡))𝜕𝑦 2− 𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡) (+𝜕 2 𝛹R(𝑥, 𝑦, 𝑡 + Δ𝑡)
𝜕𝑥 2(2.41)

𝜕 2 𝛹R(𝑥, 𝑦, 𝑡 + Δ𝑡))
𝜕𝑦 2− 2𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡)𝐵ext𝑥− 2𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)𝐵ext 𝑥𝜕𝛹R(𝑥, 𝑦, 𝑡 + Δ𝑡)𝜕𝑦𝜕𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)
)) /4
𝜕𝑦

このときの境界条件は
𝜕𝛹R (𝑥, 𝑦, 𝑡)=0𝜕𝑥𝜕𝛹I (𝑥, 𝑦, 𝑡)=0𝜕𝑥𝜕𝛹R (𝑥, 𝑦, 𝑡)+ 𝐵ext 𝑥𝛹I (𝑥, 𝑦, 𝑡) = 0𝜕𝑦𝜕𝛹I (𝑥, 𝑦, 𝑡)+ 𝐵ext 𝑥𝛹R (𝑥, 𝑦, 𝑡) = 0𝜕𝑦𝜕𝑉 (𝑥, 𝑦, 𝑡)=0𝜕𝑥
𝜕𝑉 (𝑥, 𝑦, 𝑡)
𝐽𝑦
=−
𝜕𝑦
𝜎

(2.42)(2.43)(2.44)(2.45)(2.46)(2.47)である。また、初期条件は
𝛹 (𝑥, 𝑦, 𝑡) = cos𝜃 + i sin𝜃𝑉 (𝑥, 𝑦, 𝑡) − 𝐽𝑦/𝜎である。
2 次元横磁界における計算条件のパラメータを表 2.1 に示す。
28

(2.48)(2.49)

表 2.1:2 次元横磁界下における計算条件のパラメータ
試料サイズ

10×30

電流密度𝐽

0.2

磁束密度𝐵

0.4

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1(2)3 次元横磁界
ベクトルポテンシャル𝑨は 2 次元の時と同様に考えて
𝑨 = 𝐵ext 𝑥𝒊𝑦(2.50)となる。オーダーパラメータ𝛹は1
𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
+(
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2
+ 2𝐵ext 𝑥𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)− (𝐵ext 𝑥 )2 𝛹R(𝑥, 𝑦, 𝑧, 𝑡)𝜕𝑦+ 𝛹R(𝑥, 𝑦, 𝑧, 𝑡)− (𝛹R(𝑥, 𝑦, 𝑧, 𝑡)2 + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2 )𝛹R (𝑥, 𝑦, 𝑧, 𝑡))

29(2.51)

𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)= 𝛹I(𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (−𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2+(

𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2

1
𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)− 2𝐵ext 𝑥− (𝐵ext 𝑥 )2 𝛹I(𝑥, 𝑦, 𝑧, 𝑡)𝜕𝑦(2.52)+ 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
2
1− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) + 𝛹I(𝑥, 𝑦, 𝑧, 𝑡)2 ) 𝛹I(𝑥, 𝑦, 𝑧, 𝑡))2

𝛹R(𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)+(

1
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)
𝜕𝑥 2

+

1
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)𝜕𝑦 21
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)
+
)
+
2𝐵
𝑥
ext
𝜕𝑧 2
𝜕𝑦
1
1
− (𝐵ext 𝑥 )2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) + 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2
2
2
1− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2

1+ 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2 ) 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡))2となる。また、スカラーポテンシャル𝑉は30(2.53)

𝑉 (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)= (𝑉(𝑥 + Δ𝑙, 𝑦, 𝑧, 𝑡) + 𝑉 (𝑥 − Δ𝑙, 𝑦, 𝑧, 𝑡)+ 𝑉 (𝑥, 𝑦 + Δ𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦 − Δ𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 + Δ𝑙, 𝑡)+ 𝑉(𝑥, 𝑦, 𝑧 − Δ𝑙, 𝑡)
Δ𝑙 2
𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)−(𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) (
𝜎
𝜕𝑥 2
+

𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
+
)
𝜕𝑦 2
𝜕𝑧 2

− 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) (+(2.54)

𝜕 2 𝛹R(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
𝜕𝑥 2

𝜕 2 𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
+
)
𝜕𝑦 2
𝜕𝑧 2

− 2𝛹R(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑥− 2𝛹I(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑥

𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝜕𝑦𝜕𝛹I(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
)) /6
𝜕𝑦

このときの境界条件は
𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑥𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑥𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡)+ 𝐵ext 𝑥𝛹I (𝑥, 𝑦, 𝑧, 𝑡) = 0
𝜕𝑦(2.55)(2.56)(2.57)

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)+ 𝐵ext 𝑥𝛹R (𝑥, 𝑦, 𝑧, 𝑡) = 0
𝜕𝑦(2.58)

𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑧

(2.59)

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑧

(2.60)

𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡) 𝜕𝑉 (𝑥, 𝑦, 𝑧, 𝑡)
=
=0
𝜕𝑥
𝜕𝑧
𝜕𝑉 (𝑥, 𝑦, 𝑧, 𝑡)
𝐽𝑦
=−
𝜕𝑦
𝜎
31

(2.61)(2.62)

である。また、初期条件は
𝛹 (𝑥, 𝑦, 𝑧, 𝑡) = cos𝜃 + i sin𝜃𝑉 (𝑥, 𝑦, 𝑧, 𝑡) = −𝐽𝑦/𝜎である。

(2.63)(2.64)

3 次元横磁界における計算条件のパラメータを表 2.1 に示す。
表 2.2:3 次元横磁界下における計算条件のパラメータ
試料サイズ

10×30×10

電流密度𝐽

0.2

磁束密度𝐵

0.4

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1(1)3 次元縦磁界
ベクトルポテンシャル𝑨は電流に平行に磁界を与える縦磁界では𝐵ext 𝑧𝒊𝑥 𝐵ext 𝑥𝒊𝑍−22となる。これを用いるとオーダーパラメータ𝛹は𝑨=(2.65)1
𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
+(
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2
𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)− 𝐵ext 𝑥
𝜕𝑥𝜕𝑧
2
2
𝐵ext 𝑧𝐵ext 𝑥) +() ) 𝛹R(𝑥, 𝑦, 𝑧, 𝑡) + 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)− ((22+ 𝐵ext 𝑧− (𝛹R(𝑥, 𝑦, 𝑧, 𝑡)2 + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2 )𝛹R (𝑥, 𝑦, 𝑧, 𝑡))

32(2.66)

𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)= 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (−𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2+(

𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2

− 𝐵ext 𝑧

1𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)
𝜕𝑥

𝐵ext 𝑥1𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)𝜕𝑧

𝐵ext 𝑧 2
𝐵ext 𝑥 2) +() ) 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
− ((
2
2
2
1
− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2 ) 𝛹I (𝑥, 𝑦, 𝑧, 𝑡))233(2.67)

𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)+(

1
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)
𝜕𝑥 2

+

1
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 2 Δ𝑡)𝜕𝑦 21
𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑡 + Δ𝑡)
2
+
)
+
𝐵
𝑧
ext
𝜕𝑧 2
𝜕𝑥
𝜕𝛹I(𝑥, 𝑦, 𝑡 + Δ𝑡)𝜕𝑧
2
𝐵ext 𝑥𝐵ext 𝑥 2
1) +() ) 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
− ((
2
2
2
− 𝐵ext 𝑥

1+ 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2
2
1− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2

1+ 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)2 ) 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡))2となる。また、スカラーポテンシャル𝑉は

34(2.68)

𝑉 (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)= (𝑉(𝑥 + Δ𝑙, 𝑦, 𝑧, 𝑡) + 𝑉 (𝑥 − Δ𝑙, 𝑦, 𝑧, 𝑡)+ 𝑉 (𝑥, 𝑦 + Δ𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦 − Δ𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 + Δ𝑙, 𝑡)+ 𝑉(𝑥, 𝑦, 𝑧 − Δ𝑙, 𝑡)
Δ𝑙 2
𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)−(𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) (
𝜎
𝜕𝑥 2
+

𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) 𝜕 2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
+
)
𝜕𝑦 2
𝜕𝑧 2

− 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡) (+𝜕 2 𝛹R(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
𝜕𝑥 2(2.69)

𝜕 2 𝛹R (𝑥, 𝑦, 𝑡 + Δ𝑡) 𝜕 2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
+
)
𝜕𝑦 2
𝜕𝑧 2

− 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑧𝜕𝛹R(𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
𝜕𝑥+ 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑥

𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝜕𝑧− 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑧

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
𝜕𝑥+ 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)𝐵ext 𝑥

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + Δ𝑡)
)) /6
𝜕𝑧

このときの境界条件は
𝜕𝛹R(𝑥, 𝑦, 𝑧, 𝑡) 𝐵ext 𝑧+
𝛹I (𝑥, 𝑦, 𝑧, 𝑡) = 0𝜕𝑥2𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝐵ext 𝑧−
𝛹R (𝑥, 𝑦, 𝑧, 𝑡) = 0𝜕𝑥2𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑦

(2.70)(2.71)(2.72)

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡)=0𝜕𝑦

(2.73)

𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝐵ext 𝑥−
𝛹I (𝑥, 𝑦, 𝑧, 𝑡) = 0𝜕𝑧2(2.74)

𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝐵ext 𝑥+
𝛹R (𝑥, 𝑦, 𝑧, 𝑡) = 0𝜕𝑧2(2.75)35𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡) 𝜕𝑉 (𝑥, 𝑦, 𝑧, 𝑡)
=
=0
𝜕𝑥
𝜕𝑧
𝜕𝑉 (𝑥, 𝑦, 𝑧, 𝑡)
𝐽𝑦
=−
𝜕𝑦
𝜎

(2.76)(2.77)である。また、初期条件は
𝛹 (𝑥, 𝑦, 𝑧, 𝑡) = cos𝜃 + i sin𝜃𝑉 (𝑥, 𝑦, 𝑧, 𝑡) = −𝐽𝑦/𝜎である。

(2.78)(2.79)

3 次元縦磁界における計算条件のパラメータを表 2.1 に示す。
表 2.3:3 次元縦磁界下における計算条件のパラメータ
試料サイズ

10×30×10

電流密度𝐽

0.2

磁束密度𝐵

0.4

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1

TDGL 方程式を解く際に得られた結果から計算の安定性を確認するためにエネルギーに関して注目する。それぞれポテンシャルエネルギー𝐸P と運動エネル
ギー𝐸K として
1
𝐸P = −|𝛹 |2 + |𝛹 |4
2
𝐸K = |(−i∇ − 𝑨)𝛹 |2(2.80)(2.81)とする。
また、TDGL 方程式の計算結果を表示する際に仮定される状況を図 2.1、図 2.2に示す。図 2.1 は 2 次元の場合を仮定し、横磁界における状況を再現するため電流と磁界を垂直に印加する。図 2.2 は電流と磁界を垂直に印加し横磁界を、平行に印加し縦磁界を再現するものとする。

36

𝐵

𝐽

𝑦

𝑥
図 2.1: 2 次元横磁界の模式図

𝐵(transverse)

図 2.2: 3 次元の模式図

37

また超伝導体に人工ピンを添加した際のピンの座標は(𝑥, 𝑦, 𝑧) = (±2, ±5, ±2)としており、その様子を示した図を以下の図 2.3、図 2.4 に示す。
今回は人工ピンを添加することにより、横磁界下における場合と、縦磁界下における場合とでどのような磁束線の挙動の変化に関して比較をする。ピンの配置に関しては、磁束線が通過する箇所に配置することによって、変化に影響を与えないということがないようにした。また、ピンの数としては縦磁界下においては縦方向により多くのピンを配置することによってピン止めされる磁束線が増える可能性があるため、ピンとピンの間でどのような変化を起こすかという点に注目したい。そのため縦方向、横方向に 2 個ずつ配置することにした。
また、計算においてはピンの箇所を|𝛹 |2 = 0と強制的にしている。これはこの位置に量子化磁束を常に存在させることによって、ピンとして磁束線の動きを止める役割を果たしている。

図 2.3: 超伝導体を𝑦軸正方向から見た図

ピン

図 2.4: 超伝導体を𝑧軸正方向から見た図

38

第3章

実行結果
3.1 波動方程式
プログラムを実行した結果を以下の図 3.1 に示す。
(a)は初期状態を表していて時間が経過するにつれて初期状態の波が発散していく。(b)は境界にぶつかり反射したあとの状態である。
𝑢

𝑢

1
0.01
−1 0 1

𝑥 −8

8

0

(a)

𝑥

(b)

図 3.1:波動方程式の計算結果。(a)初期状態𝑡 = 0。(b)境界にぶつかり反射したあとの状態𝑡 = 0.2。

3.2 拡散方程式
プログラムを実行した結果を以下の図 3.2 に示す。
(a)から(c)に向かうにつれて徐々に波形が拡散している様子を示している。
𝑢
1
−1

0
(a)

1

𝑥

−2

𝑢

𝑢

0.5

0.33

0
(b)

2

𝑥 −3

0
(c)

図 3.2:拡散方程式の計算結果。𝑡 = (a)0。(b)10。(c)30

39

3

𝑥

3.3 シュレディンガー方程式
プログラムを実行した結果を以下の図 3.3 に示す。
図 3.3 の(a)を見ると、初期条件の正規分布の形が実現されていることがわかる。時間が経過して(b)、(c)の状態になると、粒子の分布の集中しており、一定の場所に粒子が集まっていることがわかる。(c)の状態からさらに時間が経過す
ると粒子の分布は境界に衝突して反射し、(d)のように徐々に時間が経過し、(a)の状態と同じようになる。さらに時間を経過させると、粒子の分布が反対側に移動し、同様の挙動を示すことが分かっている。

𝜓

(a) 𝜓

𝜓

𝑥

𝑥

(c)
0

𝜋

(b) 𝜓

𝑥

(d)

𝑥

2𝜋

図 3.3:シュレディンガー方程式の計算結果。𝑡 = (a)0, (b) 800, (c) 1600,(d) 2000。

40

3.4 横磁界下における 2 次元 TDGL 方程式
横磁界下における 2 次元 TDGL 方程式のプログラムを実行した結果を以下の
図 3.4 に示す。
図 3.4 の左図は位相を、右図はオーダーパラメータの分布を表している。オー
ダーパラメータは黒点になっている部分はオーダーパラメータが低くなっていて、色が薄い部分はオーダーパラメータが高くなっている。オーダーパラメー
タが低い部分は量子化磁束が侵入しているということを示している。

41

0

(a)

(b)

(c)

(d)
𝜋

2𝜋
1

0

図 3.4:横磁界下における 2 次元 TDGL 方程式の計算結果。𝑡 =(a)5, (b)12, (c)65,(d)300。
42

3.5 横磁界下における 3 次元 TDGL 方程式
横磁界下における 3 次元 TDGL 方程式のプログラムを実行した結果を以下の
図 3.5 に示す。
今回は 3.4 節で利用したものと同じパラメータを用いて計算を行ったため、2次元で示した場合と 3 次元で示した場合にどのような対応があるかについて議論する。したがって、3.4 節で示した結果にある、磁束線の動きに変化があった
時間に焦点を合わせている。

43

(a)

(b)

(c)

(d)

𝜋
0
2𝜋
図 3.5:横磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 =(a)5, (b)12, (c)20,(d)65。
44

3.6 縦磁界下における 3 次元 TDGL 方程式
縦磁界下における 3 次元 TDGL 方程式のプログラムを実行した結果を以下の
図 3.6、図 3.7 に示す。
図 3.6 を見ると𝑡 = 10のときは磁束線が 4 本存在しているが、時間が経過し
𝑡 = 30程度になるまでに磁束線が出ていく様子がわかる。また、45 < 𝑡 < 80の区間でオーダーパラメータ|𝛹 |2 が大きく上下している。そして徐々に 3 本の磁束線が形作られるようになった。
図 3.7 を見ると、𝑡 = 150で安定してきて磁束線は 3 本らせん構造をとっており、時間の経過とともに磁束線が回転していることが位相の変化を見るとわかる。また、𝑡 = 150の時と比べて、時間の経過によって磁束線のねじれがよりきつくなっていく様子が観察できた。また、𝑡 = 270付近で磁束線の組み換えが発
生した。

45

0

(a)

(b)

(c)

(d)
𝜋

2𝜋

図 3.6:縦磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 = (a)10, (b) 30, (c)45,(d)60。
46

(a)

(b)

(c)

(d)

0

𝜋

2𝜋

図 3.7:縦磁界下における 3 次元 TDGL 方程式の計算結果。𝑡 = (a)150, (b) 250,(c)270, (d)300。
47

3.7 ピンを添加した横磁界下における TDGL 方程式
ピンを添加した場合の横磁界下における 3 次元 TDGL 方程式のプログラムを実行した結果を以下の図 3.8 に示す。
図中の丸い粒上のものをピンとしている。磁束線の動きに関しては、ピンを入れていない場合の時と同様に磁束線の動きが安定してからは磁束線が右から左へ動いている様子を見ることができた。また、ピンを入れた個所を磁束線が通過した際に磁束線がピンに引っ張られている様子を見ることができた。磁束線がピンを通過し終わった後は磁束線の動きが少し加速していた様子が見られた。また磁束線の動きとしては周期的に同じ挙動を示していた。

(a)

(b)
0

𝜋

2𝜋

図 3.8:ピンを添加した場合の横磁界下における 3 次元 TDGL 方程式の計算結果。
𝑡 = (a)150, (b) 155。
48

3.8 ピンを添加した縦磁界下における TDGL 方程式
ピンを添加した場合の縦磁界下における 3 次元 TDGL 方程式のプログラムを実行した結果を以下の図 3.9 に示す。
ピンを追加した場合はピンとピンの間に存在する磁束線が固定されていることが観測できた。ピンを入れていない場合と同様に磁束線は 3 本で構成されている磁束線はピンで止められていて、ピンを起点として磁束線の回転が行われている。そのため 3 本の磁束線が組み変わっている様子が見られた。

49

0

(a)

(b)

(c)

(d)
𝜋

2𝜋

図 3.9:ピンを添加した場合の縦磁界下における 3 次元 TDGL 方程式の計算結果。
𝑡 = (a)150, (b) 155, (c)160, (d) 165。
50

第4章

考察
4.1 波動方程式
図 3.1(a)を見ると初期条件の𝑢(𝑥, 𝑡) = exp(−𝑥 2 )の形が実現できている。今回
𝜕𝑢

𝜕𝑢の境界条件は𝜕𝑥 (0, 𝑡) = 𝜕𝑥 (𝑥, 𝑡) = 0となっておりこれはノイマン型境界条件である。ノイマン型境界条件とは境界面に変数の勾配、すなわち変化率を与える
条件であり、波に対しては自由端反射の条件を与える。したがって、初期の波の形が境界にぶつかり反射した際に元の波の形を保持したまま反射することが正しい挙動であるということが分かる。図 3.1(b)を見ると境界に反射した後の波も谷が先行した波の形が観察することができた。したがって計算結果が正しいということがわかる。

4.2 拡散方程式
図 3.2(a)の初期状態から時間が経過するにつれて(b)のように次第に拡散している様子が見てわかる。さらに時間が経過して(c)の状態になると初期値として与えられた山はほぼ視認できない程度になり、一様に拡散しきってなだらかな山のような形まで拡散していることから計算結果が正しいと言える。

4.3 シュレディンガー方程式
今回はポテンシャル𝑉(𝑥 ) = 𝜔𝑥 2 /2としてシミュレーションを行った。この𝑉(𝑥 )の値は下に凸の 2 次関数を表していて𝑥の増加量に対応してポテンシャルも増加することが分かる。したがって、井戸型ポテンシャルのような境界に接す
るまでは一定なポテンシャルとは違い、初期値からの距離が離れることによってポテンシャルの値が上昇する。そのため図 3.3(c)を見ると、粒子の分布は境界にぶつかる前から高い存在確立を示し始めていると考えられる。

4.4 横磁界下における 2 次元 TDGL 方程式
図 3.4 の左図と右図を対応させて見ると、量子化磁束の中心部分で位相の特異点、つまり位相が大きく変化している点があるということがわかる。さらに、隣り合った磁束線間で反発力が働くことによって磁束線が平行に並んだ規則的な三角格子を形成している様子がわかる。
51

また、時間に依存するエネルギーの変化の表す図を図 4.1 に示す。ここで𝐸K は式(2.80)から、𝐸P は式(2.81)から求めることができる。縦軸はエネルギーの変化を表し、横軸は計算回数を表している。今回は時間離散幅を 0.0005 間隔にして60 万回の計算を実行したため、規格化した時間で 300 までの結果が出ている。
図 4.1 を見ると、初めに𝐸K は高い値を示していたが 0 に落ち着き安定した値になった。その後𝐸K も𝐸P も𝑡 = 65付近で一旦落ち着いていた値から変化している
ことがわかる。その後はそれぞれのエネルギーは周期的な振動を初めて安定した状態が計算終了時まで続いている。このグラフと図 3.4 を対応させて考える。
図 3.4(b)の𝑡 = 12付近で形は安定しておらず、磁束線付近のオーダーパラメータも高かったが三角格子を組んでいる様子が確認できた。図 4.1 を見ると確かにエネルギーは 0 から始まり一度安定した状態になっている。図 3.4(c)の場合は
𝑡 = 65付近から計算終了まで安定した三角格子を形成している様子が確認できていたことと、一度エネルギーが増減してから周期的に安定になったことと一致しているということがわかる。

100

energyEK0EK+ EP
EP

−100
0

100

200

time
図 4.1:2 次元横磁界におけるエネルギーの時間依存性
52300

4.5 横磁界下における 3 次元 TDGL 方程式
3 次元の場合も 2 次元の時と同様にローレンツ力に従って磁束線が動く様子を確認することができた。
また、時間に依存するエネルギーの変化の表す図を図 4.2 に示す。エネルギーの変化を 2 次元の場合と比較すると、まずエネルギーの大きさ自体が 10 倍になっている。そのほかには 2 次元のときは安定するまでに 2 度の変化が確認できたが 3 次元の場合はあまりエネルギーの増減は見られないまま安定な状態へ移行した。
図 3.4 と図 3.5 を見ると、(a)のとき 2 次元では明らかに磁束線が侵入している箇所とは別に暗みがかった様子が確認できる。これは 3 次元の場合でも磁束線が安定していない状態で各所のオーダーパラメータが低い状態になっていたことから共通した結果になっていることがわかる。ただ、2 次元より 3 次元の場
合のほうが安定するまでにかかる時間が見て取れるが、これはエネルギーに関する図 4.2 を見れば明らかなことがわかる。

1000

energyEK0EK+ EP

EP

−1000
0

100

200

300

time
図 4.2:3 次元横磁界におけるエネルギーの時間依存性

53

4.6 縦磁界下における 3 次元 TDGL 方程式
時間に依存するエネルギーの変化の表す図を図 4.3 に示す。この図を見ると、20 < 𝑡 < 120の区間内でエネルギーの激しい増減や振動が目立っている。しかし、
𝑡 = 150を越えたあたりからエネルギーが安定している。エネルギーの激しい増減や振動に関しては電流密度𝐽 = 0.20と、対破壊電流密度に近い値に設定してシミュレーションを行っているため、初期状態の影響が強く残り、安定するまでに時間がかかったのではないかと考えられる。𝑡 = 150を越えたあたりから安定した状態になっているのは図 3.7 を見ると一致していることがわかる。
また、図 3.7 に関して𝑡 = 150を越えるとエネルギーは安定しているはずなのに
𝑡 = 270付近で磁束線の結合と分離が行われていた。この運動は一見、磁束カッティングモデルの分断と再結合のモデルの intra-cutting のように見えるが、今回の計算においては磁場侵入長より細い超伝導体を仮定する、細線近似を行ったため自己磁界を考慮していない。そのため、磁束カッティングモデルでは直線の磁束線のみが存在することが期待され、force-free モデルでは磁束線の組み換えは起こらないということが期待されるはずである。これに関して、他の計算方法を用いても同様の結果が得られるということが現在わかっている[17]。

54

1000

energyEK0EK+ EP
EP

−1000
0

100

200

300

time
図 4.3:3 次元縦磁界におけるエネルギーの時間依存性

4.7 ピンを添加した場合の 3 次元 TDGL 方程式
ピンを添加した場合は横磁界、縦磁界ともにピンの影響により磁束線の動きがピンによって止められる様子が観測できた。縦磁界下においては横磁界下に比べてピンとピンの間でより磁束線が止められていたため、より多くの磁束線が止まっていたということが分かる。従って、Josephson の関係式、𝑬 = 𝑩 × 𝒗より、磁束の運動量が少なければ少ないほど電磁誘導により発生する電圧が少な
くなり、より高い臨界電流密度を期待することができる。また、今回は縦に二つ並べたピンを四方に置いたため、8 つの人工ピンを添加したことになるが、縦磁界中においては縦方向に関するピンを増やすことによってより広い範囲で磁束線をピン止めすることが見込まれるため、より高い臨界電流密度が得られる可能性がある。
55

4.8 パラメータが異なる場合の実行結果
パラメータ𝐽, 𝐵を変えて計算をした際の相図を図 4.4 に示す。図中の Break は
オーダーパラメータが低すぎて磁束線の動きが観測できない状態、Dynamic は
磁束線の組み換えなどの動きが見られる状態、Static は磁束線の動きが安定してからは磁束線の動きが見られない安定した状態を示している。図を見ると、
高電流密度、高磁界になると、磁束線の動きは安定しにくくなっていることがわかる。また、磁界を大きくしたときよりも電流密度を大きくしたときの方が磁束線の挙動に与える影響が大きいことが分かる。磁束線の動きが少なければ少ないほど臨界電流密度の向上を見込むことができる。また、過剰な電流密度や磁界を与えることによって超伝導体にかかる負荷の上限を知ることができる。

0.385
Break
Dynamic
Static

J

0.3

0.2

0.1

0
0

0.5
B
図 4.4:パラメータが異なる場合の相図

56

1

第5章

まとめ
本研究では波動方程式、拡散方程式、シュレディンガー方程式、2 次元横磁界、
3 次元横磁界、3 次元縦磁界の 6 つのシミュレーションを行った。既存の研究結果が少ない 3 次元縦磁界において、磁場侵入長よりも細い超伝導体を仮定する
細線近似と、パラメータの規格化を行うことにより、TDGL 方程式を簡易化して計算することができた。
初めに波動方程式、拡散方程式、シュレディンガー方程式を解くことで TDGL
方程式を解くために必要な計算方法、計算環境を構築し、動作確認を行った。
その結果従来通りの正しい結果が得られた。
その計算方法を利用して 2 次元横磁界、3 次元横磁界、3 次元縦磁界とステッ
プアップすることで従来報告されていた研究の内容や発展した内容まで得ることができた。
3 次元縦磁界においてはパラメータを変更して計算を行うことで磁束線の挙動ごとの相図を得ることができた。これにより、パラメータごとにどういった
磁束線の挙動が見られるかを把握することができた。
また、より現実の超伝導体に近い条件としてピンを添加した際に起きる磁束線の挙動の違いに関しての計算を実行した。それにより、横磁界と縦磁界において、縦磁界の方がより高い臨界電流密度を得ることができるという臨界電流密度の特性の差も観察することができた。
今後の課題としては電流密度や磁束密度を変えた場合にどんな現象が起こり、どんな結果が得られるかを調査して電流や磁束が超伝導体に引き起こす影響を考えることが期待できる。

57

謝辞
本研究を進めるに当たり、熱心に指導してくださった、九州工業大学大学院
情報工学研究院電子情報工学研究系エレクトロニクス分野 小田部荘司教授に御礼申し上げます。小田部教授には卒業研究に関する助言だけではなく、研究室生活や私生活に関しても気遣っていただき、励ましの言葉をかけて頂きました。また、社会人としての常識や人付き合いの仕方などに関しても指導して頂きました。深く感謝申し上げます。
また、九州工業大学大学院情報工学研究院電子情報工学研究系エレクトロニ
クス分野 木内勝准教授にも多大なご指導を頂きました。超伝導体の知識に関する基礎的なことから卒業研究に繋がるような様々な内容を丁寧に教えて頂きました。数々のサポートに深く感謝申し上げます。
九州工業大学大学院情報工学府情報システム専攻 増田嘉道氏には、数値解析の時に利用したプログラムの使い方や、本研究に関する内容を基礎から発展した内容まで懇切丁寧に指導して頂きました。そのすべてに深く感謝申し上げます。
最後に公私共々お世話になりました小田部研究室、木内研究室の皆様に感謝申し上げます。

58

参考文献[1] 堂山昌男, 山本良一編: 「超伝導材料」東京大学出版会 (1987)[2] T. Matsushita: Flux Pinning in Superconductor. second ed. Springer. (2014)[3] Yu. F. Bychkov et al. : JETP Lett. 9 (1969) 404.
[4] G. W. Cullen, R. L. Novak: Appl. Phys. Lett.4 (1964) 147.
[5] J. R. Cave et al. : J. Phys. 39 (1978) C6-614.
[6] V. S. Vyatkin et al. : IEEE Trans. on Appl. Supercond. 25 (2015) 6606207.
[7] J. R. Clem: J. Low Temp. Phys. 38 (1980) 353.
[8] E. H. Brandt: J. Low Temp. Phys. 39 (1980) 41.
[9] T. Matsushita: J. Phys. Soc. Jpn. 54 (1985) 1054.
[10] V. L. Ginzburg, L. D. Landau, Zh. Eksp: Teor. Fiz. 20 (1950) 1064.
[11] D. Y. Vodolazov: Phys. Rev. B 88 (2013) 014525.
[12] J. F. Blackburn et al. : Philosophical Magazine 80 (2000) 1455.
[13] 佐々成正, 吉田春夫: 非線形 Schrödinger 方程式に対する symplectic 数
値解法, 日本応用数理学会論文誌: Vol.10, No.2 (2000) 119–131.
[14] H. Yoshida: Construction of higher order symplectic integrators, Phys. Lett.
A150 (1990) 262–268.
[15] 応用力学研究所研究集会報告 No.24AO-S3「非線形波動研究の最前線-構造と現象の多様性-」（研究代表者 太田 泰広）Research Institute forApplied Mechanics Kyusyu University March (2013)[16] 一野祐亮: 第 75 回応用物理学会秋季学術講演会(2015) 17p-A21-15.
[17] V. Vlasko-Vlasov et al.: Phys. Rev. B 91 (2015) 014516.

59

付録(1.12)式の導出
自由エネルギーは𝐸 = 𝐸s − 𝐸n =1𝛽𝐵2|(−iℏ∇ − 2𝑒𝑨)𝛹 |2 + 𝛼 |𝛹 |2 + |𝛹 |4 +
4𝑚
2
2𝜇0

(1)

である。𝑨について変分法を適用すると、
δ𝐸
=0
δ𝑨

(2)

であるので、′
𝜕𝐸 𝜕𝐸K1|(−iℏ∇ − 2𝑒𝑨)𝛹 |2)
=
=(
𝜕𝑨
𝜕𝑨
4𝑚
実部、虚部に分解する。

𝜕𝐸K
1 𝜕{(ℏ𝛻𝛹I − 2𝑒𝑨𝛹R)2 + (−ℏ𝛻𝛹R − 2𝑒𝑨𝛹I)2 }
=
𝜕𝑨
4𝑚 𝜕𝑨
偏微分を行う。
𝜕𝐸K
1{−4𝑒𝛹R(ℏ𝛻𝛹I − 2𝑒𝑨𝛹R ) + −4𝑒𝛹I(−ℏ𝛻𝛹R − 2𝑒𝑨𝛹I )}=𝜕𝑨
4𝑚
式を展開し、共通項でまとめる。
𝜕𝐸K𝑒= − {ℏ(𝛹R𝛻𝛹I − 𝛹I 𝛻𝛹R ) − 2𝑒𝑨(𝛹R2 + 𝛹I2 )}𝜕𝑨𝑚(6)式の第 2 項目を大きさで表し、展開する。
𝜕𝐸K
𝑒
2𝑒 2
|𝛹 |2 𝑨= − ℏ(𝛹R𝛻𝛹I − 𝛹I 𝛻𝛹R ) +𝜕𝑨𝑚𝑚(7)式の第 1 項に関して−1 = i2 として、整理する。
𝜕𝐸K iℏ𝑒2𝑒 2|𝛹 |2 𝑨=･2i(𝛹R 𝛻𝛹I − 𝛹I 𝛻𝛹R) +𝜕𝑨
2𝑚𝑚
2i(𝛹R𝛻𝛹I − 𝛹I 𝛻𝛹R ) = 𝛹 ∗ 𝛻𝛹 − 𝛹 𝛻𝛹 ∗ であるため、(2)式の第 1 項目は
𝜕𝐸K iℏ𝑒 ∗2𝑒 2(𝛹 𝛻𝛹 − 𝛹 𝛻𝛹 ∗ ) +| 𝛹 |2 𝑨
=
𝜕𝑨
2𝑚
𝑚
となる。
また、自由エネルギー𝐸に関して
𝐸 = 𝐸M + 𝐸K + 𝐸P
であるので
1 2
𝐸=
𝐵 + 𝐸K + 𝐸P
2𝜇0
60

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

ここで𝐵 = rot𝑨より
𝐸=

1(rot𝑨)2 + 𝐸K + 𝐸P
2𝜇0

(12)

また、ポテンシャルエネルギー𝐸P は𝑨に依存しないため偏微分を行うと
𝜕𝐸
1
𝜕𝐸K= rot rot𝑨 +
𝜕𝐴 𝜇0
𝜕𝑨

(13)

𝐵 = rot𝑨、𝐽 = rot𝐵/𝜇0 より
𝜕𝐸
𝜕𝐸K
=𝐽+
𝜕𝐴
𝜕𝑨
ここで、エネルギーは低い値を保持するため
0=𝐽+

𝜕𝐸K
𝜕𝑨

𝜕𝐸K
𝐽=−
𝜕𝑨

(14)

(15)

(9), (15)式より
iℏ𝑒 ∗
2𝑒 2
(𝛹 𝛻𝛹 − 𝛹 𝛻𝛹 ∗ ) −|𝛹 |2 𝑨2𝑚𝑚となる。これは(1.12)式と一致する。
𝐽=−

61

(16)