mumaxの使い方④ 磁場による磁壁駆動 (+困難軸磁気異方性)

埼玉大学 大学院理工学研究科 数理電子情報専攻 電気電子物理工学PG
最終更新日: 2025/03/11 作成日: 2025/03/11
このページの概要

マイクロマグネティックシミュレーションのソフトウェアmumax の実行サンプルとして磁場による磁壁駆動を解説（その２）．

TL;DR

マイクロマグネティックシミュレーションのソフトウェアmumax の実行サンプルとして磁場による磁壁駆動を解説（その２）．

大まかには以下のような流れでシミュレーションを実行するためのスクリプトを作り，シミュレーションを実行，解析する．

  1. 空間グリッドを決める．
  2. 物質パラメタ（material parameters）として交換相互作用の強さ，Gilbert減衰定数を指定する．
  3. 容易軸磁気異方性と困難軸磁気異方性をカスタム磁場として追加する．
  4. 磁化の初期条件（initial magnetization)としてTwoDomainを用いて2つの磁区があるように指定する．
  5. 磁気励起（excitation）として外部磁場ベクトルを指定する．
  6. 実行する計算を指定する．
  7. 計算を実行する．
  8. 結果を解析する．

模型

グリッドの配置は以下の図のようにする． セルのサイズは$L_0 = 1~\mathrm{nm}$の立方体とした．

グリッドの構造．1次元的にグリッドを配置している．反磁場の効果は考えないので，形状は計算量が最小限になるようにした．（$x$方向はもう少し短くても問題なさそうではある．）
図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という画像ファイルも生成される．）

bash
1

python3 analysis.py

磁壁の移動距離と移動速度．左) 移動距離．磁壁の位置が振動しながら移動している．右) 移動速度．位置の有限差分で速度を見積もっているためガタガタしている．
図2: 磁壁の移動距離と移動速度．左) 移動距離．磁壁の位置が振動しながら移動している．右) 移動速度．位置の有限差分で速度を見積もっているためガタガタしている．

磁壁の移動距離の解析解との比較．最初の立ち上がりのために少しズレが生じているがおおよそ同じふるまいをしている．
図3: 磁壁の移動距離の解析解との比較．最初の立ち上がりのために少しズレが生じているがおおよそ同じふるまいをしている．

図3を見ると，ほとんど解析解と同じふるまいをしているように見える．実際，横に$1~\mathrm{ns}$だけ並行移動させてみると，以下のようにほとんど重なる(図4)．

磁壁の移動距離の解析解との比較で，解析解の時間を$1~\mathrm{ns}$だけシフトさせたもの．
図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の使い方③ 磁場による磁壁駆動