TL;DR ¶
マイクロマグネティックシミュレーションのソフトウェアmumax
の実行サンプルとして磁場による磁壁駆動を解説（その２）．
大まかには以下のような流れでシミュレーションを実行するためのスクリプトを作り，シミュレーションを実行，解析する．
- 空間グリッドを決める．
- 物質パラメタ（material parameters）として交換相互作用の強さ，Gilbert減衰定数を指定する．
- 容易軸磁気異方性と困難軸磁気異方性をカスタム磁場として追加する．
- 磁化の初期条件（initial magnetization)として
TwoDomainを用いて2つの磁区があるように指定する．
- 磁気励起（excitation）として外部磁場ベクトルを指定する．
- 実行する計算を指定する．
- 計算を実行する．
- 結果を解析する．
模型 ¶
グリッドの配置は以下の図のようにする．
セルのサイズは$L_0 = 1~\mathrm{nm}$の立方体とした．
図1: グリッドの構造．1次元的にグリッドを配置している．反磁場の効果は考えないので，形状は計算量が最小限になるようにした．（$x$方向はもう少し短くても問題なさそうではある．）
模型として交換相互作用と磁気異方性（容易軸が$x$方向かつ困難軸が$z$方向）と外部磁場とのZeeman結合のみがあるものを考える．
形状磁気異方性を生む双極子間相互作用は考えない．
この模型の場合，rigidな磁壁のラグランジアンが多々良源著『スピントロニクス理論の基礎』5-5に与えられており，磁場による磁壁移動速度が
$$\begin{align}
B_w = \frac{\alpha K_{\perp} S}{2 \hbar \gamma}
\end{align}$$
を境に定性的に異なるふるまいをすることが見てとれるはずである（Walker breakdown）．
これをmumax3のインプットファイルで指定するパラメタで表現すると，$\frac{1}{2 a^3} K_{\perp} S^2 = K_z$, $\frac{\hbar \gamma S}{a^3} = M_s$なので，
$$\begin{align}
B_w = \frac{\alpha K_z}{M_s}
\end{align}$$
となる．
$\alpha = 0.02$, $K_z = 5 \times 10^4~\mathrm{J/m^3}$, $M_s = 8 \times 10^5~\mathrm{A/m}$とすると，$B_w = 1.25~\mathrm{mT}$となる．
$B < B_w$では磁壁速度は一定で，その大きさは
$$\begin{align}
\dot{X} = \frac{1}{\alpha} \lambda \gamma B
\end{align}$$
で与えられている．
交換スティフネス係数を$A_{ex} = 13 \times 10^{-12}\mathrm{J/m}$とし，容易軸磁気異方性を$K_x = 5 \times 10^5\mathrm{J/m^3}$とすると磁壁の幅$\lambda$は
$\lambda = 5.0~\mathrm{nm}$となる（非常に小さい）．
上述の$\alpha$を用いると，外部磁場が$B = 0.7~\mathrm{mT}$のときは，だいたい$\dot{X} \simeq 31~\mathrm{m/s}$となる．
$B > B_w$では磁壁は磁壁の面が回転しながら運動する．
その速度は
$$\begin{align}
\dot{X} = \frac{\lambda \gamma B}{\alpha} \left[
\frac{B}{B_w}
- \frac{1}{1 + \alpha^2} \frac{\left(\frac{B}{B_w}\right)^2 - 1}{\frac{B}{B_w} + \sin 2 \omega (t - t_0)}
\right]
,\end{align}$$
ここで
$$\begin{align}
\omega = \frac{\gamma B_w}{1 + \alpha^2} \sqrt{ \left(\frac{B}{B_w}\right)^2 - 1 }
.\end{align}$$
本には誤植がある（参考：武内修先生@筑波大のページ
）．
これ以上は理論について言及しない．
詳細は本や上記のウェブページを参考にしてください．
シミュレーションの実行 ¶
mumax3のインプットファイルを以下に示す．
ここでは外部磁場として$B =2~\mathrm{mT}$を$x$方向に加えている．
$B > B_w$なので振動しながら磁壁が移動する様子が見られるはずである．
domain_wall_motion_by_magnetic_field_with_hard_axis.mx3
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
|
Nx := 256
Ny := 1
Nz := 1
L0 := 1e-9 // [m]
setgridsize(Nx, Ny, Nz)
setcellsize(L0, L0, L0)
Msat = 800e3 // [A/m]
Aex = 13e-12 // [J/m]
alpha = 0.02
// no demagnetization field
enabledemag = false
// Custom anisotropy, easy, x-axis
Kx := 500e3
ux := Normalized(ConstVector(1, 0, 0))
Cx := Const( (2 * Kx) / (Msat.Average()))
Easy := Mul(Cx, Mul( Dot(ux, m), ux))
AddFieldTerm(Easy)
AddEdensTerm(Mul(Const(-0.5),Dot(Easy,M_full)))
// Custom anisotropy, hard, z-axis
Kz := -50e3
uz := Normalized(ConstVector(0, 0, 1))
Cz := Const( (2 * Kz) / (Msat.Average()))
Hard := Mul(Cz, Mul( Dot(uz, m), uz))
AddFieldTerm(Hard)
AddEdensTerm(Mul(Const(-0.5),Dot(Hard,M_full)))
m = TwoDomain(1,0,0, 0,1,0, -1,0,0) // Néel wall along y-axis
// Remove surface charges from left (mx=1) and right (mx=-1) sides to mimic infinitely long wire. We have to specify the region (0) at the boundaries.
BoundaryRegion := 0
MagLeft := 1
MagRight := -1
ext_rmSurfaceCharge(BoundaryRegion, MagLeft, MagRight)
relax()
saveas(m, "relaxed.ovf")
ext_centerWall(0) // keep m[0] (m_x) close to zero
autosave(m, 5e-11)
tableadd(ext_dwpos) // domain wall position
tableadd(ext_dwspeed) // doman wall speed
tableautosave(5e-11)
B_ext = vector(2e-3, 0e0, 0e0)
run(50e-9)
|
これを適当なディレクトリで以下のようにして実行する．
bash
1
|
mumax3 domain_wall_motion_by_magnetic_field_with_hard_axis.mx3
|
結果の解析 ¶
実行すると，domain_wall_motion_by_magnetic_field_with_hard_axis.out/table.txtというtsvファイルが生成される．
中身は以下のようになっている．
このデータを用いて磁壁の運動を考察してみることにする．
table.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
# t (s) mx () my () mz () ext_dwpos (m) ext_dwspeed (m/s)
0 2.3283064e-10 0.06237166 0 -0 0
5.0069999664518816e-11 0.0003953867 0.062361397 0.001086315 -0 0
1.0007482779614818e-10 0.0015401286 0.062331315 0.0021472154 -0 0
1.50042769337419e-10 0.0034148954 0.062282585 0.0031835493 -0 0
2.00009584522424e-10 0.0060021933 0.06221643 0.0041963174 -0 0
2.500871998449915e-10 0.009291801 0.062133953 0.0051880125 -0 0
3.000459795603734e-10 0.01324997 0.062036756 0.00615444 -0 0
3.5008461880384803e-10 0.0100605395 0.06192559 0.0070999013 1e-09 2.8564522
4.0004970468287127e-10 0.015315989 0.061801545 0.008021881 1e-09 2.8564522
4.500637327870739e-10 0.013385222 0.06166538 0.008923115 2e-09 10.0020895
5.000773920510917e-10 0.0120575195 0.06151813 0.009803152 3e-09 19.994537
5.500845679970129e-10 0.011313793 0.061360776 0.010662414 4e-09 19.99713
6.000685696641931e-10 0.011133273 0.061193764 0.011501016 5e-09 20.006401
6.500750626266121e-10 0.01150529 0.06101787 0.012320271 6e-09 19.997404
7.000583515268673e-10 0.012405019 0.060834132 0.013119962 7e-09 20.006687
7.500849830676537e-10 0.0138273705 0.060642783 0.01390167 8e-09 19.989353
8.000633786323985e-10 0.007926185 0.060444877 0.014664462 1e-08 40.017292
︙
中略
︙
4.925002361285198e-08 0.013246283 0.04656842 0.039549176 1.005e-06 40.04999
4.930005061174001e-08 0.010913078 0.046286575 0.039849006 1.008e-06 59.967617
4.9350080528947215e-08 0.008616895 0.046003588 0.040145803 1.011e-06 59.96412
4.9400082310063046e-08 0.014153225 0.04571919 0.04043927 1.013e-06 39.998573
4.9450018982125135e-08 0.011877321 0.04543463 0.040729895 1.016e-06 60.07609
4.950004589614359e-08 0.009663751 0.045148432 0.041018307 1.019e-06 59.96772
4.9550049500693897e-08 0.015273114 0.04486055 0.041303694 1.021e-06 39.997116
4.960007718608048e-08 0.013096633 0.044571627 0.041586958 1.024e-06 59.966797
4.9650013958487057e-08 0.0108942 0.044282105 0.041867346 1.027e-06 60.07597
4.9700043170603545e-08 0.008740031 0.043990523 0.042145863 1.03e-06 59.964966
4.9750045213655395e-08 0.014392081 0.043697156 0.042421628 1.032e-06 39.998367
4.9800072276682914e-08 0.0122437775 0.04340265 0.042695653 1.035e-06 59.96754
4.985000898223395e-08 0.010054808 0.04310715 0.042967077 1.038e-06 60.07605
4.9900043955160026e-08 0.007901781 0.042809162 0.04323689 1.041e-06 59.95806
4.9950040495402195e-08 0.013535401 0.04250936 0.043504264 1.043e-06 40.00277
|
具体的には，前回同様にpandasでこのファイルを読み込み，(pandasのコマンドを使わずに）matplotlibでグラフ化する．
さらに磁壁の移動速度の解析式をscipy.integrateで数値的に積分して移動距離を見積もり，それとシミュレーション結果と比較する．
analysis.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
|
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.integrate as integrate
g = 1.760e11 # rad / T s
K = 500e3 # J/m^3
K_perp = 50e3 # J/m^3
Ms = 800e3 # A/m
A = 13e-12 # J/m
alpha = 0.02
Bw = alpha * K_perp / Ms
l = np.sqrt(A/K)
B = 0.002
def v(t):
if B < Bw:
v = g*Bw*l/alpha*(t*0+1)
else:
w = g*Bw/(1+alpha**2) * np.sqrt((B/Bw)**2-1)
v = g*Bw*l/alpha * ( B/Bw - ((B/Bw)**2 - 1)/(B/Bw + np.sin(2*w*t))/(1+alpha**2) )
return v
def load_table(directory):
table = pd.read_csv(directory+"table.txt", sep='\t')
return table
def main():
directory = "./domain_wall_motion_by_magnetic_field_with_hard_axis.out/"
table = load_table(directory)
keys = table.keys()
fig, axes = plt.subplots(1,2,figsize=(10,5))
fig.tight_layout(pad=1.5)
axes[0].set_xlabel("time [s]")
axes[0].set_ylabel("distance [m]")
axes[0].plot(table[keys[0]], table[keys[-2]], label=keys[-2])
axes[1].set_xlabel("time [s]")
axes[1].set_ylabel("velocity [m/s]")
axes[1].plot(table[keys[0]], table[keys[-1]], label=keys[-1])
fig.savefig("result_B"+str(B)+"mT.png", dpi=300)
plt.show()
t = table[keys[0]].values
n = len(t)
d = np.empty(n)
for i in np.arange(n):
d[i], _ = integrate.quad(v, 0, t[i])
fig, ax = plt.subplots(1,1,figsize=(10,5))
ax.set_xlabel("time [s]")
ax.set_ylabel("distance [m]")
ax.plot(t, table[keys[-2]], label="numerical")
ax.plot(t, d, label="analytic")
fig.legend()
fig.savefig("comparison.png", dpi=300)
plt.show()
if __name__=="__main__":
main()
|
- 注意
numpyとpandasとmatplotlibとscipyは事前にpipなど経由で（以下のようなコマンドで）インストールしておく必要がある．
bash
1
|
pip install numpy pandas matplotlib scipy
|
これを以下のように実行すると，以下のようなグラフが得られる．
（それと同時にresult_B0.002T.pngとcomparison.pngという画像ファイルも生成される．）
図2: 磁壁の移動距離と移動速度．左) 移動距離．磁壁の位置が振動しながら移動している．右) 移動速度．位置の有限差分で速度を見積もっているためガタガタしている．
図3: 磁壁の移動距離の解析解との比較．最初の立ち上がりのために少しズレが生じているがおおよそ同じふるまいをしている．
図3を見ると，ほとんど解析解と同じふるまいをしているように見える．実際，横に$1~\mathrm{ns}$だけ並行移動させてみると，以下のようにほとんど重なる(図4)．
図4: 磁壁の移動距離の解析解との比較で，解析解の時間を$1~\mathrm{ns}$だけシフトさせたもの．
「ほとんど」ということは差があることを意味するが，この差は，解析解がrigidな磁壁を仮定しているところから生じていると思われる．
実際，磁気異方性として困難軸に$|K_z| = 5 \times 10^4~\mathrm{J/m^3}$としており，容易軸の$|K_x| = 5 \times 10^5~\mathrm{J/m^3}$に比べて一桁小さいだけの場合をシミュレーションしたので，理論が要請する$|K_z| \ll |K_x|$という条件には少し無理がありそうである．
$|K_z|$をもう少し小さくしてシミュレーションをすると，より解析解と一致するふるまいが見られるのかもしれない．
前の記事: mumaxの使い方③ 磁場による磁壁駆動