🎸

Juliaによる非調和振動子の可視化

10 min read

モースポテンシャルの中を運動する粒子のSchrödinger方程式の解を可視化する.

パッケージ

Juliaのパッケージモードでadd Plots, add SpecialPolynomials, add SpecialFunctionsを実行し, 事前にパッケージをインストールしておく必要がある. また, ノート上ではusing Plots, using SpecialPolynomials, using Functionsを宣言する.

パッケージ
# using Pkg
# Pkg.add("Plots")
# Pkg.add("SpecialPolynomials")
# Pkg.add("SpecialFunctions")
using Plots
using SpecialPolynomials
using SpecialFunctions

時間に依存しないSchrödinger方程式

モースポテンシャル

V(r) = D_\mathrm{e} \left( \mathrm{e}^{-2a(r-r_e)} - 2\mathrm{e}^{-a(r-r_e)} \right)

の中で運動する質量\muの粒子のハミルトニアンは

\hat{H} = \left[ - \frac{\hbar^2}{2\mu} \frac{\mathrm{d}^2}{\mathrm{d}r ^2} + V(r) \right]

であり, これに対する時間に依存しないSchrödinger方程式

\hat{H}\psi(r) = E \psi(r)

の解は

\psi_n(\xi) = N_n z^{\lambda-n-1/2} \mathrm{e}^{-z/2} L_n^{(2\lambda-2n-1)}(\xi)\\ E_n = -D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right) - \chi \hbar \omega \left( n + \frac{1}{2} \right)^2

である. ただし, \xi=2\lambda\mathrm{e}^{-a(r-r_e)}, \omega=\sqrt{k/µ}, k=2D_\mathrm{e}a^2, \lambda=\frac{\sqrt{2mD_\mathrm{e}}}{a\hbar}, \chi=\frac{\hbar\omega}{4D_\mathrm{e}}, N_n=\sqrt{\frac{n!(2\lambda-2n-1)a}{\Gamma(2\lambda-n)}}であり, L_n^{(\alpha)}は一般化ラゲール多項式である. ラゲールの陪多項式ではないので注意されたい. また, ほぼWikipediaの通りだが, 規格化定数N_naが足りないと思われるので書き加えた.

ポテンシャル
V(r,re,De,a) = De*( exp(-2*a*(r-re)) -2*exp(-a*(r-re)) )
固有関数
# using SpecialFunctions
# using SpecialPolynomials
Γ(z) = try; gamma(z); catch; Inf; end

function ψ(n,r,re,De,a,µ)= 1.0
    λ = sqrt(2*µ*De)/(a*)
    ξ = 2*λ*exp(-a*(r-re))                                    
    s = 2*λ-2*n-1
    N = sqrt(factorial(n)s*a/Γ(s+n+1))
    return N*ξ^(s/2)*exp(-ξ/2)*basis(Laguerre{s}, n)(ξ)
end
固有値
function E(n,De,a,µ)= 1
    k = 2*a^2*De
    ω = sqrt(k/µ)
    χ =*ω/(4*De)
    return - De +*ω*(n+1/2) - χ**ω*(n+1/2)^2
end
パラメータ
re = 1.38436
De = 1.13133356259-1
k = 12.772/15.569141/2
a = sqrt(k/(2*De))
µ = 1/(1/1837.108542 + 1/1837.108542)

規格化条件

固有関数が正しく宣言できているのか見た目だけでは判断できないことが多いので, ライブラリの使い方が正しいか, バグが無いか, 文献に誤りがないか等を確認する意味も含め, 規格化条件\int |\psi_n(x)|^2 \mathrm{d}x=1を満足するか簡単に計算する. 数値積分なので完全に1にはならないが, 概ね1なら問題ない.

数値積分のルーチン
integral(func, x_min, x_max, dx) = sum(func.(x_min:dx:x_max)*dx)
規格化条件の確認
println("re = $re")
println("De = $De")
println("k = $k")
println("a = $a")
println("µ = $µ")
println("n \t<ψ|ψ>")
for n in 0:10
    println(n, "\t", integral(r->ψ(n,r,re,De,a,µ)^2, -10, 50, 0.001))
end
出力
    re = 1.38436
    De = 0.13133356259000006
    k = 0.4101703491541377
    a = 1.249623750214457
    µ = 918.554271
    n 	<ψ|ψ>
    0	1.0000000000000002
    1	1.0
    2	0.9999999999999998
    3	1.0000000000000004
    4	1.0
    5	0.9999999999999997
    6	0.9999999999999999
    7	0.9999999999999997
    8	0.9999999999999989
    9	1.0
    10	0.9999999999999999

可視化

調和振動子などの説明でよく見るグラフを描くには, 波動関数をエネルギーの高さまでシフトし, 縦に縮める必要がある. ここでは0.004\psi(x)+E_nn=0,\cdots,10まで描写した.

可視化
# 軸ラベル, 描写範囲, 凡例位置
plot(xlabel="r", ylabel="V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# ポテンシャル
plot!(-3:0.02:17, r->V(r,re,De,a), label="", lc="#000000", lw=2)

for n in 0:10
    # エネルギー(破線)
    hline!([E(n,De,a,µ)], label="", lc="#000000", ls=:dash)
    # エネルギー(実線)
    plot!(0:0.02:7, r-> E(n,De,a,µ)<V(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#000000", lw=2)
    # 波動関数
    plot!(0:0.02:7, r->0.004*ψ(n,r,re,De,a,µ)+E(n,De,a,µ), label="n=$n", lc=n+1, lw=2)
end

plot!()

調和近似

このポテンシャルの調和近似に相当する調和振動子ポテンシャルとそのハミルトニアンの固有値および固有関数は

V(r) = - D_\mathrm{e} + \frac{1}{2} k (r - r_\mathrm{e})^2\\ \psi_n(r) = A_n H_n(\xi) \exp{\left( -\frac{\xi^2}{2} \right)} \\ E_n = -D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right)

ただし, \omega:=\sqrt{k/m}, \xi=\sqrt{\frac{m\omega}{\hbar}} (r - r_\mathrm{e}), A_n=\sqrt{\frac{1}{n! 2^n} \sqrt{\frac{m\omega}{\pi\hbar}}}であり, H_nエルミート多項式である. また,

a = \sqrt{k/2D_\mathrm{e}}

である. エネルギーをモールポテンシャルと比較すると, 第3項のみの違いであることがわかる. また, 調和近似と言ったものの, むしろ先にバネ定数kを決定して上式でaを決定する方が実用的だったりする.

調和振動子のモジュール
module HO
using SpecialPolynomials

# ポテンシャル
V(r,re,De,k) = -De + 1/2*k*(r-re)^2

# 固有関数
# using SpecialPolynomials
function ψ(r,re,n,k,m)= 1
    ω = sqrt(k/m) # mu = 1 / (1/m1 + 1/m2)
    A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*)))
    ξ = sqrt(m*ω/) * (r-re)
    return A*basis(Hermite, n)(ξ)*exp(-ξ^2/2)
end

# 固有値
function E(n,De,k,m)= 1
    ω = sqrt(k/m)
    return -De +*ω*(n+1/2)
end

end

まずエネルギー準位について比較する.

エネルギーの比較
# 軸ラベル, 描写範囲, 凡例位置
plot(xlabel="r", ylabel="V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# ポテンシャル
plot!(0:0.02:7, r->V(r,re,De,a), label="Morse Potential", lc="#578FC7", lw=2)
plot!(0:0.02:7, r->HO.V(r,re,De,k), label="Harmonic Oscillator", lc="#BC1C5F", lw=2)

# エネルギー
for n in 0:10
    plot!(0:0.02:7, r-> E(n,De,a,µ)<V(r,re,De,a) ? NaN : E(n,De,a,µ), label="", lc="#578FC7")
    plot!(0:0.02:7, r-> HO.E(n,De,k,µ)<HO.V(r,re,De,k) ? NaN : HO.E(n,De,k,µ), label="", lc="#BC1C5F")
    n += 1
end

plot!()

次に波動関数について比較する.

波動関数の比較
# 軸ラベル, 描写範囲, 凡例位置
plot(xlabel="r", ylabel="V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))

# ポテンシャル
plot!(-3:0.02:17, r->V(r,re,De,a), label="", lc="#000000", lw=2)
plot!(-3:0.02:17, r->HO.V(r,re,De,k), label="", lc="#000000", lw=2, ls=:dash)

# 波動関数
for n in 0:10
    plot!(0:0.02:7, r->0.004*HO.ψ(r,re,n,k,µ)+HO.E(n,De,k,µ), label="", lc=n+1, lw=2, ls=:dash)
    plot!(0:0.02:7, r->0.004*ψ(n,r,re,De,a,µ)+E(n,De,a,µ), label="n=$n", lc=n+1, lw=2)
end

plot!()

時間に依存するSchrödinger方程式

時間に依存するSchrödinger方程式

i\hbar \frac{\mathrm{d}}{\mathrm{d}t} \varPsi(x,t) = \hat{H} \varPsi(x,t)

の解は一般に

\varPsi(x,t) = \sum_{n=0}^\infty c_n \psi_n(x) \exp\left(-i \frac{E_n}{\hbar} t \right)

である. これが解であることは単に代入すれば示せる(総和と微分演算子が交換可能である条件を確認せよ). ただし, c_nは初期条件

\varPsi(x,0) = \sum_{n=0}^\infty c_n \psi_n(x)

を定義する係数である.

時間に依存する波動関数
function Ψ(N,t,r,re,De,a,µ) # ← 大文字のΨなので注意
    sum = 0.0 + im*0.0
    for n in 0:N
        sum += c(n) * ψ(n,r,re,De,a,µ) * exp(im*E(n,De,a,µ)*t)
    end
    return sum
end

# 時間発展の描写
function draw(num)
    anim = Animation()
    N = 10
    for i in 0:100
        # 時刻
        t = 10.0*i
        # 描写範囲, タイトルなど
        plt = plot(title="t=$t", xlabel="r", ylabel="V(r)", xlims=(0.2,6.2), ylims=(-0.14,0.01), legend=:bottomright)#, size=(450,400))
        # ポテンシャル
        plot!(0:0.01:17, r->V(r,re,De,a), label="", lc="#000000", lw=2)
        # 実部 Re Ψ(x,t)
        plot!(plt, 0:0.01:7, r->real(Ψ(N,t,r,re,De,a,µ))*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)),  label="Re", lc=1)
        # 虚部 Im Ψ(x,t)
        plot!(plt, 0:0.01:7, r->imag(Ψ(N,t,r,re,De,a,µ))*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)),  label="Im", lc=2)
        # 確率密度 |Ψ(x,t)|^2
        plot!(plt, 0:0.01:7, r->abs(Ψ(N,t,r,re,De,a,µ))^2*0.05 + sum(c.(0:N) .* E.(0:N,De,a,µ)), label="P(x,t)", lc=3, lw=2)
        # コマの追加
        frame(anim, plt)
    end
    gif(anim, "TDSE_$(num).gif", fps = 20)
end

まず, 初期条件が純粋に基底状態のみ(c_0=1, c_n=0 (n>0))ならば実部と虚部が交互に入れ替わり, 確率密度としては時間変化しない. このことは有限差分法を勉強すると容易に理解できる.

c(n) = n==0 ? 1 : 0
draw(1)

次に, 基底状態と第1励起状態を等しい割合で混ぜた初期条件(c_0=c_1=0.5, c_n=0 (n>1))を考えると, 振動らしい振る舞いを見せるようになる.

c(n) = n<2 ? 0.5 : 0
draw(2)

最後に, c_0=\frac{1}{2}, c_1=\frac{1}{4}, c_2=\frac{1}{8}, \cdots, c_n=\frac{1}{2^{n+1}}とすると, 調和振動子とは異なる奇妙な動きを見ることができる.

c(n) = 2.0^(-n-1)
draw(3)

動作環境

versioninfo()
出力
    Julia Version 1.6.2
    Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
    Platform Info:
      OS: Windows (x86_64-w64-mingw32)
      CPU: Intel(R) Core(TM) i7-4650U CPU @ 1.70GHz
      WORD_SIZE: 64
      LIBM: libopenlibm
      LLVM: libLLVM-11.0.1 (ORCJIT, haswell)

参考文献

この記事の元になったJupyter Notebookのデータは下記のリンクにある.

https://gist.github.com/ohno/4a59620ed9e32ebd2ba7ef8e4868a0f5

Juliaで可視化シリーズ

https://zenn.dev/ohno/articles/3101433fbe9231
https://zenn.dev/ohno/articles/870b0c2a0af590
https://zenn.dev/ohno/articles/f849d98a7f58a9
https://zenn.dev/ohno/articles/e1103bc0d58ceb

Discussion

ログインするとコメントできます