「平面波」の作り方(量子力学シミュレーション超入門【第2回】)

ポテンシャルが0のシュレーディンガー方程式の解は、平面波って呼ばれる解になるよ。

平面波の表式(1次元)

波数 $k$、角振動数 $\omega$ としたときの位置 $x$、時刻 $t$ の波動関数 $\Psi(x,t)$ だよ。

\begin{align}
\Psi(x,t) = A e^{i k x – i \omega i } + B e^{i k x + i \omega t}
\end{align}

$k$ と $\omega$ はエネルギー $E$ と

\begin{align}
k = \sqrt{\frac{2mE}{\hbar^2}} \ , \ \omega = \frac{E}{\hbar}
\end{align}

と関係があるから、$k$ と $\omega$ は独立ではなくて

\begin{align}
\omega = \frac{\hbar k^2}{2m}
\end{align}

という分散関係って呼ばれる関係もあるね。

特別な場合の計算結果

電子のエネルギー $E = 0.25, 1.0, 4.0 [{\rm eV}] $ の平面波の実部の様子を可視化するためのデータ出力用プログラムソースだよ。

プログラムソース(C++)


////////////////////////////////////////////////////////////////////
// 自由空間中の粒子
////////////////////////////////////////////////////////////////////
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <stdio.h>

/////////////////////////////////////////////////////////////////
//物理定数
/////////////////////////////////////////////////////////////////
//プランク定数
const double h = 6.6260896 * 1.0E-34;
double hbar = h / (2.0*M_PI);
//電子の質量
const double me = 9.10938215 * 1.0E-31;
//電子ボルト
const double eV = 1.60217733 * 1.0E-19;
/////////////////////////////////////////////////////////////////
//物理系の設定
/////////////////////////////////////////////////////////////////

//電子のエネルギー
double E1 = 0.25 * eV;
double E2 = 1.0 * eV;
double E3 = 4.0 * eV;
//波数
double k1 = sqrt(2.0 * me * E1 / (hbar*hbar));
double k2 = sqrt(2.0 * me * E2 / (hbar*hbar));
double k3 = sqrt(2.0 * me * E3 / (hbar*hbar));
//角振動数
double omega1 = E1 / hbar;
double omega2 = E2 / hbar;
double omega3 = E3 / hbar;

//時間間隔
double dt = 1E-16;
//空間刻み間隔
double dx = 1E-10;

int XN = 400;
int TN = 1000;

double x_min = -20.0;
double x_max = 20.0;

/////////////////////////////////////////////////////////////////

int main() {

	std::cout << 2.0*M_PI / omega2 << " " << 2.0*M_PI / k2 << std::endl;

	//出力ストリームによるファイルオープン
	std::ofstream fout;
	fout.open("wave.txt");
	fout << "#x:位置" << std::endl;
	fout << "#y:波動関数の実部" << std::endl;
	fout << "#showLines: true true true" << std::endl;
	fout << "#showMarkers: false false false" << std::endl;
	fout << "#xrange:" << x_min << " " << x_max << " " << 2 << std::endl;
	fout << "#yrange:" << -1.2 << " " << 1.2 << " " << 0.2 << std::endl;

	//各時刻における計算を行う
	for (int tn = 0; tn < TN; tn++) {
		double t = dt * tn;
		std::cout << tn << std::endl;
		fout << "#coma:" << tn << std::endl;

		for (int ix = 0; ix <= XN; ix++) {
			double x = (x_min + (x_max - x_min) * ix / XN) * dx;
			double psi1 = cos(k1 * x - omega1 * t);
			double psi2 = cos(k2 * x - omega2 * t);
			double psi3 = cos(k3 * x - omega3 * t);
			fout << x / dx << " " << psi1 << " " << psi2 << " " << psi3 << std::endl;
		}
		fout << std::endl;
	}

	fout.close();
}


水素原子に電磁波(古典・直線偏光)を加えたときのシミュレーション結果(ラビ振動)

以前解説したハミルトニアンを用いて、水素原子の基底状態の電子に電磁波(古典・円偏光)を入射したの様子をシミュレーションしたよ。結果を示すね。

角振動数 $\omega_{12} = (E_2-E_1)/\hbar$ の電磁波を入射

次のグラフは、基底状態と第1励起状態のエネルギー準位の差に対応する電磁波を入射した結果だよ。$ \varphi_{100} $ と $ \varphi_{210}, \varphi_{21-1}, \varphi_{21+1} $ の間でラビ振動している様子が分かるね。

直線偏光の場合は、$\varphi_{210}$ にしか遷移しなかったけれども、円偏光を入射すると $ \varphi_{21-1} $ と $ \varphi_{21+1} $ にもそれぞれ25%づつ遷移するね。これらの和は、y軸方向に向いた電気双極子を表すので、円偏光と言っても、実質的には、z軸方向とy軸方向にそれぞれ電気双極子が誘起されるって感じだね。


静磁場(z軸方向)+振動磁場(y軸方向)のスピノールの時間発展(計算結果)

前回示したスピノールの時間発展の表式に基づいて、数値計算を行ったので結果を示すよ。
次のグラフは、時間依存する磁場を $\boldsymbol{B}(t) = (0, B_y\sin(\omega t), B_z)$ で、$ B_y = B_z = 0.1 [{\rm T}]$、 $\omega = 2\omega_L$ として、ルンゲ・クッタ法を用いて、時間発展を計算した結果だよ。初期状態は下向きスピン100%、上向きスピン0%です。

磁気共鳴によってスピンの向きが変化している様子が分かるね。でもラビ振動のように単純な形では無いけどね。

振動磁場の角振動数に対する遷移確率の最大値

次のグラフは、振動磁場の角振動数を $\bar{\omega} = 2\omega_L$ を基準として0から4まで変化させたときの遷移確率の最大値をプロットした結果だよ。

$\omega/\bar{\omega} = 1.0$ 過ぎから $1.2$ 手前までが、最大遷移確率が100%となっているね。次のグラフは、$\omega/\bar{\omega} =1.1$ としたときの、スピノールの時間変化の結果だよ。遷移確率が100%になっている時間があることが分かるね。


スピノールの時間発展(パウリ方程式のスピン依存)

磁場中の電子に対するスピン自由度を含めた波動関数 $\Psi(\boldsymbol{r}, t)$ は、シュレーディンガー方程式から導かれるパウリの方程式と呼ばれる

\begin{align}
i\hbar \frac{\partial \Psi(\boldsymbol{r}, t)}{\partial t} = \left[ \frac{1}{2m_e} (-i\hbar\nabla + e\boldsymbol{A})^2 + V(r) + \frac{e}{m_e} \hat{\boldsymbol{S}} \cdot \boldsymbol{B} \right] \Psi(\boldsymbol{r}, t)
\end{align}

を満たすね。この方程式を満たす波動関数は空間依存部分 $\psi(\boldsymbol{r}, t)$ とスピン依存部分 $\chi(t)$ に分離して、$\Psi(\boldsymbol{r}, t)=\psi(\boldsymbol{r}, t)\chi(t)$ と表わすことで、方程式は変数分離できるね。

\begin{align}
i\hbar \frac{\partial \psi(\boldsymbol{r}, t)}{\partial t} &\, = \left[ \frac{1}{2m_e} (-i\hbar\nabla + e\boldsymbol{A})^2 + V(r) \right] \psi(\boldsymbol{r}, t)\\
i\hbar \frac{d \chi( t )}{d t} &\, = \frac{e}{m_e}\, \hat{\boldsymbol{S}} \cdot \boldsymbol{B}\, \chi(\boldsymbol{r}, t)
\end{align}

このスピン依存部分は、スピン角運動量 $\boldsymbol{S} = (S_x, S_y, S_z)$ をパウリ行列

\begin{align}
S_x = \frac{\hbar}{2} \left( \matrix{ 0 & 1 \cr 1 & 0} \right) \ , \ S_y = \frac{\hbar}{2} \left( \matrix{ 0 & -i \cr 1 & 0} \right)\ , \ S_z = \frac{\hbar}{2} \left( \matrix{ 1 & 0 \cr 0 & -1} \right)
\end{align}

で表したときのスピノールと呼ばれる縦ベクトルで表すことができるね。スピノールを

\begin{align}
\chi( t ) = \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right)
\end{align}

と表しておいて、パウリ方程式に代入して、$\chi_{\uparrow}(t)$ と $\chi_{\downarrow}(t)$ の時間依存性を計算することができるね。

静磁場中のスピノールの時間依存性

まずは、静磁場を $\boldsymbol{B} = (0, 0, B_z)$ として、静磁場中のスピノールの時間依存性を計算してみよう。パウリ方程式のスピン部分は

\begin{align}
i\hbar \frac{d }{d t} \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) &\ = \frac{e\hbar B_z}{2m_e}\, \left( \matrix{ 1 & 0 \cr 0 & -1} \right) \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ = \frac{e\hbar B_z}{2m_e}\,\left( \matrix{ \chi_{\uparrow}(t) \cr – \chi_{\downarrow}(t)} \right)
\end{align}

となるね。これは直ちに解くことができて、

\begin{align}
\chi( t ) = \left( \matrix{\chi_{\uparrow}(0) e^{-i\omega_L t} \cr \chi_{\downarrow}(0) e^{i\omega_L t}} \right) \ , \ \omega_L = \frac{e B_z}{2m_e}
\end{align}

となるね。$ |\chi( t )|^2 = |\chi( 0 )|^2 $ を満たすから、時間が経っても初期状態から変化しないことが分かるね。

振動磁場中のスピノールの時間依存性

次は、時間依存する磁場を $\boldsymbol{B}(t) = (0, 0, B_z\cos(\omega t))$ として、振動磁場中のスピノールの時間依存性を計算してみよう。パウリ方程式のスピン部分は

\begin{align}
i\hbar \frac{d }{d t} \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) &\ = \frac{e\hbar B_z\cos(\omega t)}{2m_e}\,
\left( \matrix{ 1 & 0 \cr 0 & -1} \right) \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ = \frac{e\hbar B_z\cos(\omega t)}{2m_e}\,\left( \matrix{ \chi_{\uparrow}(t) \cr – \chi_{\downarrow}(t)} \right)
\end{align}

となるね。これも直ちに解くことができて、

\begin{align}
\chi( t ) = \left( \matrix{\chi_{\uparrow}(0) e^{-i\frac{\omega_L}{\omega} \sin(\omega t)} \cr \chi_{\downarrow}(0) e^{i\frac{\omega_L}{\omega} \sin(\omega t)}} \right) \ , \
\omega_L = \frac{e B_z}{2m_e}
\end{align}

となるね。振動磁場を与えても $ |\chi( t )|^2 = |\chi( 0 )|^2 $ を満たすから、時間が経っても初期状態から変化しないことが分かるね。

静磁場に垂直方向の振動磁場を加えたときのスピノールの時間依存性

今度は、時間依存する磁場を $\boldsymbol{B}(t) = (0, B_y\sin(\omega t), B_z)$ として、振動磁場中のスピノールの時間依存性を計算してみよう。パウリ方程式のスピン部分は

\begin{align}
i\hbar \frac{d }{d t} \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) &\ =\frac{e\hbar }{2m_e}\left[B_y\sin(\omega t)\left( \matrix{ 0 & -i \cr i & 0} \right) + B_z\left( \matrix{ 1 & 0 \cr 0 & -1} \right)\right]
\left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ =\frac{e\hbar }{2m_e}\left( \matrix{ B_z& -iB_y\sin(\omega t) \cr iB_y\sin(\omega t) & -B_z} \right) \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ = \frac{e\hbar }{2m_e}\,\left( \matrix{ B_z\chi_{\uparrow}(t)-iB_y\sin(\omega t)\chi_{\downarrow}(t) \cr iB_y\sin(\omega t) \chi_{\uparrow}(t) -B_z\chi_{\downarrow}(t)} \right)
\end{align}

と、 $\chi_{\uparrow}(t)$ と $\chi_{\downarrow}(t)$ に関する連立微分方程式が得られるね。解析的に解けるかはわからないけれども、数値的に解くことで、スピノールの時間依存性を計算することができそうだね。異常ゼーマン効果で分かれる2つの準位間のエネルギーに相当する振動磁場(電磁波)の角振動数を与えると、共鳴が起こってスピンの向きが変わるね。これは磁気共鳴現象って呼ばれるよ。これを次回の課題にするね。

円偏光を入射したときのスピノールの時間依存性

最後に、時間依存する磁場を $\boldsymbol{B}(t) = (0, B_0\sin(\omega t), B_0\cos(\omega t))$ として、円偏光を入射したときのスピノールの時間依存性を計算してみよう。パウリ方程式のスピン部分は

\begin{align}
i\hbar \frac{d }{d t} \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) &\ =\frac{e\hbar }{2m_e}\left[B_y\left( \matrix{ 0 & -i \cr i & 0} \right) + B_z\cos(\omega t)\left( \matrix{ 1 & 0 \cr 0 & -1} \right)\right]
\left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ =\frac{e\hbar }{2m_e}\left( \matrix{ B_z\cos(\omega t) & -iB_y \cr iB_y & -B_z\cos(\omega t)} \right) \left( \matrix{\chi_{\uparrow}(t) \cr \chi_{\downarrow}(t)} \right) \\
&\ = \frac{e\hbar }{2m_e}\,\left( \matrix{ B_z\cos(\omega t)\chi_{\uparrow}(t)-iB_y\chi_{\downarrow}(t) \cr iB_y \chi_{\uparrow}(t) -B_z\cos(\omega t)\chi_{\downarrow}(t)} \right)
\end{align}

と、 $\chi_{\uparrow}(t)$ と $\chi_{\downarrow}(t)$ に関する連立微分方程式が得られるね。解析的に解けるかはわからないけれども、数値的に解くことで、スピノールの時間依存性を計算することができそうだね。円偏光でもスピンの向きは変化しそうだね。これも次回の課題にするね。


水素原子に電磁波(古典・円偏光)を加える場合のハミルトニアン

直線偏光の古典電磁波を水素原子に束縛された電子のハミルトニアンは以前示したね。今度は、円偏光の場合を導出するよ。偏光を考慮したベクトルポテンシャルは

\begin{align}
\boldsymbol{A} = \left\{ \matrix{ 0 \cr A_{y0} \cos(kx-\omega t + \phi_y) \cr A_{z0} \cos(kx-\omega t+ \phi_z)} \right.
\end{align}

と与えると、電場と磁場は次のようになるね。

\begin{align}
\boldsymbol{E} = -\frac{\partial \boldsymbol{A}}{\partial t} =\left\{ \matrix{ 0 \cr -\omega A_{y0} \sin(kx-\omega t + \phi_y) \cr -\omega A_{z0} \sin(kx-\omega t+ \phi_z)} \right.
\end{align}
\begin{align}
\boldsymbol{B} = \nabla \times \boldsymbol{A} =\left\{ \matrix{ 0 \cr k A_{z0} \sin(kx-\omega t + \phi_z) \cr -k A_{y0} \sin(kx-\omega t+ \phi_y)} \right.
\end{align}

これを電磁場中の電子のハミルトニアン

\begin{align}
\hat{H} = \hat{H}_0 + \frac{e}{m_e}\,\boldsymbol{A}\cdot \boldsymbol{p} + \frac{e}{m_e} \hat{\boldsymbol{S}}\cdot \boldsymbol{B} + \frac{e^2 \boldsymbol{A}^2}{2m_e}
\end{align}

に代入すると完成だね。

電磁波による状態遷移の数値計算

時間依存性を計算するために、ハミルトニアンの時間依存部分 $\hat{V}(t)$ の空間積分の項を見てみよう。

\begin{align}
\int \varphi_{n’l’m’\sigma’}^* \hat{V}(t)\varphi_{nlm\sigma} dV &\ = \int \varphi_{n’l’m’\sigma’}^* \left[
-i\frac{e}{\hbar} [\hat{H}_0,\boldsymbol{r} ]\cdot\boldsymbol{A} + \frac{e }{m_e} \hat{\boldsymbol{S}}\cdot\boldsymbol{B} + \frac{e^2 \boldsymbol{A}^2}{2m_e} \right] \varphi_{nlm\sigma} dV \\
&\ = \int \varphi_{n’l’m’\sigma’}^* \left[
-ie \frac{ E_{n’} – E_n }{\hbar}\,\boldsymbol{r}\cdot\boldsymbol{A} + \frac{e }{m_e} \hat{\boldsymbol{S}}\cdot\boldsymbol{B} + \frac{e^2 \boldsymbol{A}^2}{2m_e} \right] \varphi_{nlm\sigma} dV
\end{align}

スピンと2次の項を無視して、さらに長波長近似(電磁波の波長が原子サイズより十分に大きいとする近似)を行うと、

\begin{align}
\int \varphi_{n’l’m’}^* \hat{V}(t)\varphi_{nlm} dV &\ \simeq -ie \int \varphi_{n’l’m’}^* \left[
\frac{ E_{n’} – E_n }{\hbar}\,(yA_y +zA_z) \right] \varphi_{nlm} dV \\
&\ = \frac{-ie}{2}\,\frac{ E_{n’} – E_n }{\hbar} \left[A_{y0}\cos(\omega t – \phi_y) \int \varphi_{n’l’m’}^* y \varphi_{nlm} dV + A_{z0}\cos(\omega t – \phi_z) \int \varphi_{n’l’m’}^* z \varphi_{nlm} dV\right]
\end{align}

となるね。入射波を円偏光とするには、 $A_{y0} = A_{z0} = A_{0}$ かつ $\phi_y = 0 , \phi_z = \pm \pi/2 $ とすれば良いので、これを代入すると、

\begin{align}
\int \varphi_{n’l’m’}^* \hat{V}(t)\varphi_{nlm} dV &\ \simeq -ie \int \varphi_{n’l’m’}^* \left[
\frac{ E_{n’} – E_n }{\hbar}\,(yA_y +zA_z) \right] \varphi_{nlm} dV \\
&\ = \frac{-ieA_0}{2}\,\frac{ E_{n’} – E_n }{\hbar} \left[\cos(\omega t) \int \varphi_{n’l’m’}^* y \varphi_{nlm} dV \pm \sin(\omega t) \int \varphi_{n’l’m’}^* z \varphi_{nlm} dV\right]
\end{align}

となるね。このハミルトニアンを元に次回は水素原子の基底状態にいる電子に円偏光電磁波を与えてみるよ。


水素原子に電磁波(古典・直線偏光)を加えたときのシミュレーション結果(ラビ振動)

以前解説したハミルトニアンを用いて、水素原子の基底状態の電子に電磁波(古典・直線偏光)を入射したの様子をシミュレーションしたよ。結果を示すね。

角振動数 $\omega_{12} = (E_2-E_1)/\hbar$ の電磁波を入射

次のグラフは、基底状態と第1励起状態のエネルギー準位の差に対応する電磁波を入射した結果だよ。$ \varphi_{100} $ と $ \varphi_{210} $ の間でラビ振動している様子が分かるね。

角振動数 $\omega_{12} = (E_2-E_1)/\hbar$ と $\omega_{23} = (E_3-E_2)/\hbar$の電磁波を入射

次のグラフは、基底状態と第1励起状態のエネルギー準位の差に対応する電磁波と、第1励起状態と第2励起状態ののエネルギー準位の差に対応する電磁波を同時に入射した結果だよ。$ \varphi_{100} $ と $ \varphi_{210} $ に加えて、$ \varphi_{100} $ から直接遷移することができない $ \varphi_{300} $ と $ \varphi_{320} $ にも励起しているね。このように、直接遷移が許されない状態が加わっても、存在確率が周期的に振動するのはちょっと意外だね。

次回は、円偏光の電磁波を入射するよ。


水素原子に静電場を急激に加えたときのシミュレーション結果

前回解説した水素原子に静電場を急激に加えたときの緩和時間シミュレーションの結果を示すよ。次の図は初期状態として $\varphi_{200}$ 100%の状態に $ E_z = 10^{9}[\rm{V/m}]$ の静電場を急に加えたときの $\varphi_{200}$ と $\varphi_{210}$ の存在確率の時間経過だよ。$\varphi_{200}$ 100%の状態と $\varphi_{210}$ 100%の状態が一定の周期で交互に現れるね。

考察:単振動的な運動をする理由

$\varphi_{200}$ と $\varphi_{210}$ は静電場中ではエネルギーが高い一方で、その50%づつの混合状態が一番エネルギーが低いのだよね。つまり、下の図で示したとおり、エネルギーの高い初期状態 $\varphi_{200}$ 100% からスタートして、エネルギー低い方に状態が変化して行くけれども、一番低いところで止まらずに反対の $\varphi_{210}$ 100% の状態まで変化していっているね。もし、エネルギーが散逸するメカニズムがあれば、最低エネルギーに落ち着くよね。きっと。