物理でよくある歳差運動のアニメーションを作ってみた【Julia】

お前は誰？

私は物性物理学を研究している大学院生です。

なんで歳差運動のムービーを作成するの？

学会発表をするとき、スライドに少しでもかっこいい素材を載せたくてムービーを作成してみました。私が研究でよく使用する方程式は、物理的に「スピンの歳差運動」であると解釈できます。非物理学徒のために簡単に説明すると、スピンとはベクトル（矢印）のことだと考えて頂いて構いません。この方程式の物理的意味をわかりやすく理解するには、スピン歳差運動の様子を視覚的に示すのがベストかなと考えた次第です。

本題

本記事の主役

私はLLG方程式という時間発展方程式を使用しています。簡単のため、1つのスピン $ \boldsymbol{S}(t) $ しか存在しない状況を考えます。スピンがベクトルなのは、3つの成分を持つからです。

\boldsymbol{S}(t)
= [S_x(t), S_y(t), S_z(t)]

さて、このスピンの時間発展を表す方程式であるLLG方程式は以下のように書き下せます。

\frac{d\boldsymbol{S}(t)}{dt}
= \frac{1}{\hbar} \boldsymbol{S} \times \boldsymbol{B}(t) + \frac{\alpha}{\left|\boldsymbol{S}(t)\right|} \boldsymbol{S}(t) \times \frac{d\boldsymbol{S}(t)}{dt}

ここで、$ \boldsymbol{B} $ は物理系に加えられた磁場です。$ \hbar, \alpha $は定数です。各項の説明を行うと…

• 左辺 $ \frac{d\boldsymbol{S}(t)}{dt} $…これはスピンに時間発展を表します。スピンが変化する方向を向いています。
• 右辺第1項 $ \boldsymbol{S} \times \boldsymbol{B}(t) $…これは歳差運動を表します。スピンと磁場の外積がどの方向を向くのか？を考えるとわかります。
• 右辺第2項 $ \boldsymbol{S}(t) \times \frac{d\boldsymbol{S}(t)}{dt} $…これは緩和を表します。歳差運動の回転軸方向を向くベクトルです。

これらをJuliaでプログラムしたコードが以下になります。

ソースコード

using GLMakie

#---------- ベクトルの設定 ----------#
# 外部磁場ベクトルの定義
eb = [0, 0, 1.0] # 磁場ベクトル
spin(θ, ϕ) = [sin(θ)cos(ϕ), sin(θ)sin(ϕ), cos(θ)] # 角度で指定されるスピンの関数
M = spin(π/4, 0.0) # 初期時刻でのスピン

# 交叉積の計算
v1 = cross(B, M) # 右辺第1項
v2 = cross(M, v1) # 右辺第2項

#---------- アニメーションの作成 ----------#
# 図のサイズの設定
fig_width = 600 * 2
fig_hight = 550 * 2

# 図の作成
fig = Figure(size=(fig_width, fig_hight))
ax = Axis3(fig[1, 1],
    protrusions=(0, 0, 0, 0),
    viewmode=:fitzoom,
    elevation= 0,
    azimuth= 0,
    aspect = :data
)

# 描画範囲の設定
GLMakie.xlims!(ax, -1.0, 1.0)
GLMakie.ylims!(ax, -1.0, 1.0)

# ベクトルの描画
ar_b = arrows!(
        ax,
        [Point3f(0.0, 0.0, 0.0)],
        [Vec3f(eb)],
        fxaa=true,
        color=:red,
        linewidth=0.05,
        arrowsize=Vec3f(0.15, 0.15, 0.2),
        quality=30,
    )

ar_m = arrows!(
        ax,
        [Point3f(0.0, 0.0, 0.0)],
        [Vec3f(M)],
        fxaa=true,
        color=:orange,
        linewidth=0.05,
        arrowsize=Vec3f(0.15, 0.15, 0.2),
        quality=30,
    )

ar_v1 = arrows!(
        ax,
        [Point3f(M)],
        [Vec3f(v1 ./ 2)],
        fxaa=true,
        color=:green,
        linewidth=0.05,
        arrowsize=Vec3f(0.15, 0.15, 0.2),
        quality=30,
    )

ar_v2 = arrows!(
        ax,
        [Point3f(M)],
        [Vec3f(v2 ./ 2)],
        fxaa=true,
        color=:blue,
        linewidth=0.05,
        arrowsize=Vec3f(0.15, 0.15, 0.2),
        quality=30,
    )

hidedecorations!(ax)
hidespines!(ax)

# フレーム数とフレームレートの設定
frameEnd = 300
frames = 1:frameEnd
framerate = 60

# アニメーションの記録
record(fig, "LLGeq.mp4", frames; framerate = framerate) do frame
    # フレームに応じたφとθの計算
    φ = 8π * (frame / length(frames))
    θ = (π / 4) * (1 - frame / length(frames))
    
    # スピンの計算
    M = spin(θ, φ)
    
    # 交叉積の計算
    v1 = cross(B, M)
    v2 = cross(M, v1)

    ar_m[2] = [Vec3f(M)]

    if frame == frameEnd
        ar_v1[1] = [Point3f(M)]
        ar_v1[2] = [Vec3f([0.0, 0.0, 0.0])]

        ar_v2[1] = [Point3f(M)]
        ar_v2[2] = [Vec3f([0.0, 0.0, 0.0])]

    else
        ar_v1[1] = [Point3f(M)]
        ar_v1[2] = [Vec3f(v1 ./ 2)]

        ar_v2[1] = [Point3f(M)]
        ar_v2[2] = [Vec3f(v2 ./ 2)]
    end

end

これを実行して得られたムービーが以下です。

• 赤い矢印：磁場
• オレンジの矢印：スピン
• 緑の矢印：歳差運動の方向
• 青の矢印：緩和の方向

ということで、スピンが磁場を回転軸としながら歳差運動をして、最後は緩和する様子が確認できましたね。arrowsは凡例Legnedに対応していなそうだったので、とりあえずこれで完成ということで。