EbiMaru’s diary

ラマン分光とか機械学習とか研究している京大生のブログ

調和振動子の古典・量子

そうだ.量子論で一次元調和振動子を解こう

 というわけで確率密度関数を求めます.アニメはpyてょnしか書けないので,それで作ります.

古典力学

調和振動子ハミルトニアンは近所の佐藤少年でも知ってる

\displaystyle{
H = \frac{p^2}{2m} + \frac{m \omega^2}{2}x^2.
}

(m:質量,\omega:振動数).解き方は幼稚園児でもわかるので省略.

\displaystyle{
x(t) = x(0) \cos \omega t + \frac{p(0)}{m \omega} \sin \omega t
}
\displaystyle{
p(t) = -m \omega x(0) \sin \omega t + p(0)\cos \omega t
}

量子力学

面倒くさいなあ.まずは状態ケットの時間変化を見るか(Schrödinger描像にしときますわ.あと表記の都合上,波動関数が時間に依存するかどうかは引数の中身で区別しています).

\displaystyle{
H \left| n \right> = \left( \frac{1}{2} +n  \right) \hbar \omega \left| n \right>
}

だから,

\displaystyle{
\mathcal{U(t)}  \left| n \right> =  \exp \left[ -i  \frac{H t}{\hbar}\right] \left| n \right>= \exp \left[ -i  \left( \frac{1}{2} + n \right) \omega t \right] \left| n \right>
}

(\mathcal{U}(t):時間発展演算子)です.エネルギー固有関数が求まれば,波動関数も求まることがわかります. \sigma = \sqrt{\hbar/ m \omega}を長さの規格化定数とします.それらを求める前にどのくらいの長さか考えてみましょう.水素分子のばね定数が510 N/m程度

http://www.comp.tmu.ac.jp/qc/toudai/2007-11-9.pdf, p. 6.

mは換算質量(m = m_\mathrm{hydrogen}/2)であることを考慮して水素分子の振動数は

\displaystyle{
\omega \approx \sqrt{510/(1.674 \times 10^{-27} / 2)} = 7.81 \times 10^{14} \hspace{5pt} \mathrm{s}^{-1}
}

程度です.また,\sigma

\displaystyle{
\sigma = \sqrt{\hbar/\sqrt{km}}  \approx  \sqrt{1.055 \times 10^{-34} / \sqrt{1.674\times10^{-27} \times 510 /2}} 
}
\displaystyle{
= 1.27 \times 10^{-11} \hspace{5pt} \mathrm{m} 
}

で,大体0.1 Åだから,そんなものかあといった感じです.
 エネルギー固有関数と波動関数を求めていきましょう.u_n(x')を位置表示のエネルギー固有関数,v_n(p')を運動量表示のエネルギー固有関数とします.また,\psi_n(x',t)を位置表示の波動関数\phi_n(p',t)を運動量表示の波動関数とします.

\displaystyle{
\sqrt{\frac{m \omega}{2 \hbar}} \left< x' \left| a \right|0 \right> = \sqrt{\frac{m \omega}{2 \hbar}} \left< x' \left| \left( x + \frac{ip}{m \omega} \right) \right|0 \right> = 0
}

( a は消滅演算子)から

\displaystyle{
\left( x' + \sigma^2 \frac{d}{dx'} \right) \left< x' | 0 \right> = 0
}

この微分方程式を解くと基底状態のエネルギー固有関数が求められ

\displaystyle{
u_0(x') = \left< x' | 0 \right> = \left( \frac{1}{\pi \sigma^2 } \right)^{1/4} \exp \left[ - \frac{1}{2} \left( \frac{x'}{\sigma}\right)^2 \right]
}

です. また

\displaystyle{
   \left< x' | 1 \right> = \left< x' | a^\dagger| 0 \right>
}

から

\displaystyle{
  u_1(x') = \left< x' | 1 \right> =\left( \frac{1}{\sqrt{2} \sigma } \right) \left( x' - \sigma^2 \frac{d}{dx'} \right) \left< x' | 0 \right> = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} \frac{x'}{\sigma} \exp \left[ - \frac{1}{2} \left( \frac{x'}{\sigma} \right)^2 \right]
}

まあ,エルミート多項式使ったほうが早いんですけどね.さて,波動関数がもとまったので粒子の存在確率\rho_xを求めます.初期状態が重ね合わせ状態じゃないと振動しないので

\displaystyle{
   \left| \alpha;t=0 \right> = \frac{1}{\sqrt{2}} \left| 0 \right> + \frac{1}{\sqrt{2}} \left| 1 \right>
}

とします.ゴリゴリ計算して,

\displaystyle{
   \rho_x(x',t) = \left| \frac{1}{\sqrt{2}} \psi_0(x',t) + \frac{1}{\sqrt{2}} \psi_1 (x',t) \right|^2 
}
\displaystyle{
=  \frac{1}{ \sigma \sqrt{\pi}} \ \exp \left[ - \left( \frac{x'}{\sigma} \right)^2 \right] \left[  \frac{1}{2}  +  \left( \frac{x'}{\sigma} \right)^2 + \sqrt{2} \frac{x'}{\sigma} \cos \omega t \right]
}

です.確かにx'積分すれば,1になりますね.次は運動量の確率密度\rho_p(p,t)を求めましょう.運動量表示のエネルギー固有関数v_0(p')v_1(p')を求めます.簡単に計算できて(←言ってみたかっただけ)

\displaystyle{
v_0 (p') = \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \exp \left( \frac{- ip' x'}{\hbar} \right) \psi_0 (x')
}
\displaystyle{
= \sqrt{ \frac{\sigma}{ \hbar \sqrt{\pi}}} \exp \left[ - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right]
}

v_1(p')は少し面倒ですが,

\displaystyle{
v_1 (p') = \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \exp \left( \frac{- ip' x'}{\hbar} \right) \psi_1 (x')
}
\displaystyle{
= \frac{1}{\sqrt{2 \pi \hbar}} \int dx' \left[  \exp \left\{ - \frac{1}{2} \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right)^2 - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right\}  \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right) - \frac{ip'\sigma}{\hbar} \exp \left\{ - \frac{1}{2} \left( \frac{x'}{\sigma} + \frac{i p' \sigma}{\hbar} \right)^2 - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right\} \right]
}
\displaystyle{
= \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} \frac{- i p' \sigma}{\hbar}\exp \left[ - \frac{1}{2} \left( \frac{p' \sigma}{\hbar} \right)^2 \right]
}

なので,運動量表示の波動関数も求められます.確率密度関数\rho_p(p',t)

\displaystyle{
   \rho_p(p',t) = \left| \frac{1}{\sqrt{2}} \phi_0(p',t) + \frac{1}{\sqrt{2}} \phi_1 (p',t) \right|^2 
}
\displaystyle{
  = \frac{\sigma}{\hbar \sqrt{\pi}} \exp \left[ - \left( \frac{p' \sigma}{\hbar} \right)^2 \right] \left[ \frac{1}{2} + \left( \frac{p'\sigma}{\hbar} \right)^2 - \sqrt{2} \frac{p'\sigma}{\hbar}  \sin \omega t  \right] 
}

割と対称的できれいですわね(もっと簡単に求める方法があったら教えてください><) さて,アニメーションだ. \hbar=\sigma=1,\omega = 2 \piとします.結果は以下のとおりです
f:id:EbiMaru:20190721204339g:plain
ミュー空間(と言っていいのかな...)でのカラーマップはこんな感じ
f:id:EbiMaru:20190721204351g:plain
なんかあってる気がしないな笑.期待値はちゃんと楕円軌道になっていると信じたい.
 この調子で高エネルギーの初期条件の運動も求めていきます.さて,ここであらかじめ作っておいた準位0~3の固有関数を用意します(笑). 計算の都合上,X' = x/\sigmaP' = p \sigma /\hbar として,

\displaystyle{
u_0(X') = \left( \frac{1}{\pi \sigma^2 } \right)^{1/4} \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_1(X') = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} X' \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_2(X') = \left( \frac{4}{\pi \sigma^2} \right)^{1/4} \left[ X'^2 - \frac{1}{2} \right] \exp \left[ - \frac{1}{2} X'^2 \right]
}
\displaystyle{
  u_3(X') = \left( \frac{4}{9 \pi \sigma^2} \right)^{1/4} \left[ X'^3 - \frac{3}{2} X' \right] \exp \left[ - \frac{1}{2} X'^2 \right]
}

積分

\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] = \sqrt{2 \pi} \exp \left[ - \frac{1}{2} P'^2\right] 
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X' = i \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left( - P' \right)
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X'^2 = \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left( - P'^2 + 1 \right)
}
\displaystyle{
\int dX' \exp \left[ - \frac{1}{2} \{ X' + iP' \}^2 - \frac{1}{2} P'^2 \right] X'^3 =i \sqrt{2 \pi}  \exp \left[ - \frac{1}{2} P'^2\right] \left(  P'^3 - 3P'  \right)
}

を用いてこれに対する運動量表示の固有関数を求め

\displaystyle{
v_0 (P') = \sqrt{ \frac{\sigma}{ \hbar \sqrt{\pi}}} \exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_1(P') = i \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} (-P')\exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_2(P') = \sqrt{ \frac{2 \sigma}{ \hbar \sqrt{\pi}}} \left( - P'^2 + \frac{1}{2}  \right) \exp \left[ - \frac{1}{2} P'^2 \right]
}
\displaystyle{
v_3(P') = i \sqrt{ \frac{2 \sigma}{ 3 \hbar \sqrt{\pi}}} \left(  P'^3  -\frac{3}{2} P' \right) \exp \left[ - \frac{1}{2} P'^2 \right]
}

さて,同じ重みで重ね合わせされた,準位が1つだけ違う状態の重ね合わせの初期状態での確率密度関数

\displaystyle{
\rho_{X}^{(12)} (X',t) =  \frac{2}{ \sigma \sqrt{\pi}} \ \exp \left(- X'^2 \right) \left[  X'^2 +  \left(  X'^2 - \frac{1}{2} \right)^2 + 2X' \left(  X'^2 - \frac{1}{2} \right) \cos \omega t \right]
}
\displaystyle{
\rho_{X}^{(23)} (X',t) =  \frac{2}{ \sigma \sqrt{\pi}} \ \exp \left(- X'^2 \right) \left[ \left(  X'^2 - \frac{1}{2} \right)^2 + \frac{1}{3} X'^2 \left( X'^2 -\frac{3}{2} \right)^2 +  \frac{2}{\sqrt{3}}X' \left(  X'^2 - \frac{1}{2} \right) \left( X'^2 - \frac{3}{2} \right) \cos \omega t \right]
}
\displaystyle{
\rho_{P}^{(12)} (P',t) =  \frac{2 \sigma}{\hbar \sqrt{\pi}} \exp \left( - P'^2 \right) \left[ P'^2 + \left( P'^2 - \frac{1}{2}\right)^2 - P' \left( P'^2 - \frac{1}{2} \right)^2 \sin \omega t \right]
}
\displaystyle{
\rho_{P}^{(23)} (P',t) =  \frac{2 \sigma}{\hbar \sqrt{\pi}} \exp \left( - P'^2 \right) \left[ \left(  P'^2 - \frac{1}{2} \right)^2 + \frac{1}{3} P'^2 \left( P'^2 -\frac{3}{2} \right)^2 -  \frac{2}{\sqrt{3}}P' \left(  P'^2 - \frac{1}{2} \right) \left( P'^2 - \frac{3}{2} \right) \sin \omega t \right]
}

アニメにするとこんな感じです.

f:id:EbiMaru:20190721204410g:plain
確率密度関数(初期状態:1 and 2)

f:id:EbiMaru:20190721204428g:plain
確率密度関数のカラーマップ(初期状態:1 and 2)

f:id:EbiMaru:20190721204446g:plain
確率密度関数(初期状態:2 and 3)

f:id:EbiMaru:20190721204457g:plain
確率密度関数のカラーマップ(初期状態:2 and 3)

やはり初期エネルギーが大きいと,節が増えていくんですね.
 次に,期待値と偏差を考えていきます.

\displaystyle{
 a = \frac{1}{\sqrt{2}} \left( X + iP \right), \quad a^\dagger = \frac{1}{\sqrt{2}} \left( X - iP \right)
}

で,

\displaystyle{
 X = \frac{1}{\sqrt{2}} \left(  a^\dagger + a  \right),\quad  P = \frac{i}{\sqrt{2}} \left( a^\dagger - a \right)
}

です,初期状態を\left| \alpha;t=0 \right> = \left( \left| n \right> + \left| n+1 \right> \right)/\sqrt{2}とすると

\displaystyle{
 \left| \alpha;t \right> = \frac{1}{\sqrt{2}}\left( \exp \left[ - \frac{1}{2} (n + 1) \omega t \right] \left| n \right> + \exp \left[ - \frac{1}{2} (n + 2) \omega t \right] \left| n+1 \right> \right)
}

から,

\displaystyle{
 \left< X \right> = \sqrt{\frac{n+1}{2}} \cos \omega t , \quad \left< P \right> = -\sqrt{\frac{n+1}{2}} \sin \omega t
}

また,

\displaystyle{
 X^2 = \frac{1}{2} \left( a + {a^\dagger}^2 + a a^\dagger + a^\dagger a \right), \quad P^2 = \frac{1}{2} \left( a a^\dagger + a^\dagger a - a^2 + {a^\dagger}^2 \right)
}

なので,

\displaystyle{
 \left< X^2 \right>= n + 1 ,  \quad \left< P^2 \right> =  n + 1
}

となり,偏差は

\displaystyle{
 \Delta X = \frac{1}{2}\sqrt{(n+1) (3 - \cos 2 \omega t )}, \quad \Delta P =  \frac{1}{2}\sqrt{(n+1) (3 + \cos 2 \omega t )}
}

です.xの広がりは\sigma \propto 1/\sqrt{mk} に,pの広がりは\sigma^{-1}\propto\sqrt{mk}に比例し,いずれも\hbarに比例するので,古典極限は案の定,\hbar \rightarrow 0であり,ばね定数に対する質量の比が小さいほど位置が,大きいほど運動量が不確かになることが分かります.