🎻

Juliaで可視化：調和振動子

2022/01/02に公開

Julia

Julia言語を用いて調和振動子のSchrödinger方程式の解を可視化する.

パッケージ

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

# using Pkg
# Pkg.add("Plots")
# Pkg.add("SpecialPolynomials")
using Plots
using SpecialPolynomials

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

1次元の調和振動子ポテンシャル

V(x) = \frac{1}{2} k x^2

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

\hat{H} = \left[ - \frac{\hbar^2}{2m} \frac{\mathrm{d}^2}{\mathrm{d}x ^2} + \frac{1}{2} m \omega^2 x^2 \right]

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

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

の解は

\psi_n(x) = A_n H_n(\xi) \exp{\left( -\frac{\xi^2}{2} \right)} \\ E_n = \hbar \omega \left( n + \frac{1}{2} \right)

である. ただし, $\omega:=\sqrt{k/m}$ , $\xi=\sqrt{\frac{m\omega}{\hbar}}x$ , $A_n=\sqrt{\frac{1}{n! 2^n} \sqrt{\frac{m\omega}{\pi\hbar}}}$ であり, $H_n$ はエルミート多項式である. エルミート多項式の一般形は微分を含み, 関数として宣言するのが面倒だが, SpecialPolynomials.jlで利用可能である.

ポテンシャル

V(x) = 1/2*k*x^2

固有関数

# using SpecialPolynomials
function ψ(x,n,k,m)
    ℏ = 1
    ω = sqrt(k/m)
    A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ)))
    ξ = sqrt(m*ω/ℏ) * x
    return A*basis(Hermite, n)(ξ)*exp(-ξ^2/2)
end

固有値

function E(n,k,m)
    ℏ = 1
    ω = sqrt(k/m)
    return ℏ*ω*(n+1/2)
end

規格化条件

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

数値積分のルーチン

function integral(func, min, max, dx)
    sum = 0.0
    for x in min:dx:max
        sum += func(x)*dx
    end
    return sum
end

規格化条件の確認

k = 1.0
m = 3.0

# println("k = $k")
# println("m = $m")
println("n \t<ψ|ψ>")
for n in 1:10
    println(n, "\t", integral(x->ψ(x,n,k,m)^2, -10, 10, 0.1))
end

出力

n 	<ψ|ψ>
1	1.0
2	0.9999999999999999
3	1.0
4	1.0
5	0.9999999999999999
6	1.0000000000000002
7	1.0000000000000009
8	1.0
9	1.0000000000000007
10	1.0000000000000007

ビリアル定理

ビリアル定理も検証に使いやすい. 調和振動子では $\langle T \rangle = \langle V \rangle$ が成り立つため, $E = \langle T \rangle + \langle V \rangle = 2\langle V \rangle$ が成り立つ.（クーロンポテンシャルでは $2\langle T \rangle = -\langle V \rangle$ なので注意）

ビリアル定理の確認

# println("k = $k")
# println("m = $m")
println("n \tE \t\t\t2<V>")
for n in 1:10
    println(n, "\t", E(n,k,m), "\t", integral(x->V(x)*ψ(x,n,k,m)^2, -10, 10, 0.1)*2 )
end

出力

n 	E 			2<V>
1	0.8660254037844386	0.8660254037844384
2	1.4433756729740643	1.4433756729740643
3	2.02072594216369	2.0207259421636907
4	2.598076211353316	2.5980762113533142
5	3.1754264805429413	3.175426480542942
6	3.7527767497325675	3.752776749732567
7	4.330127018922193	4.330127018922195
8	4.907477288111819	4.907477288111821
9	5.484827557301444	5.484827557301447
10	6.06217782649107	6.062177826491078

可視化

調和振動子ポテンシャルの説明でよく見るグラフを描くには, 波動関数をエネルギーの高さまでシフトし, 縦に縮める必要がある. ここでは $0.2\psi(x)+E_n$ を $n=0,1,2$ まで描写した.

# パラメータ
k = 1.0
m = 3.0

# 軸ラベル, 描写範囲, 凡例位置
plot(xlabel="x", ylabel="V(x)", xlims=(-3.5,3.5), ylims=(0,1.8), legend=:bottomright)#, size=(450,400))

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

for i in 0:2
    # エネルギー（破線）
    hline!([E(i,k,m)], label="", lc="#000000", ls=:dash)
    # エネルギー（実線）
    plot!([-sqrt(2/k*E(i,k,m)), sqrt(2/k*E(i,k,m))], fill(E(i,k,m),2), label="", lc="#000000", lw=2)
    # 波動関数
    plot!(-3.5:0.02:3.5, x->0.2*ψ(x,i,k,m)+E(i,k,m), label="n=$i", lc=i+1, lw=2)
end

savefig("TISE.svg")

plot!()

時間に依存する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 Ψ(x,t,N,k,m)
    sum = 0.0 + im*0.0
    for n in 0:N
        sum += c(n) * ψ(x,n,k,m) * exp(im*E(n,k,m)*t)
    end
    return sum
end

時間発展の描写

function draw(num)
    anim = Animation()
    for i in 0:100
        # 時刻
        t = 0.5*i
        # 描写範囲, タイトルなど
        plt = plot(xlims=(-3,3), ylims=(-0.5,0.8), title="t=$t", xlalbel="x")
        # ポテンシャル
        plot!(-3:0.01:3, x->V(x), label="V(x)", lc="#000000")
        # 実部 Re Ψ(x,t)
        plot!(plt, -5:0.02:5, x->real(Ψ(x,t,15,k,m)),  label="Re", lc=1)
        # 虚部 Im Ψ(x,t)
        plot!(plt, -5:0.02:5, x->imag(Ψ(x,t,15,k,m)),  label="Im", lc=2)
        # 確率密度 |Ψ(x,t)|^2
        plot!(plt, -5:0.02:5, x->abs(Ψ(x,t,15,k,m))^2, 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のデータは下記のリンクにある.