スピントロニクス理論の基礎/X-6 の変更点

更新 
[[スピントロニクス理論の基礎/5-6]]

#contents

教科書４章〜５章で用いられた１次元スピン結晶のモデルを
コンピュータシミュレーションしてみたい。

* まずは運動方程式 [#oaa25af4]

連続近似をする前の式からハミルトニアンを求めれば、

&math(H&=\sum_{\bm r}\left[\left\{-J_0\sum_{\bm a}\bm S(\bm r)\cdot\bm S(\bm r+\bm a)\right\}-\frac{1}{2}KS_z(\bm r)^2+\frac{1}{2}K_\perp S_y(\bm r)^2+\gamma B \hbar S_z(\bm r)\right]\\&=\sum_{\bm r}\left[\bm S(\bm r)\cdot\left\{\gamma B \hbar\bm e_z-J_0\sum_{\bm a}\bm S(\bm r+\bm a)\right\}-\frac{1}{2}KS_z(\bm r)^2+\frac{1}{2}K_\perp S_y(\bm r)^2\right]\\
);

ここから、

&math(
\frac{\PD H}{\PD \bm S(\bm r)}&=-J_0\left\{\sum_{\bm a}\bm S(\bm r+\bm a)\right\}-K S_z(\bm r)\bm e_z+K_\perp S_y(\bm r)\bm e_y+\gamma B \hbar\bm e_z
);

により、

&math(
\dot{\bm S}&=\frac{1}{\hbar}\frac{\PD H}{\PD \bm S}\times \bm S-\frac{1}{\hbar}\frac{\alpha}{S}\bm S\times\dot{\bm S}\\&=\frac{1}{\hbar}\left[-J_0\left\{\sum_{\bm a}\bm S(\bm r+\bm a)\right\}-K S_z(\bm r)\bm e_z+K_\perp S_y(\bm r)\bm e_y+\gamma B \hbar\bm e_z+\frac{\alpha}{S}\dot{\bm S}\right]\times \bm S
);

あるいは &math(\bm B,\bm K,\bm K_\perp); をベクトルとして、

&math(
\dot{\bm S}=\frac{1}{\hbar}\left[-J_0\left\{\sum_{\bm a}\bm S(\bm r+\bm a)\right\}-\bm K\cdot \bm S(\bm r)+\bm K_\perp \cdot\bm S(\bm r)+\hbar\gamma \bm B+\frac{\alpha}{S}\dot{\bm S}\right]\times \bm S
);

を得、これが単純立方格子におけるスピンの運動方程式となる。

以下では教科書と同様に、

- z軸が容易軸
- y軸が困難軸
- z軸に平行に外部磁場 B

である。

* シミュレーション (基本) [#a6861289]

- sx[i], sy[i], sz[i] : i番目のスピンの &math(\bm S); の x, y, z 成分
- dhx[i], dhy[i], dhz[i] : i番目のスピンの &math(\PD H/\PD\bm S); の x, y, z 成分
- dsx[i], dsy[i], dsz[i] : i番目のスピンの &math(\dot {\bm S}); の x, y, z 成分
- dt : シミュレーション時間間隔

とすれば、シミュレーションの１回のイタレーションを

 LANG:C
 // 有効磁場を計算
 dhx[i] = - J * (sx[i-1]+sx[i+1]) + A * dsx[i]
 dhy[i] = - J * (sy[i-1]+sy[i+1]) + A * dsy[i] + Kp * sy[i]
 dhz[i] = - J * (sz[i-1]+sz[i+1]) + A * dsz[i] - K  * sz[i] + B
 
 // トルクを計算
 dsx[i] = dhy[i] * sz[i] - dhz[i] * sy[i]
 dsy[i] = dhz[i] * sx[i] - dhx[i] * sz[i]
 dsz[i] = dhx[i] * sy[i] - dhy[i] * sx[i]
 
 // 時間発展
 sx[i]+= dt * dsx[i]
 sy[i]+= dt * dsy[i]
 sz[i]+= dt * dsz[i]

とすれば良いことになる。~
ただし J, A, Kp, K, B はそれぞれ対応する物理パラメータの次元と係数を少しごまかした物。

* シミュレーション (注意点) [#m0fbd358]

上記計算を繰り返すと、dt を有限に取っているせいで &math(\bm S); 
の大きさが元の値から徐々にずれてしまう。

そこで、何回かイタレーションするたびに、&math(\bm S); 
を正規化して正しい大きさに戻してやることにする。

 LANG:C
 // ノルムの逆数を計算
 c = ss[i] != 0 ? ss[i]/sqrt(sx[i]^2+sy[i]^2+sz[i]^2) : 1
 
 // 正規化
 sx[i] *= c
 sy[i] *= c
 sz[i] *= c

- ss[i] : i番目のスピンの大きさ

* 見られた現象 [#e753c83c]

計算では全てのスピンの大きさを１として、
ランダムな初期値からの時間発展を追ってみた。

散逸を小さくすると計算エラーが蓄積して発散する傾向にあることが分かった。
&math(x,y,z);ではなく&math(\theta,\phi);で記述した方が良いのかもしれない。
散逸がある程度大きければ、何とかそれらしい結果が得られるみたい。

- 困難軸異方性のあるときの性質
-- 静止磁壁には２種類ある
--- x成分が正になる物と負になる物
-- 慣性が働く
--- 外部磁場を切ってもスピンのy軸成分に溜まったエネルギーが無くなるまで進み続ける
-- 磁場の弱いときは並進し、磁場の強いときは振動しながら並進する
- 困難軸異方性のあるときの磁壁間相互作用
-- x成分の正負で相互作用が異なる
--- 同じ符号の磁壁間には引力が働く
--- 異なる符号の磁壁間には斥力が働く
-- 外部磁場の中では異なる符号の磁壁がペアになった「磁壁対」が安定構造になる
--- 強い磁場では振動によりペアが生成される前に打ち消し合ってしまう
--- 弱い磁場では符号が決定してからペアが生成されるために磁壁対が自己形成する

などが見られた。

* 計算例 [#kdff7d39]

** 困難軸異方性あり、外部磁場なし [#j93ae1cd]

 J = 0.1
 K = 0.01
 Kp= 0.005
 A = 0.1
 B = 0.000
 dt = 0.2

&attachref(spin1.mp4,,535x575); &attachref(spin3000.png,,525x525);

- Sy がゼロでない磁壁は運動エネルギーを持ち、摩擦でエネルギーを失うまで動き続ける
- 最終的に一定の厚さを持つたくさんの磁壁が残り、「磁区」が形成される
-- 磁壁の厚さは一定だが、磁区の厚さはランダムである
- 最終的に残った磁壁は Sx の符号により正・負に分類できる
- 緩和の過程で、隣り合う磁壁が結合して消える様子が見られる
- 隣り合う磁壁は同じ符号を持つときだけ結合して消滅する
- 異なる符号を持つ磁壁は互いに反発し合うため結合することはない
- 磁壁が消滅する場合など、エネルギーが失われる際には大きな歳差運動が起きて、摩擦項経由でエネルギーが散逸する

** 容易軸異方性を変えてみる [#b610d18d]

K や J の大きさを変えると磁壁の厚さが変化する。

 J = 0.1
 K = 0.001 //  <= 0.01
 Kp= 0.005
 A = 0.1
 B = 0.000
 dt = 0.2

&attachref(spin2.mp4,,535x575); &attachref(spin4500.png,,525x525);

- &math(\lambda=a \sqrt{K/J_0}); が増加して、磁壁が厚くなった
- 磁壁の厚さの変化と同時に、磁壁間に働く引力・斥力相互作用の及ぶ範囲も変化する
- 符号の同じ磁壁は引力相互作用を感じ、融合すると消滅する
- 符号の異なる磁壁は斥力相互作用を感じ、反発し合って距離が広がる

** 容易軸異方性をさらに変えてみる [#b58fee87]

K の値を、0.001 -> 0.003 -> 0.01 -> 0.003 -> 0.001 -> 0.0003 と変化させた。

&attachref(spin3.mp4,,535x575); &attachref(,,525x525);

- &math(\lambda=a \sqrt{K/J_0}); に比例して、磁壁の厚さが変化する
- 符号の異なる磁壁は斥力相互作用を感じ、反発し合って距離が広がる

** 外部磁場で磁壁を加速、磁壁同士の玉突きによる運動量のやりとり [#s0bd0288]

- &math(K=0.003); に戻して 30s → 磁壁を薄くして間隔を開ける
- &math(p>780); の黒枠で囲った領域のみ &math(B=-0.001); として磁壁を加速
- 外部磁場を切ると同時に &math(A=0.01); として摩擦を減らし、慣性による運動を見た。
-- 摩擦が小さいと計算誤差の蓄積が気になるため、ここでは &math(\Delta t=0.05); としている

&attachref(spin4.mp4,,535x575); 

- 外部磁場により徐々に磁壁が加速される
- 磁場を切った後も慣性＋摩擦により磁壁は徐々に減速しながら運動を続ける
- 符号の異なる磁壁間に働く斥力により、磁壁の玉突きが見られる
-- 動いていた磁壁が停止して
-- 止まっていた磁壁が同じ速度で動き出す
- 次々に玉突きが起き、最後には左端の磁壁まで運動量が伝わった

** 磁壁対 [#e7493dd7]

弱い外部磁場の元で、符号の異なる磁壁が押付けあわされて磁壁対を形成する。

その後磁場を切ると磁壁同士の斥力により磁壁は互いに離れる方向へ動き出す。(5:00以降)

&attachref(spin5.mp4,,535x575); 

- 外部磁場により隣り合う磁壁が接近する
- 符号の異なる磁壁間に働く斥力によりバウンドする
- 徐々にエネルギーが失われ、安定な「磁壁対」が形成される
- 磁場が無くなると斥力によりペアになった磁壁は反対方向へと押し戻される

** 強い磁場の下では [#b9d9ae37]

強い磁場の下では磁壁の運動は振動を伴うようになる。

完全に位相が逆でない限り、隣り合う磁壁間に斥力は働かず、融合して消えるため、
強い磁場の下で磁壁対は自然形成されない。

ただし、一旦形成された磁壁対はそのような強い磁場の下でも安定である。

* ここに上げた動画 [#wf4540f8]

+ 数値計算を Igor のマクロで行う
+ 計算途中のグラフを Igor の SavePICT コマンドで png として保存する
+ ffmpeg -i "spin%d.png" -sameq spin.mp4 で mp4 にする
+ FLAVER 3.0 で表示

の手順で作成した。

Windows7 では IE, Chrome, Firefox で動作確認できた。

Android では再生はされたがコマ落ちがひどかったので、
もう少し良いコーデックがないかを探すべきかも。

* 時間発展の精密化 [#i509acaa]

&math(\bm S); に &math(\Delta t\dot{\bm S}); を足す部分、
もう少しうまく足すことで誤差を減らせないか検討してみる。

&math(\bm S); の時間変化として &math(\Delta t\dot{\bm S}); を加える計算では、

&attachref(vector.png);

この図に見るように、

- 本来円弧方向に &math(|\Delta t\dot{\bm S}|); の長さ進むべきが、
垂直方向に進んでしまうため中心角の変化量がずれてしまう
- 進んだ先は必ずスピンベクトルの大きさが大きくなってしまう。

という２つの問題が生じる。

これを何とかしたい。

** 中心角について [#y685977e]

計算では簡単のため &math(|\bm S|=1); としておく。

すると本来の角度変化は、弧の長さが &math(|\Delta t\dot{\bm S}|); なので、ラジアンで表せばそのまま

&math(|\Delta t\dot{\bm S}|);

で与えられるのだが、単に

&math(\bm S+\Delta t\dot{\bm S});

としてしまうと、角度変化は

&math(\arctan |\Delta t\dot{\bm S}|);

になってしまう。

&math(\arctan x=x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+O(x)^9);

に注意しつつ、角度変化量を正しくするには、

&math(
\delta \bm S=\frac{|\Delta t\dot{\bm S}|}{\arctan |\Delta t\dot{\bm S}|}\Delta t\dot{\bm S}
\sim \frac{1}{1-\frac{|\Delta t\dot{\bm S}|^2}{3}}\Delta t\dot{\bm S}
\sim \left(1+\frac{|\Delta t\dot{\bm S}|^2}{3}\right)\Delta t\dot{\bm S}
);

とすればよい。

** 長さについて [#u6225e10]

&math(\bm S\rightarrow\bm S+\delta \bm S); 

による中心角変化は、

&math(\theta=\arctan |\delta \bm S|);

である。この &math(\theta); を使えば

&math(|\bm S+\delta \bm S|=\frac{1}{\cos\theta}); 

となるが、

&math(\frac{1}{\cos\theta}=\sqrt{1+\tan^2\theta}=\sqrt{1+|\delta \bm S|^2});

を使えば、&math(\bm S+\delta \bm S); を正規化するには、

&math(\frac{\bm S+\delta \bm S}{|\bm S+\delta \bm S|}
=\frac{\bm S+\delta \bm S}{\sqrt{1+|\delta \bm S|^2}}
=\left(1-\frac{|\delta \bm S|^2}{2}\right)(\bm S+\delta \bm S)); 

とすれば良いことが分かる。

&math(
|\delta \bm S| = \left(1+\frac{|\Delta t\dot{\bm S}|^2}{3}\right)|\Delta t\dot{\bm S}|
);

&math(
|\delta \bm S|^2 = \left(1+\frac{|\Delta t\dot{\bm S}|^2}{3}\right)^2|\Delta t\dot{\bm S}|^2\sim |\Delta t\dot{\bm S}|^2
);

なので結局、

&math(\frac{\bm S+\delta \bm S}{|\bm S+\delta \bm S|}
=\left(1-\frac{|\Delta t\dot{\bm S}|^2}{2}\right)(\bm S+\delta \bm S)); 

** まとめると [#n6c754a0]

&math(
&\delta \bm S = \left(1+\frac{|\Delta t\dot{\bm S}|^2}{3}\right)\Delta t\dot{\bm S}\\
&\bm S\leftarrow\left(1-\frac{|\Delta t\dot{\bm S}|^2}{2}\right)(\bm S+\delta \bm S)\\
\\
&\therefore
\bm S\leftarrow\left(1-\frac{|\Delta t\dot{\bm S}|^2}{2}\right)\bm S
+\left(1-\frac{|\Delta t\dot{\bm S}|^2}{6}\right)\Delta t\dot{\bm S}
);

という計算をすればよいことになる。

** 補正の効果 [#e1995e04]

&math(|\Delta t \dot{\bm S}|); の大きさに対する効果を見た。

&attachref(スピン計算誤差１.png);

まずまず効果はありそうだ。

足し算直後に規格化をした場合、角度誤差のみが効くけれど、
そちらへの効果のみを見ても、かなりの効果がある。

&attachref(スピン計算誤差２.png);

** コードとしては [#d7c688d8]

 LANG:C
 // 有効磁場を計算
 dhx[i] = - J * (sx[i-1]+sx[i+1]) + A * dsx[i]
 dhy[i] = - J * (sy[i-1]+sy[i+1]) + A * dsy[i] + Kp * sy[i]
 dhz[i] = - J * (sz[i-1]+sz[i+1]) + A * dsz[i] - K  * sz[i] + B
 
 // トルクを計算
 dsx[i] = dhy[i] * sz[i] - dhz[i] * sy[i]
 dsy[i] = dhz[i] * sx[i] - dhx[i] * sz[i]
 dsz[i] = dhx[i] * sy[i] - dhy[i] * sx[i]
 
 // 時間発展
 ds2[i] = dt^2 (dsx[i]^2+dsy[i]^2+dsz[i]^2)
 
 sx[i] = (1-ds2[i]/2) * sx[i] + (1-ds2[i]/6) * dt * dsx[i]
 sy[i] = (1-ds2[i]/2) * sy[i] + (1-ds2[i]/6) * dt * dsy[i]
 sz[i] = (1-ds2[i]/2) * sz[i] + (1-ds2[i]/6) * dt * dsz[i]

とすればよい。

** でも、あんまり意味ないかな？ [#lf6ea163]

&math(\Delta t); の最大値はこの誤差では決まってこないようで、
補正にはあまり意味が無いかもしれないと思えてきた。

&math(\Delta t\dot{\bm S}<0.001); にして数十回に１回くらい？
規格化する方が計算コストを下げられるのかも。

* 長距離相互作用について [#f2d32755]

このモデルではスピン・スピン相互作用は再隣接間に限られている。

磁気力は長距離力なので、特に強磁性や反強磁性ではより長距離の相互作用を入れた方が良い場合がありそうに思える。

&math(\frac{\PD H}{\PD \bm S}); を計算するときに遠くからの寄与を入れれば良いのだが、
普通にやったら計算量が &math(\Theta(n)); から &math(\Theta(n^2)); に上がってしまうので、
遠くの効果は少し平均して &math(\Theta(n\log n)); 程度に抑えたり・・・できるのかな。

* コメント・質問 [#dbe28f3c]

#article_kcaptcha
Counter: 6853 (from 2010/06/03), today: 1, yesterday: 5