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

ポテンシャルが0のシュレーディンガー方程式の解である平面波は前回示したね。この平面波を重ね合わせることで、任意の波形を作ることができるよ。今回は、波数分布をガウス分布とした波束の作り方を解説するよ。波束は異なる波数(速度)の平面波の重ね合わせで作ることができるよ。中心波数 $k_0$ としたときの表式は次のとおりだよ。

\begin{align}
\psi(x,t)= \frac{1}{ \sqrt{2\sigma \pi L}}\int dk\, \exp\left[ ik(x-x_0) -i\omega t-\left(\frac{ k-k_0 }{2\sigma} \right)^2 \right]
\end{align}

$k_0$ は波束が進む速度 $v$ と、 $v_0 =p_0 / m_e = \hbar k_0 / m_e$ の関係があるよ。

波束の運動の計算結果

波束の中心エネルギー: $E_0 = 10[{\rm eV}]$
空間スケール: $10^{-11}[{\rm m}]$ //横軸の値
時間スケール: $10^{-16}[{\rm s}]$ //動画1コマの時間間隔

波束の中心エネルギー: $E_0 = 0[{\rm eV}]$
空間スケール: $10^{-11}[{\rm m}]$ //横軸の値
時間スケール: $4\times 10^{-15}[{\rm s}]$ //動画1コマの時間間隔

ちなみに、どんな波束も時間とともに広がっていくよ。その理由は分散関係 $\omega$ が $k$ に比例しないからだよ。

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


////////////////////////////////////////////////////////////////////
// 【第3回】「波束」の作り方
////////////////////////////////////////////////////////////////////
#define _USE_MATH_DEFINES
#include <iostream>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <string>
#include <complex>

/////////////////////////////////////////////////////////////////
//物理定数
/////////////////////////////////////////////////////////////////
//光速
const double c = 2.99792458E+8;
//真空の透磁率
const double mu0 = 4.0*M_PI*1.0E-7;
//真空の誘電率
const double epsilon0 = 1.0 / (4.0*M_PI*c*c)*1.0E+7;
//プランク定数
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;
//複素数
const std::complex<double> I = std::complex<double>(0.0, 1.0);
/////////////////////////////////////////////////////////////////
//物理系の設定
/////////////////////////////////////////////////////////////////
//空間分割数
const int Nx = 500;
//空間分割サイズ
double dx = 1.0E-10;
//重ね合わせの数
const int N = 200;
//パルスの幅
double delta_x = 1.0E-8;
double sigma = 2.0*sqrt(2.0*log(2.0)) / (delta_x);
//波数の間隔
double dk = 30.0 / (delta_x * double(N + 1));
//電子波のエネルギー
const double E0 = 10.0 * eV;
//波数の中心
double k0 = sqrt(2.0 * me * E0 / pow(hbar, 2));
//角振動数の中心
double omega0 = hbar / (2.0*me) * pow(k0, 2);

//計算時間の幅
const int ts = 0, te = 300;
//時間間隔
double dt = 1.0 * 1.0E-16;
/////////////////////////////////////////////////////////////////
const int precision_N = 4;

int main() {
	//出力ストリームによるファイルオープン
	std::ofstream fout;
	fout.open("wave.txt");
	fout << "#x:位置" << std::endl;
	fout << "#y:確率振幅" << std::endl;
	fout << "#legend: 実部 虚部 絶対値" << std::endl;
	fout << "#showLines: true true true" << std::endl;
	fout << "#showMarkers: false false false" << std::endl;
	fout << "#xrange:" << -Nx / 2 << " " << Nx / 2 << " " << Nx / 10 << std::endl;
	fout << "#yrange:" << -0.30 << " " << 0.30 << " " << 0.1 << std::endl;

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

		for (int nx = 0; nx <= Nx; nx++) {
			double x = dx * (nx - Nx/2);
			std::complex<double> Psi = std::complex<double>(0.0, 0.0);
			for (int jz = 0; jz <= N; jz++) {
				double k = (k0 + dk * double(jz - N / 2));
				double omega = hbar / (2.0*me) * pow(k, 2);
				Psi += exp(I*(k*x - omega * t_real)) * exp(-1.0 / 2.0 * pow((k - k0) / sigma, 2));
			}
			Psi = Psi / double(N);
			fout << std::setprecision(precision_N);
			fout << x / dx << " " << Psi.real() << " " << Psi.imag() << " " << abs(Psi) << std::endl;
		}
		fout << std::endl;
	}
	fout.close();
}


「平面波」の作り方(量子力学シミュレーション超入門【第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%になっている時間があることが分かるね。


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

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

角振動数 $\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% の状態まで変化していっているね。もし、エネルギーが散逸するメカニズムがあれば、最低エネルギーに落ち着くよね。きっと。


スピン―軌道相互作用を考慮した水素原子に磁場1テスラ加えたときのエネルギー準位

前回解説した内容を踏まえて、外部から加える磁場を1テスラとしたときの主量子数1から3までのエネルギー準位の計算結果を示すよ。水素原子の場合、外部磁場が1テスラでスピン―軌道相互作用の効果とゼーマン効果がどちらも $10^{-5}[{\rm eV}]$ 程度となるよ。

主量子数1(2状態)

主量子数2(8状態)

主量子数3(18状態)の上半分

主量子数3(18状態)の下半分