1.1はじめに
超伝導とは、ある物質を非常に低い温度まで冷却すると電気抵抗がゼロになる現象である。1911年、オランダのLeiden大学の物理学者であるHeike Kamerlingh Onnesによって水銀が4.2Kの低温で電気抵抗がゼロになる現象が発見された。この物質のこの状態を超伝導状態と呼ばれ、今後、超電導物理学の研究が始まった。 1933年、MeissnerとOchsenfeldが超伝導体の磁気特性を研究していた時、超伝導体に外部磁界を加えても超伝導内部の磁束密度はゼロのまま変化しないという特性が発見された。この特性をマイスナー効果(完全反磁性)という。1935にはF. LondonとH. London兄弟によって、超伝導の2つの基本的な特性である電気ゼロ抵抗とマイスナー効果をよく説明したロンドン理論が提案された。 1950年にはV.L. GinzburgとL.D. Landauによって新しい現象論的理論の説明も提案された。これがGinzberg-Landau理論である。 G-L理論に導入されたバラメーターにとって、超伝導体は第１種超伝導体（）と第２種超伝導体（）に分類され、第２種超伝導体の電磁特性についても説明された。1957にはBardeen、Cooper、Schrifferらの超伝導の微視的理論（略してBCS理論）により超伝導発現機構が説明された。この理論は、物性物理学や物理学全体の発展に大きな影響を与えるだけでなく、超伝導研究を新たな段階へと導かれる。ついに 1986年、ドイツの科学者Johannes Georg BednorzとKarl Alex MuellerがLa-Ba-Cu-O高温酸化物超伝導体を作成し、その臨界温度Tcは35Kに達した、当時の多くの研究者は超伝導転移温度が30Kを超えることはできないという限界を超えた。この大きな発見は、超電導の新しい研究分野である高温酸化物超伝導体を切り開くと同時に、世界中の研究者による高温酸化物超電導材料の研究の急増をもたらした。
1.2 第一種超伝導体と第二種超伝導体
超伝導体は電気抵抗ゼロ、マイスナー効果の2つの性質を持っている。通常はマイスナー効果によって超伝導内部の磁束密度はゼロに保たれているが、外部磁場H。を高くしていき、臨界磁場H,を超えると、マイスナー効果が失われ常伝導状態となる。超伝導体は超伝導状態から常伝導状態へと遷移する状態の違いによって第一種超伝導体、第二種超伝導体の2種類に分類できる。第一種超伝導体も第二種超伝導体もHeが低い領域ではマイスナー効果を示す。第一種超伝導体はH以上のH 。をかけるとマイスナー効果を完全に失い、常伝導状態へと遷移する。対して第二種超伝導体は下部臨界磁界H c1以上のH eをかけたとき、超伝導体が破壊され磁束が一部侵入する。しかし、超伝導状態は完全には破壊されず、超伝導状態と常伝導状態が混在した状態となる。超伝導状態と常伝導状態が混在している、H 。がH ←以上で上部臨界磁場Hc2未満における状態を混合状態と呼ぶ。Heを高くすると常伝導状態の範囲が増え、入り込む磁束線の数が増える。H。がH 2以上となると超伝導状態は完全に破壊され、常伝導状態に遷移する。このように、第一種超伝導体と第二種超伝導体は混合状態の有無で分類されている。
1.3
(Time-Dependent) Ginzburg-Landau 方程式
ここでは、超伝導体の数値解析に用いた時間依存 Ginzburg-Landau(Time-Dependent Ginzburg-Landau)方程式について説明する。
1.3.1 Ginzburg-Landau(GL)方程式
Ginzburg-Landau(G-L)理論は、1950年にV. L. GinzburgとL. D. Landauによって提唱された超伝導現象を説明する現象論である[1]。G-L理論は磁界が存在するときの超伝導が共存する場合の相転移を取り扱ったもので、特に第２種超伝導体の磁気特性を良く記述することが知られている。以下では、G-L理論を成り立たせる仮定とそれによって求められる方程式を説明する。[2] G-L理論では、まず超伝導状態の秩序の程度を表す量として複素数であるオーダーバラメーターΨ = |Ψ|exp(iφ)を定義する。そしてΨは、次の関係を満たすと仮定する。 |Ψ|² ∞ n − (1.1)ここで、n sは超伝導電子密度である。超伝導状態の自由エネルギーE sはn sに依存しているから(1.1)式より、|Ψ|²の関数である。よって、E 5は|Ψ|²のべき展開で表すことが出来る。
ここで、Enは常伝導状態の自由エネルギー、αおよびβはそれぞれべき展開した際の１次と２次の展開係数であり、T < Tcではα < 0,β > 0である。２次の項までの展開であるのは|Ψ|²が転移点近傍(T ≒ Tc)において十分に小さいからである。次に、磁界の存在がΨの空間的変化に寄与することを考慮し、(1.2)式に磁界のエネルギー密度と運動エネルギー密度を加算する。
ここで、μ oは真空中の透磁率、Aはベクトルポテンシャル、m*は超伝導電子の質量、あはディラック定数、e*は超伝導電子の電荷量、iは虚数単位である。 ΨとAは(1.3)式のEが最小なるように決定される。よってΨの共役複素数Ψ*とAについて

5

うになる。
(1.13)式、(1.14)式、(1.15)式から得られる∂Es  ∂Esをすべて(1.12)式に代入すると以下のようになる。成分をまとめると次のようになる。
ションを持つ項すなわち(1.3)式でここで、E sのうち微分されたA xを含むのは、ロトいうところの右辺第４項だけであるので、それ以外の項は０となり、(1.19)式は以下のようになる。
∂A x ∂A x ∂A xによる偏微分を行うと次のようになる。
 ∂x ' ∂y ' ∂z
8
rot(rot A) | x
Po y成分、z成分も(1.21)式と同様に変形できる。
9 E ДЕ ç ∂Esをすべて(1.18)式に代P (1.21)式、(1.22)式、(1.23)式から得られる∇・ AVA ∂∇ A x ∂VAz入すると以下のようになる。

1.3.2 Time-Dependent G-L (TDGL) 方程式
Time-Dependent G-L(TDGL)方程式は、G-L方程式に時間依存性を付与したものである。 (1.4)及び(1.5)式に対して時間発展する場合を考えると、以下の2式となる。
ここで、γとνはそれぞれΨとΛの時定数である。さらにゲージ変換を与えると次のように9なる。
ここで、Vはスカラーポテンシャルである。(1.28)、(1.29)式を(1.26)、(1.27)式それぞれに代入すると最終的に次のようになる。

1.3.3
TDGL方程式の簡易化
本研究ではTDGL方程式を数値解析するが、(1.30)、(1.31)式をそのまま解くのは困難であるため、細線近似と規格化による2つの簡易化を行う。本研究の数値解析では、Aは時間に対して一定に与えるから(1.31)式の左辺第1項の時間偏微分が0となる。次に(1.31)、(1.32)に対して規格化を行う。超伝導体のコヒーレンス長ζと磁界侵入長λは熱力学的臨界磁界をH cとして以下の2式のように表す。
そして、以下の式に記す規格化を行う。

10

(1.30)式の規格化を行う。(1.30)式の左辺第1項は、次のようになる。
(1.30)式の左辺第2項は、次のようになる。
(1.30)式の左辺第3項は、次のようになる。 11
(1.30)式の左辺第4項は、次のようになる。 (1.39)式から(1.42)式をまとめると、(1.30)式は最終的に次のようになる。 (1.31)式も同様に簡易化する。(1.31)式左辺第一項は、Aが一定として次のようになる。 (1.31)式左辺第2項は次のようになる。
(1.31)式左辺第3項は次のようになる。
(1.31)式左辺第4項は次のようになる。
(1.45)式から(1.48)式をまとめると、(1.31)式は最終的に次のようになる。
ここで、Jは超伝導体に流れる電流であり、(1.52)式は電流の発散を表す式である。(1.31)式の両辺で▼との内積を取ると、(1.51)式、(1.52)式から、次のような式が得られる。 (1.49)
いこで、
よび(1.53)式を解くことによって超伝導体の解析を行う。

1.4 アフィン数値積分法(Affine Integrator; AFI)
アフィン変換数値積分法(AFI)[3][4][5]とは、近年提案された陽的数値積分法である。拡散方程式のような放物型偏微分方程式や、ゲージ場存在下のSchrödinger方程式やTime- Dependent Ginzburg-Landau(TDGL)方程式の数値積分のために考案された。数値的安定性に関して、AFIはラプラシアン項に関して無条件安定であり、方程式を構成する項全体に関して高い数値的安定性を有する。エネルギー保存に関しては、AFIは線形保存系(Schrödinger方程式やTDGL方程式)において全エネルギーが厳密に保存される性質を持つ。 AFIは空間に関して離散化するための格子をチェッカーボード状に分解することから導かれるアフィン変換対で構成される。数学的に単純なスキームであり、数値的安定性やエネルギー保存に関する理論解析が容易である。また必要とされる記憶領域は最小であり、付加的なワーキングメモリは必要とされないこともAFIの特徴である。
1.4.1  AFIのイメージ
図1-1の左側に空間に関する離散化のためのチェッカーボード格子を示す。基礎方程式はこのチェッカーボード格子によって、まずは空間に関して離散化される。右側にチェッカーボード格子に基づいて構成された高次元常微分方程式の概形を示す。黒丸格子点の集合をΩ集合、白丸格子点の集合をP集合と呼ぶことにする。 Ω集合およびP集合上で標本化された従属変数の時間発展は図1-1の右側に示された方程式のプロ
15
ック行列が支配する。[6]図1-2に常微分方程式を時間に関して離散化するようすを示す。このように、AFIスキームはＱ集合上で定義された従属変数の時間発展とＰ集合上で定義された従属変数の時間発展を交互に繰り返していくものである。図1-1図1-2時間に関する離散化のイメージ
1.4.2    1次元空間においてAFI法の時間発展方程式の導出
1次元空間でのAFIの数値積分スキームを導出する。まずは空間に関する離散化を行う。本研究で用いる基礎方程式はTDGL方程式である(1.55)式とする。
1次元空間では(1.56)式が(1.55)式で適用されており、ベクトルポテンシャルは固定されていると考える。
ベクトルポテンシャルはリンク変数wi,jによって次のように実装される。
ここでx i, x iはそれぞれ格子点i, jのx座標で、hは格子点間隔である。格子点iと格子点jが接続していない場合はw i j = 0と約束する。簡単のため、周期的境界条件を採用する。空間に関する離散化のために導入した格子およびリンク変数の状況を図1-3に示す。図1-3
(1.55)式を空間に関して離散化する。ψ i (t) = ψ(x i ,t),α i = α(x i ),β i = β(x i )とおくと(1.58)式で表される連立微分方程式が得られる。
奇数の格子点の組(図1-3の黒丸に対応)と偶数の格子点の組(図1-3の白丸に対応)の順に並べているが、図1-1の右側のイメージに従って表記している。この連立微分方程式は以下の(1.59)式で表せる。 d ψ 3 Y dt ψ 2
ここで、
とおく。 ∑ iはψ iの推定値である。理由は後述する。また、
とおく。また、
とすると、(1.59)式は(1.64)式で表せる。
時間推進演算子Åを用いると(1.64)式は(1.65)式で表される。また、 Åは(1.66)で定義される。
ただし、
である。また、Wij = Wjiである。次に、(1.65)式を時間に関して離散化することを考える。時間刻み幅をτ、(q,p) = (q(t),p(t))、(q',p') = (q(t + τ),p(t + τ))として時間発展方程式の(1.70)式が得られる。

18
LTS分解とは、exp x (A + B)の形のものを(1.72)式のように積に分解する方法の1つである。LTS分解でどのように分解できるかは後述する。また、Oは誤差オーダー部分である。
従って、指数関数演算子exp (ξ Âq)とexp (ξ Âp)の行列表現が得ると、数値積分スキーム(AFI)を構成できる。まず、ここで演算子Âǫをψi(i ∈ Q)に1回作用させる。
次に、演算子Âǫをゆ¡(i ∈ Q)に2回作用させる。
従って、演算子Âǫをψi(i ∈ Q)にm回(mは正の整数)作用させると
となる。
19

が得られる。ただし、
であるので、
ここで、Iは単位行列であり、
とおいた。 (1.80)式と(1.81)式はそれぞれ指数関数演算子exp(ξ Âρ )とexp(ξ Âρ )の行列表現を表しており、数値積分スキーム(AFI)が得られたことになる。次に、(1.60)式の非線形項βi|φi|²で、ψではなく推定値φ iが用いられる理由を説明する。推定値は(1.83)式から(1.86)式で計算される。
注意する点として、ψi(i ∈ Q)の価はψi(i ∈ P)で推定され、ψi(i ∈ P)の価はψi(i ∈ Q)で推定されることである。このことにより、指数関数演算子の行列表現構成が容易になる。例を挙げると、指数関数演算子exp(ξÂo)の行列表現を得るために、(1.75)式においてÂoをベクトルqの要素に何回も作用させた。その際、σi,i ∈ Qはベクトルqの要素に依存せず、ベクトルpの要素に依存することが重要である。つまり、時間推進演算子Âoにとってσi,i ∈ Qは定数となる。このことが行列表現の導出を可能としている。
1.4.3    2次元空間においてAFI法の時間発展方程式の導出
2次元空間でのAFI法の数値積分スキームを導出する。基礎方程式は1次元空間で用いたものと同様に(1.55)式を用いる。ただし、
である。ここで、格子点やリンク変数を指定するためにダブルインデックスを用いる。そうすることでダプルインデックスを格子点やリンク変数を空間的な位置に対応させることができる。図1-4にダプルインデックスの定義を示す。また、リンク変数はx成分のw iと、y成分のwζ」に分けて定義される。境界条件はx方向、y方向ともに周期的境界条件とする。 . X 2次元空間での周期的境界条件下の図1-4格子点とリンク変数のダブルインデックスでの定義まず空間に関して離散化すると(1.88)式が得られる。
とおく。また、 σi,j = 4 + Ui,jh²とおく。ここまで格子点とリンク変数をダブルインデックスで定義したが、図1-4において格子点を(1,1)、(1,2)、…、(4,3)、(4,4)と定義している点をシングルインデックス1、2、…、 15、16と定義したものを考えると、1次元空間で導出した数値積分スキームの(1.80)式と(1.81)式、時間発展方程式の(1.76)式がそのまま用いることができる。数値積分スキームは図1-4の黒丸に対して(1.80)式、白丸に対して(1.81)式が対応する。具体例として格子点(1,1)での時間発展方程式を(1.92)式に示す。なお、従属変数ψιιが時間発展した後のものをψ′11とする。(1.76)式を適用するが、使用されている格子点とリンク変数に関してはダプルインデックスで定義した表記を使用している。ほかの格子点についても同様に(1.76)式を用いて導出可能である。
(1.95)
これを一般化すると、2次元空間におけるAFI法での時間発展方程式の(1.97)式を得ることができる。

1.5 可逆現象
ここでは、超伝導体内の磁束線の可逆運動について説明する。[1]超伝導体内の不可逆現象はピン・ポテンシャル付近での磁束線の運動、その運動がピン・ポテンシャル内に限られ、臨界状態モデルによる記述に合わなくなってくることが予想される。ここで、ある領域の磁束線がビン・ポテンシャル内の平衡状態にあるとする。単位体積当たりの磁束線が受ける力は以下のようになる。

23

(1.102)
ここで、uは平衡位置からずらす変位量、Labuschバラメーターαιは磁束線が感じたピン・ポテンシャル。磁束線の変位方向をx軸方向として、磁束密度bの連続の式を時間に関して積分する。
ここで、Bは磁束密度の平衡値である。この変化により生じるLorentz力の変化分は以下のようになる。
Lorentz力とピン力との釣り合いにより(1.102)~(1.105)式からuを消去する。 F + FL = 0
したがって、超伝導体がx ≥ 0を占める場合には、磁束密度は以下のようになる。
ここで、はCampbellの交流磁場の侵入深さと呼ばれる量である。
臨界状態から増磁に転じた後の磁束分布の変化を図1-5(a)に示す。一方、臨界状態モデルでは図1-5(b)に示す。 (1.105) (a) (b)図1-5初期状態からの変位量uが小さい範囲ではでビン力密度は以下のようになる。 F = −αι u + J c Bしたがって、Lorentz力密度を得られる。 (1.109)
この二つの力の釣り合いの式から再び(1.107)式の解が得られる。ここで、ピン力密度-変位特性を以下のように表すことができる。
(1.112)
(1.103)式を用いて力の釣り合いの式からuを消去するとの非線型微分方程式が得られる(図1-6参照)。

1.6 粒子群最適化(Particle swarm optimization; PSO)
ここでは、粒子群最適化いついて説明する。[8]粒子群最適化法は1995年にJ. KennedyとR. Eberhartにより開発された、母集団型の確率的なメタヒューリスティクスの最適化手法である。PSOは鳥や魚などの群れを成して移動する生物の行動を模倣したアルゴリズムである。これらの群れを成す生物を粒子(Particle)としてモデル化し、粒子は最適化問題における解候補を表す。粒子群を多次元解空間に配置し、各粒子が他粒子との相互作用により探索を行うことで問題に対する近似解を求めることができる。粒子群におけるi番目の粒子(粒子i)は、速度ベクトルvi、位置ベクトルxi、自身の探索過程における最良の位置ベクトル(Personal best solution;  pbest)pbiを持つ。また、粒子は粒子群の探索過程における最良の位置ベクトル(Global best solution; gbest)gbを共有する。Dは最適化問題の次元数を表す。
(1.116)
(1.118) tステップ目の粒子におけるj次元目の速度と位置の更新式を表すことをできる。
ここで、wは慣性速度定数、c 1 、c 2はpbest、gbestそれぞれに対する加速度係数（Cognitive coefficient;  c 1 、Social coefficient;  c 2 )、 r 1 、 r 2は[0,1]の一様乱数を表す。粒子は図1-7および更新式に示すように自身の速度とpbest、gbestの影響を受けて探索を行い、最終的に全ての粒子がgbestに到達すると粒子は停止するため、探索が終了する。また、粒子モデルを図1-7に示す。
1.7研究の目的
理論上では、超伝導体の可逆運動は予測されたが、現実には量子化磁束線の可逆運動の観測が困難である。本研究では、AFI(Affine Integrator)数値積分法を用いて、時間依存Gizburg-Landau(TDGL)方程式を解き、量子化された磁束線を可視化し、数値計算を行い、理論的な磁束可逆運動を二次元で再現することを目的として行う。

2.1実装方法
実装するにあたって、統合開発環境でありプログラミング言語でもあるProcessingで実装を行う。Processingは、アーティストよるコンテンツ制作作業のために、詳細な設定を行う関数を排除しており、非プログラマーでも容易にプログラムを実装できる。本研究では、2次元超伝導体をシミュレーションし、超伝導体内の量子化磁束線を可視化する。
2.2    2次元超伝導体
ここでは、定義域は矩形のx ∈ [0,Lx] × y ∈ [0,Ly]であるとする。従属変数ψ(x,y,t)が従う基礎方程式は次のようになる。また、周期的境界条件を採用する。
ただし、格子点間隔をk、格子点の数をN、定義域の長さをLとする。定義域の長さLと関係は以下のようになる。
ここで、周期的境界条件として、以下のようにする。 x方向(左端と右端)は次のようになる。
y方向(上端と下端)は次のようになる。
初期条件は次のようにする。

2.3 規格化
ここではTDGL方程式の規格化について説明する。 [9]
2.3.1  空間に関する規格化
ここではTDGL方程式から導く。

2.3.2
磁場の規格化
マクスウェルの方程式による

2.3.3
ーダーパラメーターと侵入長に関する規格化
オ
ベクトルポテンシャルに関する方程式(2.10)の両辺に以下の式
29

2.3.4
可逆運動の非線型微分方程式の規格化
ただし、ここで規格化オーダーパラメーターψを5でスケール変換して、Ψを無次元化にする。 (2.17)式の右辺において、Im[]が無次元量で式全体も無次元なのでその係数も無次元量である。したがってここで侵入長λ [m]は式(2.10)に基づき次式で定義される。また、可逆運動の非線型微分方程式によるここでは、磁場密度b、x方向の変位量、超伝導体の幅2dを規格化する。規格化された可逆運動の非線型微分方程式は以下のようになる。
2.4  計算条件
2.4.1　臨界状態からの可逆運動
可逆現象をより顕著に観察し、且つ臨界状態モデルと比較するために、外部磁場B,を変えて、超伝導体を臨界状態から正弦波交流磁場に加える(図2-1参照)。図2-1
2.4.2 ゼロ磁場からの可逆運動
超伝導体は臨界状態以外には、顕著な可逆運動は確かめるため。外部磁場B2を変えて、超伝導体をゼロ磁場から正弦波交流磁場に加える(図2-2参照)。図2-1外部磁場の経時変化

3.1   臨界状態からの可逆運動のシミュレーション
まず、超伝導体に200個ピンを入る、規格化した直流磁場Bz = −4を加える。臨界状態になるため、待機時間t = 10000ステップを待つ。続いて、正比例に直流磁場Bz = 2まで増磁する。B z = 2のところに、待機時間t = 10000ステップを待つ。最後は、超伝導体にB z = 1.5 + 0.5 sin ωtの規格化した磁界を加えて、ωtがπ/2、5π/8、6π/8、7π/8、πのとき、y軸方向のB xを全部加算して、Bmaxに規格化する。磁束密度は侵入深さによって徐々に減少する。ωtの増大によって減磁する。図3-1参照。図3-2図3-3は、ωtが6π/8、7π/8、πの時の式(1.114)の理論式と数値計算による磁束密度の差の分布を比較して示す。理論計算は直線で表す。その結果、磁束線が顕著な可逆運動するという理論予測に近い。
3.2ゼロ磁場からの可逆運動のシミュレーション
3.2.1   ピンが200個の場合
超伝導体に 200 個ビンを入る、周期T = 5000ステップの規格化した正弦波交流磁界Bz = 2 sin ωtを加えて、ωtがπ/8、2π/8、3π/8、4π/8のとき、y軸方向のB xを全部加算して、B maxに規格化する。磁束密度は侵入深さによって徐々に減少する。ωtの増大によって増磁する。ゼロ磁場からの可逆運動ので、0からの磁束密度の差は磁束密度と同じ。図3-4参照。図3-5は、ωtがπ/8、2π/8、3π/8、4π/8の時の式(1.114)の理論式と数値計算による磁束密度の差の分布を比較して示す。理論計算は直線で表す。その結果、ωtが3π/8から、理論計算と一致していない。可逆運動がゼロ磁界でより顕著になっている。直流磁界が印加されると、可逆運動よりも不可逆運動が強調されている。また、x > 5のときは0ではなく、少し磁場が入れた原因は磁束クリープと考えている。
3.2.2 ピンが100個の場合
超伝導体に100個ピンを入る、周期T = 5000ステップの規格化した正弦波交流磁界B z = 2 sin ωtを加えて、ωtがπ/8、2π/8、3π/8、4π/8のとき、y軸方向のB xを全部加算して、B maxに規格化する。磁束密度は侵入深さによって徐々に減少する。ωtの増大によって増磁する。ゼロ磁場からの可逆運動ので、0からの磁束密度の差は磁束密度と同じ。図3-6参照。超伝導体内の磁束密度と侵入深さの関係図3-6図3-7は、ωtがπ/8、2π/8、3π/8、4π/8の時、ビンの数による超伝導体内の磁束密度と侵入深さの関係。ビンを100個入るの場合は直線、ビンを200個入るの場合は点。その結果、ωtが3π/8から、同じ侵入深さで、ビンを100個入るの場合の磁束密度が高い。ビンの数は多いほど可逆運動の磁束クリップの影響が小さくなる。
本研究では、TDGL方程式をとくことで、AFIを用いて量子化磁束線を二次元で可視化した。磁束線の可逆運動と理論計算による予測の一致性を確かめた。臨界状態からの可逆運動のシミュレーションを行った。臨界状態になるまでは、結構時間がかかるので、直流磁界B, = −4からB, = 2までの範囲で、非常に長い待機時間tが設定された。臨界状態になるから、交流磁界を加えて、磁束密度と理論計算に比較した。その結果、磁束密度は折り曲げる点がなく、ゆっくり減少することて、顕著な可逆運動するという理論予測に近い。ゼロ磁場からの可逆運動のシミュレーションを行った。超伝導体に交流磁界を加えることで、臨界状態ではないときの可逆運動を理論計算に比較した。その結果、不可逆的の臨界状態モデルに違い、可逆運動が特に顕著ですが、直流磁界が印加されると、不可逆運動がより際立っている。最後に、異なるピン留め密度をシミュレーションし、可逆運動に対する影響を検討しました。同様に、直流磁場がゼロの状態で、交流磁場が印加した。100本のビンがある超伝導体において、磁束クリープは200本のビンがある超伝導体よりも多いことが確認されました。
学会発表1.温謙, 上田天馬, 小田部荘司,松野哲也「AFI数値積分法を用いた2次元の超電導体内の量子化磁束線の可逆運動の研究」第105回低温工学・超電導学会研究発表会，東京都江戸川区タワーホール船堀,令和5年6月2. 温謙, 上田天馬, 小田部荘司, 松野哲也 ｢AFI数値積分法を用いた2次元の超伝導体内の量子化磁束線の可逆運動の解析」第 71 回応用物理学会春季学術講演会,東京都市大学,オンライン,令和６年３月