FDTD法によるTDGLシミュレーション	
リンク変数法によるシミュレーション，および	  オブジェクト指向プログラミングによる込み入った数式の実装	
有明高専 電子情報工学科

松野哲也	

応用物理学会２０１５秋（９月１３日（日）12:00-­‐13:30） 	  インフォーマルミーティング 「量子化磁束動力学シミュレーション研究会 」	
1	

FDTD法によるTDGLシミュレーション	●	  FDTD： 時間領域有限差分法	  （ 偏微分方程式	  →	  時間発展差分方程式 ）	  	  ●	  FEM	  ： 有限要素法	  （ 偏微分方程式	  →	  変分問題	  →	  線形方程式 ）	

2	

ix A

Time dependent GL equation

TDGL方程式の導出	
ensity is defined as
ϵGL = ϵK + ϵP + ϵB ,ϵKϵPϵB

1=|Dψ|2 , D ≡ ∇ − igA,2
η=(|ψ|2 − 1)2 ,21=|∇ × A|2 .
2iation:δ≡∂−∇·∂,

δ=∂−∇·∂

3	, i = x, y, z∂t+ ∇φ = −δA,e TDGL equations will be obtained as!
" us note that the equations (25) and (alar potential.
Let∂
δϵGL
τ+ igφ ψ = −,(25transformation:∂tδψ
TDGL方程式	  ∂AδϵGL ∂χ 時間依存ギンツブルグランダウ	  + ∇φ φ=−→− φ −,(26方程式	  A −→ A∂t+ ∇χ,,ψ−→
ψexp(igχ).
δA(	  Time	  dependent	  Ginzburg-­‐∂t)	
tential. Let us note that the equations (25)Landau	  andequa@on	  (26) areinvarian
Lrmation:equations is obtained as follows:!

"∂∂χ→A+φ −→− =,τ ∇χ, +igφ φ ψ∂t∂t∂A as follows:ions is obtained+ ∇φ ="∂t1+ igφ ψ =2

1 22ψD−→ψψ+exp(igχ).
(27η(1 − |ψ| )ψ,2#$gIm ψDψ − ∇ × ∇ × A.

(∇ − igA)2 ψ + η(1 − |ψ|2 )ψ,#$B+ ∇φBoundarycondition= gIm ψ(∇ − igA)ψ − ∇ × ∇ × A.

A(284	(29∂t+ ∇φ = −δA,

he scalar potential. Let us note that the equations (25) and (26)auge transformation:

ゲージ対称性	

equation
expressed
byvariables
link varian expressedby link∂χA −→ A + ∇χ, φ −→ φ −, ψ −→ ψ exp(igχ).
∂tTDGLequations
is obtained
asgaugefollows:is as follows:
equations
oftheφ-zero
the φ-zero! gauge" is as follows:∂1 22τ ∂ψ + igφ 1
ψ =2 D ψ + η(1 − |ψ|)ψ,21τ ∂t2 = D ψ2+2η(1 − |ψ| )ψ,=D∂tψ
#$
2 − |ψ| )ψ,∂A+ η(1− ∇ × ∇ × A.
2 ∂A ∂t + ∇φ =! gIm ψDψ"

! = gIm"ψDψ − ∇ × ∇ × A.
∂t ψDψ − ∇ × ∇ × A.
= gIm
dix B Boundary conditiony discretized version of the above eqs.(14) and (15)a rectangle
the two-dimensional
simulation
area. Thereare ldderversion
ofastheabove eqs.(14)and (15)
can ber parameter5	なぜ φゼロゲージ，なぜリンク変数法か？	
on expressedby link variablesof the φ-zero gauge is as follows:ψ
1FDTD法と相性が良い．	
22ﬂux	  ﬂow	  シミュレーションを	  =D ψ + η(1 − |ψ| 行うにはリンク変数法が必須．	
)ψ,t2!
"A
陽的数値積分を組める．	
= gIm
ψDψ − ∇ × ∇ × A.
tE	  一定は	  A	  の時間微分が一定	  であることを意味するから．	  
つまり，t	  →	  ∞	  で	  A	  は有界ではない．	

並列化が容易	

ed version 高速化	  of the above eqs.(14) and (15) can bei+1,ji−1,ji−1,ji,ji,j+1i,j−1i,j−16	リンク変数法	
Ny

格子点の上で	  オーダーパラメータを定義	

.

格子点と格子点の間で	  ベクトルポテンシャルを定義	

.
.

21

12. . .

Nx7	Appendix CEquations with the link v

yj+1 link variables as
Let us defineUxi,the
リンク変数：	i, j+1

Lz i, jUy i, ji, j

Ux i, jUy i+1, j

Uxi,j,k = exp(−ihgAi,j,kx ),Uyi,j,k = exp(−ihgAi,j,ky ),hxi+1, jUzi,j,k = exp(−ihgAi,j,k),zhwhere h is thelattice spacing. By using these link variables, w
discrete version of the equations (28) and (29).
The dynamics of the order parameters are described as follows

Ny + 1Ny.

∂ψ i,j,k1
i−1,j,k i−1,j,ki,j,k i+1,j,ki,j,k i,j+1τ
=(Uψ+Uψ+Uψ
xxy2∂t2hi,j,k−1 i,j,k−1+Uzi,j,k ψ i,j,k+1 + U zψ
− 6ψ i,j,k ) + η(1.
.

210012. . .

Nx

Nx + 1ψUx

8Uy8	n1Dx ψi,j=1 i+1,j−igA= Vi,jx ni,j(Vi,ji,jxi+1,jψi,jx ψ),
i+1,j1i,ji+1,ji+1,ji+1,ji,jψ i,j )ψ−→V(Vψ−V∂x∂x−→DV(Vψ−Vψ)x
Dx ψx −→x Vxx h (Vxxψ− Vxx ψ )xhh
1i,j1i,ji+1,ji+1,ji,j1
共変微分演算子（２次元）	
i,ji+1,ji+1,ji,jhis,=we can(Vobtainthediscretizedversionofthei+1,ji+1,ji,j=(VVψ−ψ)ψxx Vxx− ψψ ) − ψ )x=Vx h (Vhh1 i,ji+1,j11 i,j i,j1 i,j i+1,ji+1,j i+1,ji,j i,ji,ji,ji+1,ji,jDψ−→V(Vψ−V)=(Uψ−ψ),= x (Ux=ψx−x ψxxψψx),(U−ψ),
hh hh11i,j i+1,j i+1,ji,ji,j22 i,j 1i−1,ji−
1i,ji+1,ji+1,ji−1,ji−1,ji,i+1,ji+1,ji−1,ji−=(VVψ−ψ)
D
ψ−→V
(V+Vψ
−→DxxVψx −→(Vx Vxx ψ2x (Vxx+ Vxψ ψ + Vxx− 2Vψx2h hh2
h111 i,j i−1,ji−1,j1 i,j i+1,ji+1,ji,ji−1,ji−1,ji−1,ji,ji,ji+1,ji−1,j=(Uψ−ψ),=+U
ψ−2=(Uψ+Uψ−2ψ)
xx=
(Uψ+U
ψ
−xx
222xxhhhhi,j 12i+1,j i+1,ji−1,j iD
ψ−→V(Vψ+Vψplains
theabovediscretizations.
xxxxe
above
discretizations.
xplains the above discretizations.
h21 i,j i+1,ji−1,j i−1,j9	Uxi, j+1
リンク変数の積は経路積分に対応する．	

Uy i, jUy i+1, ji, jdefine a variable Li,jusing
linkvariablesasz
Ux i, jLi,jz≡=≃i,j i+1,j i,j+1 i,jUx Uy U x U yi,ji,j+1i,j
exp[−ihg(Ai+1,j−A−A+Ayyxx )]exp[−ih2 g(∇ × A)i,jz ],i,jz is the z-component of the rotation of the vector potenti10	) or (xi + h/2, yj + h/2). That is, the variable Li,jcorrespozα=1,2

1 i,j,k i−1,j,k i,j,k i,j,k−1i, j, k− 2 (Lz Lz
Lx Lz z−Lx1),
hi,j,k3次元では	!
∂Uzi,j,k i,j,k i,j,k+1
2 i,j,k 2= −ig(1 + κ |ψα | )Im[ψ α Uz ψα]∂tα=1,2
1 i,j,k i,j−1,k i,j,k i−1,j,k− 2 (Lx LxLy i, j, k Ly Ly i, j, k − 1),
hyi,j,ki,j,k xp” variables Li,j,k,L,L
are defined asxyz
Lz i, j, kLi,j,kx=

Li,j,k=yLi,j,k=zi,j,k i,j+1,k i,j,k+1 i,j,kUy Uz
UyUz ,i,j,k i,j,k+1 i+1,j,k i,j,kUz Ux
UzUx ,i,j,k i+1,j,k i,j+1,k i,j,kUx UyUx
Uy .

physical meaning of the loop variable is expponential of thepa
11	

i,ji,j+1i,j= exp[−ihg(Ai+1,j−A−A+Ayyxx )]2i,j≃exp[−ihg(∇×A)],z
２次元	Lz i-1, j

Lz i, js the z-component of theLz i, jrotation of the vector por (xi + h/2, yj + h/2). That is, the variable Li,jz corr
j-1i, j(xi + h/2, yj +i, jh/2). See
Lz i, Fig.3(a).
an obtain the expressions asi,j i,j−1
Lz Lzi,j i−1,jLz Lz

≃≃

3i,j
exp[−ih g(∇ × ∇ × A)x ],exp[−ih3 g(∇ × ∇ × A)i,jy ],(12) are defined at the sites (xi +h/2, yj ) and (xi , yje same positions for Uxi,j and Uyi,j , respectively.
12	L equation expressed by link variables

GL equations of the φ-zero gauge is as follows:

2TDGL equation expressed by link variables∂ψ1∂t2∂ψ ! 1 2"∂Aτ=D ψ + η(1 − |ψ|2 )ψ,= gImψDψ− ∇ × ∇ × A.
∂t
2
∂t!
"
∂A2 φ-zero gauge is 2as follows:
The set of TDGLτ equations= ofDtheψ + η(1 − |ψ| )ψ,∂t((14(= gIm ψDψ − ∇ × ∇ × A.

(15tially
discretized version of the above eqs.(14) and (15) can be obtainedSemi-­‐discre@zed	  equa@on	  (	  2次元	  )	Hence, the spatially discretized version of the above eqs.(14) and (15) can be obtained as

1 i,j i,j i+1,ji−1,j i−1,ji,j−1 i,j−1i,j i,j+1i,j∂ψ
1i−1,ji,j−1= τ 2 (Ux =ψ+U
ψ+U
ψ
+U
ψ
−4ψ
i,j i+1,ji−1,jyi,j i,j+1 yi,j−1i,j )x(Uψ+U
ψ
+U
ψ
+U
ψ
−4ψ)
xxyy2h ∂t22hi,j 2+η(1 − |ψ +η(1| )ψ−i,j|ψ, i,j |2 )ψi,j ,(16(∂Uxi,j1 i,j−1 i,ji,j i,j i+1,j 1i,ji,ji,j−1i,ji+1,ji,ji,j Lz − 1)Ui,jxi,j ,= −igIm[ψU]U− 2 (LL(17(x ψ − ]Ux (Lz= −igIm[ψUψ−1)U,xxzzx∂th2hi,j∂Uy1 i,j i−1,ji,j i,j i,j+1i,ji,j=i,j −igIm[ψUψ]U−(LL−1)U,(181i,jyyzzy
2 i−1,ji,j i,j+1i,ji,j∂th= −igIm[ψ Uy ψ]Uy − 2 (Lz Lz− 1)Uy ,(hy using link variables andariables and∂Uxi,j∂Ai,jxi,j−1i,ji,jLzLi,jz −113	x
Uyi,j,kUzi,j,k==xexp(−ihgAi,j,ky ),exp(−ihgAi,j,k),
z(43)

(44)

equa@on	  (	  3次元	  hereSemi-­‐discre@zed	  h is the latticespacing.
By)	
using these link variables, we can get the spatially
iscrete version of the equations (28) and (29).
The dynamics of the order parameters are described as follows:∂ψ i,j,k1
i−1,j,k i−1,j,ki,j−1,k i,j−1,ki,j,k i+1,j,ki,j,k i,j+1,kτ=(Uψ+Uψ+Uψ+Uψ
xyy∂t2h2 xi,j,k−1 i,j,k−1i,j,k i,j,k+1i,j,k+Uψ
+U− 6ψ) + η(1 − |ψ i,j,k |2 )ψ i,j,k .
(45)
The dynamics of thez link variablesz are ψ
described
as follows:∂Uxi,j,k
1 i,j,k i,j,k−1 i,j,k i,j−1,ki,j,k i,j,k i+1,j,k= −igIm[ψ Ux ψLz Lz− 1),8] − 2 (Ly Ly∂th∂Uyi,j,k1 i,j,ki,j,ki,j,k i,j,k−1= −igIm[ψ Uyi,j,k ψ i,j+1,k ] − 2 (Lz Li−1,j,kL− 1),zx Lz∂th∂Uzi,j,k1 i,j,ki,j,ki,j,k i−1,j,k= −igIm[ψ Uzi,j,k ψ i,j,k+1 ] − 2 (Lx Li,j−1,kL− 1),xy Ly∂th(4(4(4i,j,ki,j,khere the ”loop” variables Li,j,k
are defined asx , Ly , Lzi,j,k+1Li,j,k= Uyi,j,k Uzi,j+1,k U yxLi,j,ky=i,j,kUz,i+1,j,k i,j,kUzi,j,k Uxi,j,k+1 U z
Ux ,(414	(5境界条件	
境界での磁束密度指定	
(b)

(a)

hj = jB + 1

j = jB + 1hj = jBi=0i=1

j = jB

i=2(c)

i = Nx - 1 i = Nx i = Nx + 1(d)j = Ny + 1i = iB

j = Nyj=2
j = Ny - 1
j=1i = iB

j=015	

数値シミュレーション	• 動画：２次元	  横磁界	  • 動画：３次元	  横磁界，縦磁界，リング	  • リアルタイムデモ	  ２次元(40x80)，３次元(30x30x30)	  	  	  ※	  order	  parameter	  値に応じた色で可視化	  16	

コーディング（実装）	

17	表現に類似性をもたせる		 	 	 	 	 a	 =	 a	 +	 b・c・d	 +	 e・f	 	 	 	 	 	 	 a.equal(	 a.plus(	 b.mul(c).mul(d)	 ).plus(	 e.mul(f)	 );	 	 	 	 	 Complex	 plus(	 Complex	 z	 ){	 	 	 	 	 return	 new	 Complex(	 x	 +	 z.x,	 y	 +	 z.y	 );	 	 }	 	 

18	

「複素数クラス」を定義		 Complex	 mul(Complex	 z){	 	 	 	 	 double	 x1	 =	 z.x;	 	 	 	 	 double	 y1	 =	 z.y;	 	 	 	 	 double	 x2	 =	 x;	 	 	 	 	 double	 y2	 =	 y;	 	 	 	 	 return	 new	 Complex(x1*x2	 -	 y1*y2,	 x1*y2	 +	 y1*x2);	 	 	 }	 19	

i,j,ki,j,kp” variables Li,j,k,L,L
are defined asxyz

Li,j,kx=

Li,j,k=yLi,j,k=zi,j,k i,j+1,k i,j,k+1 i,j,kUy Uz
UyUz ,i,j,k i,j,k+1 i+1,j,k i,j,kUz Ux
UzUx ,i,j,k i+1,j,k i,j+1,k i,j,kUx UyUx
Uy .

physical meaning of the loop variable is expponential of the pa
Complex	 Lz(Complex[][][]	 Ux,	 Complex[][][]	 Uy,	 Complex[][][]	 Uz,	 int	 i,	 int	 j,	 int	 k){	 
ntial.	 	 	 return	 Ux[i][j][k].mul(Uy[i+1][j][k]).mul(Ux[i][j+1][k].conj()).mul(Uy[i][j][k].conj());	 
}	 

320	

The dynamics of the order parameters are described as follows:∂ψ i,j,k1
i−1,j,k i−1,j,ki,j−1,k i,j−1,ki,j,k i+1,j,ki,j,k i,j+1,kτ=(Uψ+U
ψ
+U
ψ
+U
ψ
xyy
∂t
2h2 x
i,j,k−1 i,j,k−1+Uzi,j,k ψ i,j,k+1 + U zψ− 6ψ i,j,k ) + η(1 − |ψ i,j,k |2 )ψ i,j,k .
(45//////////	 
void	 do_Laplacian(Complex[][][]	 psi,	 Complex[][][]	 Ux,	 Complex[][][]	 Uy,	 Complex[][][]	 Uz,	 double	 dt,	 double	 D){	 8	 	 for(int	 k	 =	 1;	 k	 <=	 Nz;	 k++){	 	 	 	 	 for(int	 j	 =	 1;	 j	 <=	 Ny;	 j++){	 	 	 	 	 	 	 for(int	 i	 =	 1;	 i	 <=	 Nx;	 i++){	 	 	 	 	 	 	 	 	 psi_s[i][j][k].equal(psi[i][j][k].plus(	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 	 ((psi[i+1][j][k].mul(Ux[i][j][k])).plus(psi[i-1][j][k].mul(Ux[i-1][j][k].conj()))	 	 	 	 	 	 	 	 	 	 	 	 	 	 .plus(psi[i][j+1][k].mul(Uy[i][j][k])).plus(psi[i][j-1][k].mul(Uy[i][j-1][k].conj()))	 	 	 	 	 	 	 	 	 	 	 	 	 	 .plus(psi[i][j][k+1].mul(Uz[i][j][k])).plus(psi[i][j][k-1].mul(Uz[i][j][k-1].conj()))	 	 	 	 	 	 	 	 	 	 	 	 	 	 .minus(psi[i][j][k].mul(6.0)	 )	 )	 	 	 	 	 	 	 	 	 	 	 	 	 	 .mul(dt	 *	 	 D	 /	 (h	 *	 h)	 )	 	 	 	 	 	 	 	 	 	 	 	 	 )	 	 	 	 	 	 	 	 	 );	 
ラプラシアン演算の実装		 	 	 	 	 	 }	 	 	 	 	 }	 	 	 }	 	 	 for(int	 k	 =	 1;	 k	 <=	 Nz;	 k++){	 	 	 	 	 for(int	 j	 =	 1;	 j	 <=	 Ny;	 j++){	 	 	 	 	 	 	 for(int	 i	 =	 1;	 i	 <=	 Nx;	 i++){	 	 	 	 	 	 	 	 	 psi[i][j][k].equal(psi_s[i][j][k]);	 	 	 	 	 	 	 }	 	 	 	 	 }	 	 	 }	 }	 //////////	 
21	

境界条件の実装		 ///	 constraints	 by	 the	 applied	 field	 and	 current	 	 	 ///	 y-z	 plane	 (	 Bz-Uy-component	 )	 set	 Ja	 =	 0	 when	 longitudinal	 	 	 for(int	 k	 =	 1;	 k	 <=	 Nz;	 k++){	 	 	 	 	 for(int	 j	 =	 1;	 j	 <=	 Ny-1;	 j++){	 	 	 	 	 	 	 ex.setPolar(1.0,	 	 h*h*(Ba	 +	 Ja*SSx/2.0));	 	 	 	 	 	 	 Uy[1][j][k].equal(Uy[2][j][k].mul(Ux[1][j+1][k].conj()).mul(Ux[1][j][k]).mul(ex));	 	 	 	 	 	 	 ex.setPolar(1.0,	 -h*h*(Ba	 -	 Ja*SSx/2.0));	 	 	 	 	 	 	 Uy[Nx][j][k].equal(Uy[Nx-1][j][k].mul(Ux[Nx-1][j+1][k]).mul(Ux[Nx-1][j][k].conj()).mul(ex));	 	 	 	 	 }	 	 	 }	 	 	 ///	 z-x	 plane	 (	 Bz-Ux-component	 )	 omit	 when	 transverse	 	 	 for(int	 i	 =	 1;	 i	 <=	 Nx-1;	 i++){	 	 	 	 	 for(int	 k	 =	 1;	 k	 <=	 Nz;	 k++){	 	 	 	 	 	 	 ex.setPolar(1.0,	 -	 h*h*Ba);	 	 	 	 	 	 	 Ux[i][1][k].equal(	 (Uy[i+1][1][k].conj()).mul(Ux[i][2][k]).mul(Uy[i][1][k]).mul(ex)	 );	 	 	 	 	 	 	 ex.setPolar(1.0,	 	 	 h*h*Ba);	 	 	 	 	 	 	 Ux[i][Ny][k].equal(	 Ux[1][Ny-1][k].mul(Uy[i+1][Ny-1][k]).mul(Uy[i][Ny-1][k].conj()).mul(ex)	 );	 	 	 	 	 }	 	 	 }	 	 	 ///	 y-z	 plane	 (	 J=(0,0,Jb)-Uz-component	 )	 longitudinal	 	 	 for(int	 k	 =	 1;	 k	 <=	 Nz-1;	 k++){	 	 	 	 	 for(int	 j	 =	 1;	 j	 <=	 Ny;	 j++){	 	 	 	 	 	 	 ex.setPolar(1.0,	 h*h*Jb*SSx/2.0);	 	 	 	 	 	 	 Uz[1][j][k].equal(	 Ux[1][j][k].mul(Uz[2][j][k]).mul(Ux[1][j][k+1].conj()).mul(ex)	 );	 	 	 	 	 	 	 Uz[Nx][j][k].equal(	 (Ux[Nx-1][j][k].conj()).mul(Uz[Nx-1][j][k]).mul(Ux[Nx-1][j][k+1]).mul(ex)	 );	 	 	 	 	 }	 	 	 }	 	 	 ///	 z-x	 plane	 (	 J=(0,0,Jb)-Uz-component	 )	 longitudinal	 	 	 for(int	 i	 =	 1;	 i	 <=	 Nx;	 i++){	 	 	 	 	 for(int	 k	 =	 1;	 k	 <=	 Nz-1;	 k++){	 	 	 	 	 	 	 ex.setPolar(1.0,	 h*h*Jb*SSy/2.0);	 	 	 	 	 	 	 Uz[i][1][k].equal(	 Uy[i][1][k].mul(Uz[i][2][k]).mul(Uy[i][1][k+1].conj()).mul(ex)	 );	 	 	 	 	 	 	 Uz[i][Ny][k].equal(	 (Uy[i][Ny-1][k].conj()).mul(Uz[i][Ny-1][k]).mul(Uy[i][Ny-1][k]).mul(ex)	 );	 	 	 	 	 }	 	 	 }	 22	

今後の課題	• 各種物理量の表示	  • アルゴリズムの改良	  • 並列化（高速化）	  	  • ピンニング	  • Two-­‐component	  系	  23