スピントロニクス理論の基礎/X-6 の履歴(No.5)
更新教科書4章〜5章で用いられた1次元スピン結晶のモデルを コンピュータシミュレーションしてみたい。
まずは運動方程式†
連続近似をする前の式からハミルトニアンを求めれば、
&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 );
を得て、これがスピンの運動方程式となる。
- z軸が容易軸
- y軸が困難軸
- z軸に平行に外部磁場 B
シミュレーション (基本)†
- sx[i], sy[i], sz[i] : i番目のスピンの の x, y, z 成分
- dhx[i], dhy[i], dhz[i] : i番目のスピンの の x, y, z 成分
- dsx[i], dsy[i], dsz[i] : i番目のスピンの の x, y, z 成分
- dt : シミュレーション時間間隔
とすれば、シミュレーションの1回のイタレーションを
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 はそれぞれ対応する物理パラメータの次元と係数を少しごまかした物。
シミュレーション (注意点)†
上記計算を繰り返すと、dt を有限に取っているせいで の大きさが元の値から徐々にずれてしまう。
そこで、何回かイタレーションするたびに、 を正規化して正しい大きさに戻してやることにする。
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番目のスピンの大きさ
見られた現象†
計算では全てのスピンの大きさを1として、 ランダムな初期値からの時間発展を追ってみた。
散逸を小さくすると計算エラーが蓄積して発散する傾向にあることが分かった。 ではなく で記述した方が良いのかもしれない。 散逸がある程度大きければ、何とかそれらしい結果が得られるみたい。
- 困難軸異方性のあるときの性質
- 静止磁壁には2種類ある
- x成分が正になる物と負になる物
- 慣性が働く
- 外部磁場を切ってもスピンのy軸成分に溜まったエネルギーが無くなるまで進み続ける
- 磁場の弱いときは並進し、磁場の強いときは振動しながら並進する
- 静止磁壁には2種類ある
- 困難軸異方性のあるときの磁壁間相互作用
- x成分の正負で相互作用が異なる
- 同じ符号の磁壁間には引力が働く
- 異なる符号の磁壁間には斥力が働く
- 外部磁場の中では異なる符号の磁壁がペアになった「磁壁対」が安定構造になる
- 強い磁場では振動によりペアが生成される前に打ち消し合ってしまう
- 弱い磁場では符号が決定してからペアが生成されるために磁壁対が自己形成する
- x成分の正負で相互作用が異なる
などが見られた。
計算例†
困難軸異方性あり、外部磁場なし†
J = 0.1 K = 0.01 Kp= 0.005 A = 0.1 B = 0.000 dt = 0.2
- Sy がゼロでない磁壁は運動エネルギーを持ち、摩擦でエネルギーを失うまで動き続ける
- 最終的に一定の厚さを持つたくさんの磁壁が残り、「磁区」が形成される
- 最終的に残った磁壁は Sx の符号により正・負に分類できる
- 緩和の過程で、隣り合う磁壁が結合して消える様子が見られる
- 隣り合う磁壁は同じ符号を持つときだけ結合して消滅する
- 異なる符号を持つ磁壁は互いに反発し合うため結合することはない
- 磁壁が消滅する場合など、エネルギーが失われる際には大きな歳差運動が起きて、摩擦項経由でエネルギーが散逸する
容易軸異方性を変えてみる†
J = 0.1 K = 0.001 // <= 0.01 Kp= 0.005 A = 0.1 B = 0.000 dt = 0.2
- が増加して、磁壁が厚くなった
- 符号の同じ磁壁は引力相互作用を感じ、融合すると消滅する
- 符号の異なる磁壁は斥力相互作用を感じ、反発し合って距離が広がる
容易軸異方性をさらに変えてみる†
K の値を、0.001 -> 0.003 -> 0.01 -> 0.003 -> 0.001 -> 0.0003 と変化させた。
- に比例して、磁壁の厚さが変化する
- 符号の異なる磁壁は斥力相互作用を感じ、反発し合って距離が広がる
ここに上げた動画†
- 数値計算を Igor のマクロで行う
- 計算途中のグラフを Igor の SavePICT コマンドで png として保存する
- ffmpeg -i "spin%d.png" -sameq spin.mp4 で mp4 にする
- FLAVER 3.0 で表示
の手順で作成した。
Windows7 では IE, Chrome, Firefox で動作確認できた。
Android では再生はされたがコマ落ちがひどかったので、 もう少し良いコーデックがないかを探すべきかも。
時間発展の精密化†
に を足す部分、 もう少しうまく足すことで誤差を減らせないか検討してみる。
の時間変化として を加える計算では、
この図に見るように、
- 本来円弧方向に の長さ進むべきが、 垂直方向に進んでしまうため中心角の変化量がずれてしまう
- 進んだ先は必ずスピンベクトルの大きさが大きくなってしまう。
という2つの問題が生じる。
これを何とかしたい。
中心角について†
計算では簡単のため としておく。
すると本来の角度変化は、弧の長さが なので、ラジアンで表せばそのまま
で与えられるのだが、単に
としてしまうと、角度変化は
になってしまう。
に注意しつつ、角度変化量を正しくするには、
&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} );
とすればよい。
長さについて†
による中心角変化は、
である。この を使えば
となるが、
を使えば、 を正規化するには、
&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));
まとめると†
&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} );
という計算をすればよいことになる。
補正の効果†
の大きさに対する効果を見た。
かなり効果はありそうだ。
足し算直後に規格化をした場合、角度誤差のみが効くけれど、 そちらへの効果のみを見ても、下のように絶大な効果がある。
足し算時に補正を入れつつ、 を少し大きめに取れるなら、計算コスト的にも得になりそう?
コードとしては†
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]
とすればよい。
長距離相互作用について†
このモデルではスピン・スピン相互作用は再隣接間に限られている。
磁気力は長距離力なので、特に強磁性や反強磁性ではより長距離の相互作用を入れた方が良い場合がありそうに思える。
を計算するときに遠くからの寄与を入れれば良いのだが、 普通にやったら計算量が から に上がってしまうので、 遠くの効果は少し平均して 程度に抑えたり・・・できるのかな。