平成 27 年度
修士学位論文

時間依存 Ginzburg-Landau 方程式を用いた
縦磁界下の超伝導体内の磁束線に関する研究

増田 嘉道(学籍番号:14676130)

九州工業大学大学院 情報工学府
先端情報工学専攻 電子情報工学分野

指導教員:小田部 荘司 教授
平成 28 年 2 月 9 日

目次
第 1 章 序章...........................................................................................................................................1
1.1 はじめに......................................................................................................................................1
1.2 縦磁界効果..................................................................................................................................2
1.2.1 磁束カッティングモデル ..................................................................................................3
1.2.2 force-free モデル...................................................................................................................3
1.3 (Time-Dependent) Ginzburg-Landau 方程式 ..............................................................................6
1.4 ガウス=ザイデル法 ...................................................................................................................91.5 微分方程式の数値解法 .............................................................................................................91.5.1 オイラー法 ..........................................................................................................................91.5.2 シンプレクティック法 ................................................................................................... 10
1.6 本研究の目的 ...........................................................................................................................11
第 2 章 計算手法................................................................................................................................ 13
2.1 計算方法................................................................................................................................... 13
2.2 計算内容................................................................................................................................... 24
2.2.1 横磁界における 2 次元 TDGL 方程式 ........................................................................... 24
2.2.2 横磁界における 3 次元 TDGL 方程式 .......................................................................... 24
2.2.3 縦磁界における 3 次元 TDGL 方程式 .......................................................................... 24
第 3 章 計算結果................................................................................................................................ 26
3.1 横磁界における 2 次元 TDGL 方程式 ................................................................................. 26
3.2 横磁界における 3 次元 TDGL 方程式 ................................................................................. 30
3.3 縦磁界における 3 次元 TDGL 方程式 ................................................................................. 34
第 4 章 考察........................................................................................................................................ 42
4.1 横磁界における 2 次元 TDGL 方程式 ................................................................................. 42
4.2 横磁界における 3 次元 TDGL 方程式 ................................................................................. 43
4.3 縦磁界における 3 次元 TDGL 方程式 ................................................................................. 45
第 5 章 まとめ.................................................................................................................................... 46
5.1 結言........................................................................................................................................... 46
5.2 今後の課題............................................................................................................................... 46
謝辞 ...................................................................................................................................................... 47
参考文献.............................................................................................................................................. 48

I表目次
表 3.1:2 次元横磁界におけるパラメータ ...................................................................................... 27
表 3.2:3 次元横磁界におけるパラメータ ...................................................................................... 31
表 3.3:磁束線が動的な計算の 3 次元縦磁界におけるパラメータ............................................. 35
表 3.4:磁束線が静的な計算の 3 次元縦磁界におけるパラメータ............................................. 39

II

図目次
図 1.1:超伝導体に対し磁界と電流を平行に印加した状態 ............................................................2
図 1.2:磁束カッティングの模式図 ....................................................................................................3
図 1.3:磁束構造の歪み ........................................................................................................................4
図 1.4:磁束線の運動 ............................................................................................................................5
図 1.5:超伝導円柱におけるらせん状磁束フロー ............................................................................5
図 1.6:オーダーパラメータΨと超伝導体内における局所的な磁束密度 b の関係 ...................9
図 2.1:フローチャート図 ................................................................................................................. 15
図 2.2:2 次元の模式図 ...................................................................................................................... 25
図 2.3:3 次元の模式図 ...................................................................................................................... 25
図 3.1:2 次元横磁界におけるエネルギーの時間依存性 .............................................................. 27
図 3.2:2 次元横磁界における計算結果 (t=(a)4,(b)12) .................................................................. 28
図 3.3:2 次元横磁界における計算結果 (t=(a)150,(b)160,(c)170,(d)180) .................................. 29
図 3.4:3 次元横磁界におけるエネルギーの時間依存性 .............................................................. 31
図 3.5:3 次元横磁界における計算結果 (𝑡 = (a)8, (b)20)........................................................... 32
図 3.6:3 次元横磁界における計算結果 (𝑡 = (a)150, (b)160, (c)170, (d)180) ....................... 33
図 3.7:磁束線が動的な計算の 3 次元縦磁界におけるエネルギーの時間依存性 .................... 35
図3.8: 磁 束 線 が 動 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)11, (b)33, (c)46, (d)50) ................................................................................................... 36

図3.9: 磁 束 線 が 動 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)150, (b)250, (c)270, (d)300) ......................................................................................... 37

図 3.10:磁束線が静的な計算の 3 次元縦磁界におけるエネルギーの時間依存性 .................. 39
図3.11: 磁 束 線 が 静 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)8, (b)23, (c)50, (d)97) ...................................................................................................... 40

図3.12: 磁 束 線 が 静 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)150, (b)160, (c)170, (d)180)........................................................................................... 41

図 4.1:図 3.3(b)のフーリエ解析....................................................................................................... 42
図 4.2:(a)三角格子サンプルと(b)そのフーリエ解析 .................................................................... 43
図 4.3:エネルギーの時間依存性の比較 ......................................................................................... 44
図 4.4:2 次元横磁界における計算結果 .......................................................................................... 44III

第1章 序章
1.1 はじめに1908 年にオランダの Leiden 大学の Kamerlingh Onnes が世界で初めてヘリウムの液化に成功し、1911 年に液体ヘリウムを用いて水銀の電気抵抗を測定する過程において、4.2 K 以下
において電気抵抗が測定不能なほどに小さくなることを発見した。さらに、1933 年には Fritz
Walther Meissner によって超伝導体の完全反磁性(Meissner 効果)が発見された。この 2 つの現象を超伝導現象と言い、一定条件下で超伝導現象を示す物質を超伝導体と言う。超伝導体
は電気抵抗が 0 であることから大きな電流を損失なく通電することができると期待されていた。しかし、超伝導体はある磁界を超えると超伝導状態を維持できなくなることがわかった。超伝導状態から常伝導状態に代わる温度を臨界温度𝑇c、磁界を臨界磁界𝐵c と言う。超伝導状態になる条件はあるが、電気抵抗が 0 となることは非常に有用性が高いため、現在も研究は続けられている。
その後も超伝導現象の発現については不透明なままだった。しかし、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が液体窒素の常圧下での温度(77.3 K)を超えるYBa2 Cu3O7−𝛿 やBi2 Cr2 Ca2 Cu3O10+δ(𝛿は微小量)などの酸化物超伝導体が発見された。これらは銅酸化物超伝導体または高温超伝導体と呼ばれ、液体ヘリウムよりも安価な液体窒素中で超伝導状態となることから応用が期待されている。
超伝導体は磁気的な振る舞いの違いから第一種超伝導体と第二種超伝導体に分類される。
第一種超伝導体では臨界温度𝑇c以下の温度で超伝導状態となり完全反磁性(Meissner 効果)を示すが、外部から磁界を加えていくと、ある磁界を超えた所で Meissner 効果および超伝導状態は消失する。この境となる磁界を臨界磁界𝐵c という。一方で、第二種超伝導体では 𝑇c以下の温度で、下部臨界磁界𝐵c1と呼ばれる磁界までは Meissner 効果を示し、その磁界を超えても、第一種超伝導体とは異なり超伝導体内に一定の磁束の侵入を許すことで超伝導状態を維持する。そしてさらに磁界が増加すると上部臨界磁界𝐵c2と呼ばれる磁界で超伝導状態が消滅する。
現在工学利用されている超伝導体は主にNbTiやNb3 Sn等の金属系超伝導体(𝑇cが 10∼20 K程度)であり、冷媒に液体ヘリウムを使用する必要があるため冷却コストが非常に高い。そのため、高い𝑇cを持ち、比較的安価な液体窒素での冷却が可能な銅酸化物超伝導体が、金属系超伝導体の代替素材として注目されており、送電ケーブルや高磁界発生超伝導マグネッ
トなどへの応用が期待されている。

1

1.2 縦磁界効果
縦磁界とは、超伝導体に流している電流に対し、平行に磁界を印加した状態のことを言う。一般的に、超伝導体に流れる電流によって生じる自己磁界を含めて磁界は垂直方向となるが、図 1.1 に示すように、磁界と電流が平行な場合、垂直な場合とは異なる様々な現象が観測される。以下、その時の例を示す。

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

𝐽

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

2

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

(a)

(b)
(1)

(2)

(3)

図 1.2:磁束カッティングの模式図[1]。(a)は inter-cutting、(b)は intra-cutting の例である。

3

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

𝑱=

1
1
1
(𝒊𝐵 × ∇)𝐵 + 𝐵 rot𝒊𝐵
rot𝑩 = −
𝜇0
𝜇0
𝜇0

(1.3)

と表される。この電流は 3 つの成分から、すなわち、磁束密度の大きさ𝐵の勾配による成分、磁束の曲がりによる成分および force-free 電流成分である。第 1 の成分は磁気圧に、第 2 の成分は線長力に寄与する。すなわち、これらの成分は共に Lorentz 力に寄与する。磁気圧に関する成分は(1.3)式の第 1 項にのみ含まれ、一方、force-free 電流成分は第 2 項にのみ含まれる。そして線張力に関係する成分は両方の項から生じる。
以上のように、一般に電流が流れた時には、磁束構造は歪みをもつ。図 1.3 に各種の基本的な歪みを示す。(a)は磁束線の勾配、(b)は磁束の曲げ、(c)は force-free 電流が流れるときのせん断歪みである。それぞれにおける矢印方向にトルクが働くと考えられる。
ここで、実際の実験について考える。縦磁界下における臨界状態は force-free 電流の駆動
トルクとピン力が釣り合っている状態である。図 1.4 のような磁束線構造を考えると、臨界角δ𝜃cを超えた磁束線は、それを緩和させようと𝛿𝜃を減少させる𝑣1 方向に動く。さらに、

∂𝜃
𝜕𝜃
= 𝑣2
𝜕𝑡
𝜕𝑦

(1.4)

という、定常状態の条件を満たすためには𝑣1 の回転運動の他に、𝑣2 の並進運動も必要となる。
このとき、円柱状超伝導体の場合に期待される磁束線の動きを図 1.5 に示す。(a)は円柱の中心を通る磁束線の運動の方向を示したものであり、(b)は円柱の長さ方向の各部分での磁束線の運動の方向を示したものである。こうした磁束線運動は一様な円柱内の直線的な磁束線フローを円柱の軸の周りにねじった形となっている。
(a)

(b)

(c)

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

4

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

(b)

(a)

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

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

𝛽
|𝛹|4
2

(1.5)

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

(1.6)

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

𝐵2
2𝜇0

(1.7)

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

(1.8)

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

(1.9)

δ𝐸
=0
δ𝑨

(1.10)

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

となり、

𝑱=

iℏ𝑒 ∗
2𝑒 2
(𝛹 ∇𝛹 − 𝛹∇𝛹 ∗ ) −
|𝛹|2𝑨
2𝑚
𝑚

を得る。

6

(1.12)

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

(1.13)

δ𝐸
∂𝑨
= −𝜈
δ𝑨
∂𝑡

(1.14)となり、γ 及び νはそれぞれ𝛹と𝑨の時定数である。ここでゲージ変換を考えると、ℏ∂𝛹∂→ (ℏ + 2i𝑒𝑉) 𝛹∂𝑡∂𝑡(1.15)

∂𝑨 𝜕𝑨
→
+ ∇𝑉
∂𝑡
𝜕𝑡

(1.16)

であり、𝑉はスカラーポテンシャルである。これを(1.13)式及び(1.14)式に適用すると、∂𝛹
1(−iℏ∇ − 2𝑒𝑨)2 𝛹 + 𝛼𝛹 + 𝛽|𝛹|2 𝛹 = 0γ(+ 2i𝑒𝑉𝛹) +∂𝑡
4𝑚(1.17)

∂𝑨
iℏ𝑒 ∗
2𝑒 2
(𝛹 ∇𝛹 − 𝛹∇𝛹 ∗ ) −|𝛹|2𝑨 = 0
𝜈 ( + ∇𝑉) + ∇ × ∇ × 𝑨 +∂𝑡2𝑚𝑚(1.18)となる。
ただし、(1.17)式及び(1.18)式をそのまま解くのは困難である。そこで、今回は 2 つの簡易化を行うものとする。1 つ目は非常に細い超伝導体を仮定することで、超伝導体全体に外部磁界が進入することを仮定する[10]。これにより、𝑨が外部磁界にのみ依存する変数となるため、(1.18)式が一定となる。例えば、後述するように電流を𝑦軸に沿って流した場合、横磁界では𝑨 = 𝐵ext 𝑥𝒊𝑦 、縦磁界では𝑨 =𝐵ext 𝑧2𝐵𝑥𝒊𝑥 − ext 𝒊𝑧 となる。2 つ目は定数が多いため、規格2化を行う[11]。コヒーレンス長𝜉と磁場侵入長𝜆は
𝜉=𝜆=ℏ√4𝑚|𝛼|√2𝑒𝜇0 𝐻c√𝑚|𝛼|(1.19)(1.20)と表すことができるため、
1𝑥→𝑥𝜉

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

2𝑒𝛾
𝑉→𝑉
|𝛼|

(1.23)

𝜆√2𝜇0 𝐻c

7

𝑨→𝑨(1.24)

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

2 2 𝐻c𝐽c = ( )3 𝜆(1.27)

を用いると、(1.24)式より規格化された電流密度𝐽は2 𝟑/𝟐 1
𝐽→𝐽
( )
3
√2 𝐽c

(1.28)となる。よって、2/3√3 ≅ 0.385を超えるような電流密度を考えると、対破壊電流を超えた電流を流していることになるため、注意する必要がある。
TDGL 方程式を解く際に、𝛹と𝑉について解くが、方程式が(1.26)式しかないため、電流の発散についての式∇∙𝑱 = 0(1.29)

を第 2 式として解く。電流には超伝導電流𝐽sと常伝導電流𝐽n の 2 つがあり、それぞれ規格化
した結果、
i
𝑱s = (𝛹 ∗ ∇𝛹 − 𝛹∇𝛹 ∗ ) − |𝛹|2 𝑨
2(1.30)𝑱n = −𝜎∇𝑉(1.31)と表される。ここで、𝛹 ∗ は𝛹の共役複素数、𝜎は常伝導状態における導電率である。これより、(1.29)式は𝑱 = 𝑱s + 𝑱n(1.32)i𝜎∇2 𝑉 = (𝛹 ∗ ∇2 𝛹 − 𝛹∇2 𝛹 ∗ ) − ∇ ∙ (|𝛹|2 𝑨)
2(1.33)

を用いると、となる。図 1.6 に常伝導核周辺における𝛹と超伝導体内における局所的な磁束密度𝑏の関係を示す。図において左端や右端のように、|𝛹|が高いところでは超伝導状態が強く、その部
分では磁束線が進入しにくい。逆に図の中央部分のような、|𝛹|が低いところでは超伝導状態が弱くなって常伝導となっているため、磁束線が侵入しやすい。

8

図 1.6:オーダーパラメータΨと超伝導体内における局所的な磁束密度 b の関係[1]1.4 ガウス=ザイデル法
ガウス=ザイデル法とは、連立一次方程式を反復法で解く手法の 1 つとして知られている。
ガウス=ザイデル法は同様な反復法であるヤコビ法を改良したものである。𝑨𝒙 = 𝒃を仮定
すると、ヤコビ法では𝑖−1

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

𝑗=𝑖+1

𝑖−1

𝑛

(1.34)となっており、ガウス=ザイデル法では
𝑥𝑖𝑘+1 = (𝑏𝑖 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘+1 − ∑ 𝑎𝑖𝑗 𝑥𝑗𝑘 )𝑗=1(1.35)

𝑗=𝑖+1の形となる。反復法であるため、変化値|𝑥𝑖𝑘+1 − 𝑥𝑖𝑘 |が任意の変化値以下になるまで繰り返し
でこの計算を行う。この 2 つの方法の間の違いは1 ≤ j ≤ i − 1で用いる𝑥が更新される前後の違いのみである。また、ヤコビ法と比べて実装した際の使用容量と計算速度の両方の点で優れている。

1.5 微分方程式の数値解法
1.5.1 オイラー法
オイラー法は 1 階常微分方程式の数値解法の 1 つである。これは微分の定義である
d𝑦𝑦(𝑥 + ℎ) − 𝑦(𝑥)= 𝑓(𝑥, 𝑦) = lim
ℎ→0
d𝑥
ℎ

9

(1.36)

が元になっており、ℎが十分に小さい時は𝑦(𝑥 + ℎ) − 𝑦(𝑥)ℎ(1.37)𝑦(𝑥 + ℎ) = 𝑦(𝑥) + ℎ 𝑓(𝑥, 𝑦)(1.38)

𝑓(𝑥, 𝑦) =すなわちと書くことができる。
オイラー法は数学的に理解しやすく、アルゴリズムとしても簡単である。しかし、数値解法としては誤差が蓄積されるため精度が悪い。ℎを小さく取ることで、計算時間はかかるが、精度を上げることは可能である。また、この時の誤差は𝑂(ℎ)である。

1.5.2 シンプレクティック法
1 粒子の運動をハミルトン形式で見ることから始める。粒子の位置を𝑞、運動量を𝑝とし、この粒子のハミルトニアンを𝐻とすると、運動方程式は
d𝑞 ∂𝐻
=
d𝑡 ∂𝑝

(1.39)

d𝑝
∂𝐻
=−
d𝑡
∂𝑞

(1.40)となる。(1.39)式と(1.40)式において、ベクトルを用いて表すと、
d 𝑝𝑝(𝑞 ) = 𝐷 (𝑞 )d𝑡(1.41)となり、ここで𝐷は作用素である。これをスカラー量の微分方程式に照らし合わせると、そ
の解は
𝑝(𝑡)
𝑝(0)
(
) = exp(𝑡𝐷) (
)
𝑞(0)
𝑞(𝑡)

(1.42)となる。
次に作用素𝐷に注目する。ハミルトニアンは運動エネルギー𝐾とポテンシャルエネルギー
𝑈の和であるため、作用素である𝐷も運動エネルギーへの作用素𝐷𝐾 とポテンシャルエネルギ
ーへの作用素𝐷𝑈 の和となる。つまり、𝐷は
𝐷 = 𝐷𝐾 + 𝐷𝑈

(1.43)のように分解できる。これを(1.42)式に適用し、初期値を𝑡とすると、
𝑝(𝑡 + 𝛥𝑡)𝑝(𝑡)() = exp(𝛥𝑡(𝐷𝑈 + 𝐷𝐾 )) ()
𝑞(𝑡)𝑞(𝑡 + 𝛥𝑡)となる。ここで𝐷𝐾 と𝐷𝑈 が非可換の場合(𝐷𝐾 𝐷𝑈 ≠ 𝐷𝑈 𝐷𝐾 )、BCH 公式より

10

(1.44)1
exp(𝛥𝑡𝐷𝑈 ) exp(𝛥𝑡𝐷𝐾 ) = exp (𝛥𝑡(𝐷𝑈 + 𝐷𝐾 ) + 𝛥𝑡 2 (𝐷𝑈 𝐷𝐾 − 𝐷𝐾 𝐷𝑈 ) + ⋯ )2(1.45)となる。今回は𝛥𝑡は十分に小さいものとし、2 次以降の項は無視するものとする。これを(1.44)に適用すると、
𝑝(𝑡 + 𝛥𝑡)𝑝(𝑡)() = exp(𝛥𝑡𝐷𝑈 ) exp(𝛥𝑡𝐷𝐾 ) ()
𝑞(𝑡)𝑞(𝑡 + 𝛥𝑡)

(1.46)となる。ここで、指数関数の部分をテイラー展開すると、
𝑝(𝑡 + 𝛥𝑡)𝑝(𝑡)() = (1 + 𝛥𝑡𝐷𝑈 )(1 + 𝛥𝑡𝐷𝐾 ) ()
𝑞(𝑡)𝑞(𝑡 + 𝛥𝑡)

(1.47)である。なお、𝐷𝐾 は𝑞のみに、𝐷𝑈 は𝑝のみに作用するため、𝐷𝑈2 及び𝐷𝐾2 は 0 となる。また、同様に注意してこれを計算すると、
𝑝(𝑡)
𝑝(𝑡)
(1 + 𝛥𝑡𝐷𝐾 ) (
)=(
)
𝑞(𝑡)
𝑞(𝑡 + 𝛥𝑡)

(1.48)𝑝(𝑡)𝑝(𝑡 + 𝛥𝑡)
(1 + 𝛥𝑡𝐷𝑈 ) ()=()
𝑞(𝑡 + 𝛥𝑡)𝑞(𝑡 + 𝛥𝑡)

(1.49)となる。この時の誤差は𝑂(𝛥𝑡 2 )となる。さらに、この計算で誤差が累積することはなく、常に𝑂(𝛥𝑡 2 )でとどまる。
上記が 1 次のシンプレクティック法に対し、さらに高次のシンプレクティック法も存在する。例えば、2 次のシンプレクティック法であれば、
11𝑝(𝑡 + 𝛥𝑡)𝑝(𝑡)() = exp ( 𝛥𝑡𝐷𝐾 ) exp(𝛥𝑡𝐷𝑈 ) exp ( 𝛥𝑡𝐷𝐾 ) ()
𝑞(𝑡)𝑞(𝑡 + 𝛥𝑡)22(1.50)となり、(1.48)式及び(1.49)式と同様に、
𝑝(𝑡)
1
𝑝(𝑡)
1
(1 + 𝛥𝑡𝐷𝐾 ) (
)=(
)
𝑞(𝑡)
2
𝑞 (𝑡 + 𝛥𝑡)
2(1.51)𝑝(𝑡)𝑝(𝑡 + 𝛥𝑡)
(1 + 𝛥𝑡𝐷𝑈 ) (
1
1
)=(
)
𝑞 (𝑡 + 𝛥𝑡)
𝑞 (𝑡 + 𝛥𝑡)2
2(1.52)𝑝(𝑡)1
𝑝(𝑡 + 𝛥𝑡)1
(1 + 𝛥𝑡𝐷𝐾 ) ()=()
𝑞(𝑡 + 𝛥𝑡)2
𝑞 (𝑡 + 𝛥𝑡)
2(1.53)のように解くこととなる。また、この時の誤差は𝑂(𝛥𝑡 3 )となり、1 次の時と異なり、(1.50)式の作用素が対称となっていることから、時間反転対称性を持っていることがわかる。

1.6 本研究の目的
Ginzburg-Landau 方程式(GL 方程式)は超伝導体の現象的特性を説明することで知られてい11

る。臨界温度やコヒーレンス長、磁場侵入長などの第二種超伝導体の特性は GL 方程式を用
い て 調 査 で き る 。 加 え て 、 GL 方 程 式 の 時 間 発 展 性 を 記 述 し た Time-Dependent
Ginzburg-Landau 方程式(TDGL 方程式)が存在し、第二種超伝導体内の磁束線の動きを調査する際に最適である[15, 16]。しかし、TDGL 方程式を用いた研究の多くは、電流に対して垂直に磁界を印加する横磁界の場合である[17, 18, 19]。例えば、横磁界下における人工ピン入りの超伝導体にどのように磁束線がピン留めされるかを 3 次元の TDGL 方程式を用いて再現したものや[22]、横磁界中の超伝導体の電圧-電流密度特性を 2 次元の TDGL 方程式を解くことで計算したものがある[23]。
一方で、電流に対して平行に磁界を印加する縦磁界中において、常磁性効果や臨界電流密度の上昇などの特徴的な縦磁界効果が観測されている[1, 8, 21]。最近では、縦磁界効果を利用した超伝導直流電力ケーブルが設計され、総輸送電流は従来の超伝導ケーブルに比べて改善された[5]。
TDGL をそのまま解くことは困難だが、簡易化して解いても、よく知られた結果を得られたという研究が存在する[9, 10]。一方で、本来の TDGL 方程式を用いて縦磁界における磁束線のシミュレーションを行った研究も存在する[11, 20]。そこで本研究では、簡易化した 3
次元の TDGL 方程式を解くことによって、縦磁界下で超伝導体内において、どのように磁束線が動いているかについて調査を行う。

12

第2章 計算手法
2.1 計算方法
本研究では Time-Dependent Ginzburg-Landau 方程式(TDGL 方程式)を解くことによって、超伝導体内における磁束線の動きの評価・考察を行う。今回の計算では、磁界侵入長よりも小さい超伝導体を仮定する。これによってベクトルポテンシャル𝑨は外部磁場のみに影響さる変数となり、電流を𝑦軸に沿って流したと仮定した場合、電流に対して磁界を垂直に与える横磁界では、
𝑨 = 𝐵ext 𝑥𝒊𝑦(2.1)と表され、電流に対して平行に磁界を与える縦磁界では、
𝑨=

𝐵ext 𝑧𝐵ext 𝑥
𝒊 −
𝒊
2 𝑥
2 𝑧

(2.2)

と表される。ここで、𝐵extは外部磁界を、𝒊𝑥 、𝒊𝑦 、𝒊𝑧 はそれぞれ𝑥軸、𝑦軸、𝑧軸方向の単位ベクトルを表している。
次に(1.26)式について考える。(1.26)式を時間𝑡について前方差分を行い、第 3 項を展開す
ると、
𝛹(𝑡 + 𝛥𝑡) − 𝛹(𝑡)+ i𝑉(𝑡)𝛹(𝑡) − 𝛻 2 𝛹(𝑡) + 2i𝑨 ∙ 𝛻𝛹(𝑡) + |𝑨|2 𝛹(𝑡) − 𝛹(𝑡)
𝛥𝑡

(2.3)

+ |𝛹(𝑡)|2 𝛹(𝑡) = 0となり、さらに実部と虚部に分解すると、
𝛹R (𝑡 + 𝛥𝑡) − 𝛹R (𝑡)− 𝑉(𝑡)𝛹I (𝑡) − 𝛻 2 𝛹R (𝑡) − 2𝑨 ∙ 𝛻𝛹I (𝑡) + |𝑨|2 𝛹R (𝑡)
𝛥𝑡

(2.4)

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

(2.5)

− 𝛹I (𝑡) + (𝛹R (𝑡)2 + 𝛹I (𝑡)2 )𝛹I (𝑡) = 0の 2 つの式で表すことができる。ここで𝛹(= 𝛹R + i𝛹I )はオーダーパラメータ、𝑉はスカラ
ーポテンシャル、𝛥𝑡は時間の離散幅を示している。さらに(2.4)式を𝛹R (𝑡 + 𝛥𝑡)、(2.5)式を
𝛹I (𝑡 + 𝛥𝑡)について整理すると、
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥𝑡(𝑉(𝑡)𝛹I (𝑡) + 𝛻 2 𝛹R (𝑡) + 2𝑨 ∙ 𝛻𝛹I (𝑡) − |𝑨|2 𝛹R (𝑡)+ 𝛹R (𝑡) − (𝛹R (𝑡)2 + 𝛹I (𝑡)2 )𝛹R (𝑡))𝛹I (𝑡 + 𝛥𝑡) = 𝛹I (𝑡) + 𝛥𝑡(−𝑉(𝑡)𝛹R (𝑡) + 𝛻 2 𝛹I (𝑡) − 2𝑨 ∙ 𝛻𝛹R (𝑡) − |𝑨|2𝛹I (𝑡)+ 𝛹I (𝑡) − (𝛹R (𝑡)2 + 𝛹I (𝑡)2 )𝛹I(𝑡))

(2.6)
(2.7)

となる。シンプレクティック法は計算精度が他の方法に比べて高く、これを利用したいが、
TDGL 方程式が散逸系であるのに対し、シンプレクティック法は非散逸系のみ対応している。
13

そこで今回の計算では、2 次のシンプレクティック法を参考にしたオイラー法で解くため、実際に解く式は
1
1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡) + 𝛥𝑡(𝑉(𝑡)𝛹I (𝑡) + ∇2 𝛹R (𝑡) + 2𝑨 ∙ 𝛻𝛹I (𝑡)
2
2

(2.8)

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

(2.9)

2
1
− (𝛹R (𝑡 + 𝛥𝑡) + 𝛹I (𝑡)2 ) 𝛹I (𝑡))2

1
1
𝛹R (𝑡 + 𝛥𝑡) = 𝛹R (𝑡 + 𝛥𝑡) + 𝛥𝑡(𝑉(𝑡)𝛹I (𝑡 + 𝛥𝑡) + 𝛻 2 𝛹R (𝑡 + 𝛥𝑡)2
2
1
1+2𝑨 ∙ 𝛻𝛹I (𝑡 + 𝛥𝑡) − |𝑨|2 𝛹R (𝑡 + 𝛥𝑡) + 𝛹R (𝑡 + 𝛥𝑡)22(2.10)

2
1
1− (𝛹R (𝑡 + 𝛥𝑡) + 𝛹I (𝑡 + 𝛥𝑡)2 ) 𝛹R (𝑡 + 𝛥𝑡))2
2の 3 式を順番に解くこととなる。
次に(1.33)式について考える。初めに左辺の離散化を行うと、
𝑛

𝜎∑
𝑘=1

𝑉(𝑥𝑘 + 𝛥𝑙, 𝑡) − 2𝑉(𝑥𝑘 , 𝑡) + 𝑉(𝑥𝑘 − 𝛥𝑙, 𝑡)
𝛥𝑙 2(2.11)= (𝛹R 𝛻 2 𝛹I − 𝛹I 𝛻 2 𝛹R − 2𝛹R 𝑨 ∙ 𝛻𝛹R − 2𝛹I 𝑨 ∙ 𝛻𝛹I )となる。𝑛は次元数(2 次元ならば𝑛 = 2、3 次元ならば𝑛 = 3)、𝛥𝑙は距離の離散幅である。これを𝑉について整理すると、𝑛𝑉(𝑡) = (∑(𝑉(𝑥𝑘 + 𝛥𝑙, 𝑡) + 𝑉(𝑥𝑘 − 𝛥𝑙, 𝑡)) −
𝑘=1

𝛥𝑙 2
(𝛹R 𝛻 2 𝛹I − 𝛹I 𝛻 2 𝛹R
𝜎
(2.12)− 2𝛹R 𝑨 ∙ 𝛻𝛹R − 2𝛹I 𝑨 ∙ 𝛻𝛹I )) /2𝑛となる。これでは時間発展がないため、左辺のみ𝛥𝑡だけ進めると、𝑛𝑉(𝑡 + 𝛥𝑡) = (∑(𝑉(𝑥𝑘 + 𝛥𝑙, 𝑡) + 𝑉(𝑥𝑘 − 𝛥𝑙, 𝑡)) −
𝑘=1

𝛥𝑙 2
(𝛹R 𝛻 2 𝛹I − 𝛹I 𝛻 2 𝛹R
𝜎
(2.13)− 2𝛹R 𝑨 ∙ 𝛻𝛹R − 2𝛹I 𝑨 ∙ 𝛻𝛹I )) /2𝑛
また、(1.33)式はラプラス方程式となるため、この方程式を解く際にはガウス=ザイデル法を用いて解く。
最後に初期条件・境界条件だが、初期条件は

14

𝛹(𝑡 = 0) = cos𝜃 + i sin𝜃(2.14)

𝑉(𝑡 = 0) = −𝐽𝑦/𝜎(2.15)𝒏 ∙ (𝛻𝛹 − i𝑨𝛹) = 0(2.16)

𝛻𝑉 = −𝑱/𝜎(2.17)とし、境界条件はする。ここで𝜃は位相であり初期値が0 ≤ 𝜃 < 2𝜋の乱数、𝒏は法線ベクトルである。さらに
電位基準は𝑦 = 0において𝑉 = 0とする。以上のことをまとめたフローチャートを図 2.1 に示す。
開始

初期条件適用(2.12)、(2.13)
1

実部：𝑡 → 𝑡 + 𝛥𝑡
2

(2.8)
虚部：𝑡 → 𝑡 + 𝛥𝑡
(2.9)
1

実部：𝑡 + 𝛥𝑡 → 𝑡 + 𝛥𝑡2

(2.10)
𝑉について解く(2.11)

𝑉の変化値が規定値以下かYes
境界条件適用(2.14)、(2.15)
図 2.1:フローチャート図
15No

次に、実際に解く問題について考える。2 次元横磁界の場合、 (2.8)～(2.10)及び(2.13)に当てはめると、
1𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
2= 𝛹R (𝑥, 𝑦, 𝑡)1
𝜕2 𝛹R (𝑥, 𝑦, 𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑡)+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑡)𝛹I (𝑥, 𝑦, 𝑡) + (
+
)
2
𝜕𝑥 2
𝜕𝑦 2
+ 2𝐵ext 𝑥(2.18)

𝜕𝛹I (𝑥, 𝑦, 𝑡)− (𝐵ext 𝑥)2 𝛹R (𝑥, 𝑦, 𝑡) + 𝛹R (𝑥, 𝑦, 𝑡)𝜕𝑦− (𝛹R (𝑥, 𝑦, 𝑡)2 + 𝛹I (𝑥, 𝑦, 𝑡)2 )𝛹R (𝑥, 𝑦, 𝑡))
𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡)= 𝛹I (𝑥, 𝑦, 𝑡)
1+ 𝛥𝑡 (−𝑉(𝑥, 𝑦, 𝑡)𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)2
𝜕2 𝛹I (𝑥, 𝑦, 𝑡) 𝜕2 𝛹I (𝑥, 𝑦, 𝑡)𝜕𝛹R (𝑥, 𝑦, 𝑡)+(+) − 2𝐵ext 𝑥
2
2
𝜕𝑥
𝜕𝑦
𝜕𝑦
− (𝐵ext𝑥)2 𝛹I (𝑥, 𝑦, 𝑡) + 𝛹I (𝑥, 𝑦, 𝑡)2
1− (𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡) + 𝛹I (𝑥, 𝑦, 𝑡)2 ) 𝛹I (𝑥, 𝑦, 𝑡))2

16(2.19)

𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑡)𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
11
𝜕2 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
2
2
+(
+
)
𝜕𝑥 2
𝜕𝑦 2
1
𝜕𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
12+ 2𝐵ext𝑥− (𝐵ext 𝑥)2 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)𝜕𝑦2(2.20)

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

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

+

𝜕2 𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡))
𝜕𝑦 2(2.21)

𝜕2 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)− 𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡) (
+
)
𝜕𝑥 2
𝜕𝑦 2
− 2𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)𝐵ext 𝑥− 2𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡)𝐵ext 𝑥

𝜕𝛹R (𝑥, 𝑦, 𝑡 + 𝛥𝑡)𝜕𝑦𝜕𝛹I (𝑥, 𝑦, 𝑡 + 𝛥𝑡)
)) /4
𝜕𝑦

となる。初期条件は、
𝛹(𝑥, 𝑦, 𝑡 = 0) = cos𝜃 + i sin𝜃(2.22)𝑉(𝑥, 𝑦, 𝑡 = 0) = −𝐽𝑦/𝜎(2.23)

境界条件は17𝜕𝛹R (𝑥, 𝑦, 𝑡)=0𝜕𝑥(2.24)

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

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

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

𝜕𝑉(𝑥, 𝑦, 𝑡)=0𝜕𝑥

(2.28)

𝜕𝑉(𝑥, 𝑦, 𝑡)= −𝐽𝑦/𝜎𝜕𝑦(2.29)である。さらに 3 次元横磁界の場合には、
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 (𝑥, 𝑦, 𝑧, 𝑡))

18(2.30)

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

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

(2.31)

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

𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
11
𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2
2
+(
+
𝜕𝑥 2
𝜕𝑦 2
1
𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)2+) + 2𝐵ext 𝑦
2
𝜕𝑧
𝜕𝑦
1
1
− (𝐵ext𝑦)2𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) + 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)2
2
2
1− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)2 ) 𝛹R (𝑥, 𝑦, 𝑧, 𝑡
2

1
+ 𝛥𝑡))
219(2.32)

𝑉(𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)= (𝑉(𝑥 + 𝛥𝑙, 𝑦, 𝑧, 𝑡) + 𝑉(𝑥 − 𝛥𝑙, 𝑦, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦 + 𝛥𝑙, 𝑧, 𝑡)+ 𝑉(𝑥, 𝑦 − 𝛥𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 + 𝛥𝑙, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 − 𝛥𝑙, 𝑡)−

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

+

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

𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)− 𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) (𝜕𝑥 2+(2.33)

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

− 2𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝐵ext𝑦− 2𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝐵ext𝑦𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝜕𝑦𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
)) /6
𝜕𝑦

となる。初期条件は、
𝛹(𝑥, 𝑦, 𝑧, 𝑡 = 0) = cos𝜃 + i sin𝜃(2.34)𝑉(𝑥, 𝑦, 𝑧, 𝑡 = 0) = −𝐽𝑦/𝜎(2.35)

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

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

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

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

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

(2.40)

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

(2.41)

𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡) 𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡)
=
=0
𝜕𝑥
𝜕𝑧

(2.42)

𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡)= −𝐽𝑦/𝜎𝜕𝑦(2.43)

境界条件は20

となる。最後に３次元縦磁界においては1
𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2
𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)
+(
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2
+ 𝐵ext 𝑧(2.44)

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

𝐵ext𝑥 2𝐵ext𝑧 2− (() +() ) 𝛹R (𝑥, 𝑦, 𝑧, 𝑡) + 𝛹R (𝑥, 𝑦, 𝑧, 𝑡)22− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡)2 + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2 )𝛹R (𝑥, 𝑦, 𝑧, 𝑡))
𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)= 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
1+ 𝛥𝑡 (−𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)2
𝜕2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) 𝜕2 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)
+(
+
+
)
𝜕𝑥 2
𝜕𝑦 2
𝜕𝑧 2
1
1
𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝜕𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)2
2− 𝐵ext 𝑧+ 𝐵ext 𝑥
𝜕𝑥𝜕𝑧
𝐵ext 𝑥 2
𝐵ext 𝑧 2− (() +() ) 𝛹I (𝑥, 𝑦, 𝑧, 𝑡) + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2
2
2
1− (𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) + 𝛹I (𝑥, 𝑦, 𝑧, 𝑡)2 ) 𝛹I (𝑥, 𝑦, 𝑧, 𝑡))221(2.45)

𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
1= 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2+ 𝛥𝑡 (𝑉(𝑥, 𝑦, 𝑧, 𝑡)𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
11
𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡) 𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2
2
+(
+
𝜕𝑥 2
𝜕𝑦 2
1
𝜕2 𝛹R (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
2+) + 𝐵ext 𝑧
𝜕𝑧 2
𝜕𝑥(2.42)− 𝐵ext 𝑥𝜕𝛹I (𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)
𝜕𝑧

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

1
+ 𝛥𝑡))
222

𝑉(𝑥, 𝑦, 𝑧, 𝑡 + 𝛥𝑡)= (𝑉(𝑥 + 𝛥𝑙, 𝑦, 𝑧, 𝑡) + 𝑉(𝑥 − 𝛥𝑙, 𝑦, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦 + 𝛥𝑙, 𝑧, 𝑡)+ 𝑉(𝑥, 𝑦 − 𝛥𝑙, 𝑧, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 + 𝛥𝑙, 𝑡) + 𝑉(𝑥, 𝑦, 𝑧 − 𝛥𝑙, 𝑡)−

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

+

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

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

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

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

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

(2.46)

境界条件は
𝛹(𝑥, 𝑦, 𝑧, 𝑡 = 0) = cos𝜃 + i sin𝜃(2.47)𝑉(𝑥, 𝑦, 𝑧, 𝑡 = 0) = −𝐽𝑦/𝜎(2.48)

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

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

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

(2.51)

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

(2.52)

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

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

境界条件は23

𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡) 𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡)
=
=0
𝜕𝑥
𝜕𝑧

(2.55)

𝜕𝑉(𝑥, 𝑦, 𝑧, 𝑡)= −𝐽𝑦/𝜎𝜕𝑦(2.56)となる。

2.2 計算内容
本研究ではエネルギーの観点から計算の安定性を確認し、オーダーパラメータについてプロットすることで磁束線について評価を行う。エネルギーについてはポテンシャルエネ
ルギー𝐹P と運動エネルギー𝐹K について確認を行い、それぞれ
1
𝐹P = −|𝛹|2 + |𝛹|42

(2.57)𝐹K = |(−i∇ − 𝑨)𝛹|2(2.58)とする。また、|𝛹|2 < 0.1において磁束線が侵入していると判断するものとする。

2.2.1 横磁界における 2 次元 TDGL 方程式
初めに、動作確認を含めて既知の問題である横磁界中における超伝導体内の磁束線の動きについて、2 次元の TDGL 方程式を用いて解く。この時に仮定される状況を図 2.2 に示す。
電流を超伝導体の長手方向に印加し、超伝導体の面に対し垂直に磁界を印加する。観測方
向としては磁界と同じ方向を見るものとする。計算においては、ベクトルポテンシャルを電流と平行に印加することによって横磁界を再現するものとする。

2.2.2 横磁界における 3 次元 TDGL 方程式
次に、2.2.1 で計算した 2 次元横磁界についての計算を 3 次元に拡張する。この時に仮定される状況も図 2.2 のようになる。

2.2.3 縦磁界における 3 次元 TDGL 方程式
最後に縦磁界中にある超伝導体について 3 次元の TDGL 方程式を解く。この時に仮定される状況を図 2.3 に示す。電流と磁界を長手方向に平行に印加することで、縦磁界を再現す
るものとする。

24

𝐵𝐽

図 2.2:2 次元の模式図

𝐽, 𝐵(longitudinal)𝐵(transverse)𝑦𝑧𝑥

図 2.3:3 次元の模式図

25

第3章 計算結果
3.1 横磁界における 2 次元 TDGL 方程式
計算条件を表 3.1 に示す。また、今回の計算におけるエネルギーの時間依存を図 3.1 に、磁束線の動きの計算結果の一部を図 3.2 及び図 3.3 に示す。図 3.2 及び図 3.3 は左側が位相𝜃、右側が超伝導電子密度|𝛹|2を示している。
今回表 3.1 に示すようなパラメータを利用した理由としては、3 次元でだが以前に別の方法でこのような計算を行っていたため、似たような条件を利用した。空間離散幅と時間離散幅に関しては何度か検証した後、安定する中である程度計算速度を保てる条件としてこのパラメータを選択した。
エネルギーについては、初めに𝐸Kは高い値を示していたが、すぐに𝐸P とともに 0 付近に落ち着いた。その後、𝐸Kは55付近に、𝐸P は－100付近で一度落ち着いた後、𝑡 = 50あたりで−110まで下がり、
𝐸Kは振幅 2 で、𝐸P は振幅 4 で周期 8.7 の振動を始め、計算終了となる𝑡 = 300まで振動し続けた。エネルギーの総和に関しては、𝐸Kが安定的だったのもあり、𝐸P によって支配されているかのような挙動を取った。エネルギーが安定した後は、−55を基準に振幅5 で振動していた。
オーダーパラメータについては、全体において初期状態が|𝛹|2 = 1となるようにしていたにも関わらず、すぐに低い値となった。この時は𝜃も整っておらず、状態は観察できなかった。その後、𝜃が整うとともに|𝛹|2 も整いはじめ、𝑡 = 4には磁束線が左から右へと流れていく様子が確認できるようになった。さらに𝑡 = 12には三角格子を組んでいる様子も確認できた。その後は計算終了となる𝑡 = 300まで、磁束線が三角格子を組みつつ、左から右へと流れ続ける様子が確認できた。

26

表 3.1:2 次元横磁界におけるパラメータ
電流

0.2

磁場

0.4

試料サイズ

10×30

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1

100

energy

EK

0
EK+EP

EP

−100

0

100

200
time

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

27

300

(a)
𝜃

(b)

0

π

2π

0

0.5

1

|𝛹|2

図 3.2:2 次元横磁界における計算結果 (t=(a)4,(b)12)。右が𝑥軸正、上が𝑦軸正、電流は上方向、磁界は手前方向。

28

𝜃

(a)

(b)

(c)

(d)

0

π

2π

0

0.5

1

|𝛹|2

図 3.3:2 次元横磁界における計算結果 (t=(a)150,(b)160,(c)170,(d)180)。右が𝑥軸正、上が𝑦軸正、電流は上方向、磁界は手前方向。

29

3.2 横磁界における 3 次元 TDGL 方程式
計算条件を表 3.2 に示すパラメータでの計算を示す。また、今回の計算におけるエネルギ
ーの時間依存性を図 3.4 に、磁束線の動きの計算結果の一部を図 3.5 および図 3.6 に示す。
図 3.5 及び図 3.6 は|𝛹|2 < 0.1において𝛹の𝜃に準じた色をプロットしている。
今回表 3.2 に示すようなパラメータを利用した理由としては、3.1 節で計算したものを 3次元に拡張し、かつ同様な計算結果を得られるかどうか確認するためである。
エネルギーについては、初めに𝐸Kは高い値を示していたが、すぐに𝐸P とともに 0 付近に落ち着いた。その後、𝑡 = 23までにかけて𝐸Kは 500 付近に、𝐸P は−1100付近に変化した。その後は𝐸Kは振幅 20、𝐸P は振幅 50 の周波数 12 で振動した状態が計算終了まで続いた。
オーダーパラメータについては、全体において初期状態が|𝛹|2 = 1となるようにしていたにも関わらず、すぐに低い値となった。この時は𝜃も整っておらず、状態は観察できなかった。その後、𝑡 = 8からは磁束線の様子が観察できるようになった。𝑡 = 20には磁束線がまっすぐになり、三角格子を組んでいる様子が確認できるようになった。その後は計算終了時まで、𝑥軸正の方向に磁束線が流れていく様子が確認できた。

30

表 3.2:3 次元横磁界におけるパラメータ
電流

0.2

磁場

0.4

試料サイズ

10×30×10

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1

1000

energy

EK
0
EK+EP

−1000

0

EP

100

200
time

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

31

300

(a)
𝜃

(b)

0

π

2π

図 3.5:3 次元横磁界における計算結果 (𝑡 = (a)8, (b)20)。赤線が𝑥軸、青線が𝑦軸、緑線
が𝑧軸、電流は𝑦軸方向、磁界は𝑧軸方向。

32

𝜃

0

(a)

(b)

(c)

(d)
π

2π

図 3.6:3 次元横磁界における計算結果 (𝑡 = (a)150, (b)160, (c)170, (d)180)。赤線が𝑥軸、青線が𝑦軸、緑線が𝑧軸、電流は𝑦軸方向、磁界は𝑧軸方向。
33

3.3 縦磁界における 3 次元 TDGL 方程式
磁束線が動的である計算の一例として、表 3.3 に示すパラメータでの計算を示す。また、今回の計算におけるエネルギーの時間依存性を図 3.7 に、磁束線の動きの計算結果の一部を図 3.8 および図 3.9 に示す。図 3.8 及び図 3.9 は|𝛹|2 < 0.1において𝛹の𝜃に準じた色をプロットしている。
今回表 3.3 に示すようなパラメータを利用した理由としては、以前行った計算の再現が行うことができるかどうか確認するためである。空間離散幅及び時間離散幅は以前の計算ではプログラム側で自動調整していたが、今回は自ら設定を行った。
エネルギーについては、初めに𝐸Kは高い値を示していたが、すぐに𝐸P とともに 0 付近に落ち着いた。その後、20 < 𝑡 < 120において激しく振動しながら𝐸Kは600付近に、𝐸P は－1200
付近に向かった。変化しきった後に一度エネルギーの増加・減少が起き、𝑡 = 150以降は計算終了までエネルギーは安定していた。
オーダーパラメータについては、全体において初期状態が|𝛹|2 = 1となるようにしていたにも関わらず、すぐに低い値となった。この時は𝜃も整っておらず、状態は観察できなかった。その後、𝑡 = 11からは磁束線の様子が観察できるようになった。13 < 𝑡 < 33にかけて
磁束線が 1 本侵出した。46 < 𝑡 < 80では超伝導体の下から 1/3 の部分で|𝛹|2が大きく上下していた。t = 86, 105,110,116, 127では超伝導体の下から 1/3 のあたりで、t = 97では超伝導体の上から 1/3 のあたりで磁束線の組み換えが発生した。その後は 3 本の磁束線でらせん構造を取り、位相が回転しながらきつくなっていく様子が確認できた。最後に𝑡 = 280において、超伝導体の上から 1/3 付近で磁束線の組み換えが発生した。

34

表 3.3:磁束線が動的な計算の 3 次元縦磁界におけるパラメータ
電流

0.2

磁場

0.4

試料サイズ

10×30×10

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1

1000

energy

EK
0
EK+EP

−1000

0

EP

100

200

300

time
図 3.7:磁束線が動的な計算の 3 次元縦磁界におけるエネルギーの時間依存性

35

𝜃
図

0

(a)

(b)

(c)

(d)
π

2π

3.8: 磁 束 線 が 動 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)11, (b)33, (c)46, (d)50)。赤線が𝑥軸、青線が𝑦軸、緑線が𝑧軸、電流及び磁界は𝑦軸方向。
36

𝜃
図

(a)

(b)

(c)

(d)

0

π

2π

3.9: 磁 束 線 が 動 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)150, (b)250, (c)270, (d)300)。赤線が𝑥軸、青線が𝑦軸、緑線が𝑧軸、電流及び磁界は𝑦軸方向。
37次に、磁束線が静的である計算の一例として、表 3.4 に示すパラメータでの計算を示す。
また、今回の計算におけるエネルギーの時間依存性を図 3.10 に、磁束線の動きの計算結果の一部を図 3.11 及び図 3.12 に示す。図 3.11 及び図 3.12 は|𝛹|2 < 0.1において𝛹の𝜃に準じた色をプロットしている。
今回表 3.4 に示すようなパラメータを利用した理由としては、表 3.3 に比べて少ない電流を流すと、磁束線の運動が安定しやすいのではないかと考えたため、同様なパラメータから電流量を半分にした。
エネルギーについては、初めに𝐸Kは高い値を示していたが、すぐに𝐸P とともに 0 付近に落ち着いた。25 < 𝑡 < 100において、𝐸Kは 480 付近で、𝐸P は-1230 付近で僅かに不安定となるが、その後は一定となった。
オーダーパラメータについては、全体において初期状態が|𝛹|2 = 1となるようにしていたにも関わらず、すぐに低い値となった。この時は𝜃も整っておらず、状態は観察できなかった。その後、𝑡 = 8からは磁束線の様子が観察できるようになった。𝑡 = 23において超伝導体内には 4 本の磁束線が存在し、そのうちの 1 本は側面から外側に伸びていた。その後、その磁束線は時計回りに運動し、𝑡 = 50において運動方向にあった磁束線との組み換えが発生した。組み換え後の磁束線の 1 本は𝑡 = 97において超伝導体の外側に移動し、残った 3 本の磁束線が超伝導体内でらせん構造を取った。その後は計算終了まで位相が回転するのみ
で、大きな変化はなかった。

38

表 3.4:磁束線が静的な計算の 3 次元縦磁界におけるパラメータ
電流

0.1

磁場

0.4

試料サイズ

10×30×10

空間離散幅

0.2

時間離散幅

0.0005常伝導導電率1

1000

energy

EK

0

EK+EP
−1000

0

EP

100

200

300

time
図 3.10:磁束線が静的な計算の 3 次元縦磁界におけるエネルギーの時間依存性

39

𝜃
図

(a)

(b)

(c)

(d)
π

0

2π

3.11: 磁 束 線 が 静 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)8, (b)23, (c)50, (d)97)。赤線が𝑥軸、青線が𝑦軸、緑線が𝑧軸、電流及び磁界は𝑦軸方向。
40

𝜃
図

(a)

(b)

(c)

(d)
π

0

2π

3.12: 磁 束 線 が 静 的 な 計 算 の 3 次 元 縦 磁 界 に お け る 計 算 結 果(𝑡 = (a)150, (b)160, (c)170, (d)180)。赤線が𝑥軸、青線が𝑦軸、緑線が𝑧軸、電流及び磁界は𝑦軸方向。
41

第4章 考察4.1 横磁界における 2 次元 TDGL 方程式
動作確認のために、既知の問題として 2 次元の TDGL 方程式を解いてみたが、知られている通り磁束線がローレンツ力にしたがって流れていく様子が確認できた。また、磁束線が流れていく際に三角格子を組んでいる様子も確認できた。図 3.3(d)の|𝛹|2のプロットについてフーリエ解析を行った結果、図 4.1 のような結果が得られた。さらに図 4.2(a)について
フーリエ解析を行った結果、図 4.2(b)のような結果が得られた。縦横比が異なるものの、図4.2(b)を 1 周期とした形状が図 4.1 に確認できるため、三角格子を組んでいると考えられる。
エネルギーに関しては、磁束線の侵入・侵出に伴うエネルギーの増加・減少はあるものの、全体的におおよそ保存されていることから、安定した計算を行うことができたと考えられる。位相の特異点と|𝛹|2 の最も低い部分、つまり磁束線の中心と一致していることがわかる。
さらに今回の計算において、すべての磁束線の中心に対し位相がちょうど2π変わっていることが分かる。

𝑦

波長𝑥

図 4.1:図 3.3(b)のフーリエ解析
42

(a)

(b)

波長
𝑦

𝑥

図 4.2:(a)三角格子サンプルと(b)そのフーリエ解析

4.2 横磁界における 3 次元 TDGL 方程式
TDGL 方程式を 3 次元に拡張した場合の動作確認として、既知の問題である横磁界の場合について解いた。その結果、2 次元の場合と同様にローレンツ力に従って磁束線が動く様子が観察された。エネルギーについて注目すると、2 次元に比べ𝑧軸方向の大きさである 10 倍
大きくなっていることが分かる。エネルギーの辞官依存性についてのグラフを図 4.3 に示す。
図 4.4 は 3.1 節と計算結果は異なるが、同様の条件で計算を行い、|𝛹|2についてのプロット
を|𝛹|2 > 0.1において白、|𝛹|2 < 0.1において黒でプロットを行った図である。この図を参考に磁束線のサイズに注目すると、超伝導体内に存在する磁束線の大きさは 2 次元の場合は
平均 17.5、3 次元の場合は平均 16.2 とおおよそ同じ大きさとなった。さらに磁束線の配置も非常に近いものとなっている。これらの点により、おおよそ同様の計算ができていると考えられる。

43

図 4.3:エネルギーの時間依存性の比較

図 4.4:2 次元横磁界における計算結果

44

4.3 縦磁界における 3 次元 TDGL 方程式
𝐽 = 0.20の計算のエネルギーに関しては、初期状態の影響が長く残り、計算が安定するまで時間がかかった。今回の場合はおおよそ𝑡 > 150において計算が安定化したことがわかる。
次に磁束線の動きについては、超伝導体内に磁束線が 3 本らせん構造をとり、らせん構造が緩和されつつ位相が回転している様子が確認できた。直線になった後は上部で 3 本の磁束線が 1 点で結合・分離することで、磁束線が再構成される様子が観測できた。この一連の運動は、方位角成分のみ侵入するという磁束カッティングとも、超伝導体内で磁束線がフローするという force-free モデルとも異なっている。今回の計算においては、電流を流した際の自己磁界を考慮していないが、その場合には磁束カッティングモデルの場合は直線のみの磁束線が存在することが期待され、force-free モデルの場合には磁束線の組み換えが起きないことが期待されるはずである。また、他の計算方法でも今回の計算結果と同様であることが報告されている。
𝐽 = 0.10の計算においては、𝐽 = 0.20の計算に比べて非常に短い時間で計算が安定した。
これは初期値が乱数となっているため、影響の残り具合に大きく差が出たと考えられる。
|𝛹|2に注目すると、計算安定後において三角格子を組んだままらせん構造を組んでいることがわかる。これは、三角格子は磁束線が安定して分布する配置として知られており、らせん構造となっているのは𝐽 = 0.20と同様な回転運動が働いているためと考えられる。
𝐽 = 0.10ではその回転運動の働きが小さく、超伝導体が弱い状態になる直方体の角周辺で捕らえられるため、過剰にねじれることはなく安定して存在している。

45

第5章 まとめ
5.1 結言本研究では、
縦磁界中にある超伝導体内の磁束線の挙動を調査するために Time-Dependent
Ginzburg-Landau 方程式(TDGL 方程式)を用いた計算を行った。
TDGL 方程式をそのまま解くことは困難なため、磁場侵入長より短い超伝導体を仮定する
ことと、パラメータを規格化するという 2 つの近似を行った。この簡易化した TDGL 方程式を実部と虚部について分解し、交互に解くものとした。
初めに既知の問題である 2 次元の TDGL 方程式で横磁界中にある超伝導体について解く
ことで、作成したプログラムの動作確認を行った。その結果、従来の計算通り、超伝導体内にある磁束線がローレンツ力にしたがって動く様子が確認できた。この結果より、TDGL
方程式が正確に解くことに成功していると判断した。
次に、2 次元の TDGL 方程式を 3 次元に拡張を行い、再度横磁界中にある超伝導体について解いた。こちらでも 2 次元の場合と同様に、磁束線がローレンツ力にしたがって動く様子が確認できたため、3 次元に拡張した TDGL 方程式を解くことに成功していると判断した。
最後に縦磁界中にある超伝導体について 2 パターンの電流について解いた。一方の電流の場合では磁束線が回転運動しながららせん構造が変化し、磁束線の組み換えが確認された。もう一方の電流では、磁束線が完全に静止し、位相のみ回転している様子が確認された。

5.2 今後の課題
本研究では、磁場は 0.4 のみ、電流は 0.1 と 0.2 のみについてしか計測できておらず、具体的にどのような条件において磁束線が動くのか、または静止するのかということについて調査できていない。そのため、様々なパターンについて調査を行い、電流と磁場についての相図を作成する必要がある。
さらに今回は直方体の超伝導体についてのみの計算であり、直方体だと形状効果によって角の部分で超伝導状態が弱くなる傾向があるため、ここで磁束線が捉えられることが多い。そこで、形状効果を打ち消すことができる円柱状の超伝導体を設定し、これについても調査を行う必要がある。

46

謝辞
本研究を進めるにあたり、熱心にご指導いただきました国立大学法人 九州工業大学 大学院情報工学研究院 電子情報工学研究系 エレクトロニクス分野 教授 小田部 荘司先生
及び准教授 木内 勝先生に深く御礼申し上げます。小田部先生及び木内先生には、研究全
般の物理知識に関するご指導のみならず、技術者の心構え、社会人としての常識、安全意
識の必要性など、多岐にわたりご教授頂きました。さらに、小田部先生には共同研究を行う際に、先方との連絡を取って頂く、実際に会って打ち合わせする際の日程調整など、共同研究において御助力頂きました。心から感謝申し上げます。
独立行政法人 産業技術総合研究所 エネルギー技術研究部門 超伝導技術グループ 主任
研究員 馬渡 康徳様には、TDGL 方程式の基礎から簡易化まで、丁寧にご指導いただき、多大なる御助力・御助言を頂きましたことを心から感謝申し上げます。
国立有明工業高等専門学校 松野 哲也先生には、TDGL 方程式を解くプログラムを作成する際に、言語の紹介やプログラムの組み方など、多大なる御助力・御助言を頂きましたことを心から感謝申し上げます。
研究活動において、数多くの御助言・御助力頂きました、九州工業大学 大学院 情報工学研究院 電子情報工学研究系 エレクトロニクス分野の先生方に厚く御礼申し上げます。
最後に、公私共々お世話になりました小田部・木内研究室の皆様に深く感謝申し上げま
す。

47

参考文献
[1] T. Matsushita: Flux Pinning in Superconductor, second ed., Springer, 2014[2] Yu. F. Bychkov et. al.: JETP Lett. 9 (1969) 404[3] G. W.Cullen,R. L. Novak: Appl. Phys. Lett.4 (1964) 147[4] J. R. Cave et. al.: J. Phys. 39 (1978) C6-614[5] V. S. Vyatkin et. al.: IEEE Trans. on Appl. Supercond. 25 (2015) 6606207[6] J. R. Clem: J. Low Temp. Phys. 38 (1980) 353[7] E. H. Brandt: J. Low Temp. Phys. 39 (1980) 41[8] T. Matsushita, J. Phys. Soc. Jpn. 54 (1985) 1054
[9] V. L. Ginzburg, L. D. Landau, Zh. Eksp: Teor. Fiz. 20 (1950) 1064
[10] D. Y. Vodolazov, Phys. Rev. B 88 (2013) 014525[11] J. F. Blackburn et. al: Philosophical Magazine 80 (2000) 1455
[12] R. Barrett et. al.:Templates for the Solution of Linear Systems, second ed., SIAM, 1994
[13] J. F. Epperson: An Introduction to Numerical Methods and Analysis, second ed., WILEY, 1989
[14] Yoshida, H.: Phys. Lett. A 150 (1990) 262[15] A. Schmid: Phys. Condens. Mat. 5 (1966) 302[16] M. Tinkham: Introduction to Superconductivity, second ed., McGraw-Hill,1996[17] e.g. R. S. Thompson and C.-R. Hu: Phys. Rev. Lett. 27 (1971) 1352
[18] N. Kopnin: Theory of Nonequilibrium Superconductivity, Clarendon Press, Oxford, 2001[19] R. Kato, Y. Enomoto, S. Maekawa: Phys. Rev. B 44 (1991) 6916
[20] V. Vlasko-Vlasov et. al: Phys. Rev. B 91 (2015) 014516[21] J. R. Clem: Phys. Rev. B 26 (1982) 2463.
[22]松本 圭司、松本 要:平成 23 年度日本金属学会九州支部学術講演会(2011) B5[23]一野ら:第 75 回応用物理学会秋季学術講演会(2015) 17p-A21-15

48