水素原子の外場による光電効果の数値計算方法

水素原子の電子に直線偏光の電磁波（古典）を入射したときの状態遷移は以前解説したね。今回は、入射した電磁波の振動数を高めることで生じる「光電効果」（束縛状態の電子が外場によって非束縛状態へ遷移する効果）をシミュレーションするのに必要な表式を導出するよ。スピンと外場の２次の項を無視したハミルトニアンは、

\begin{align}
\hat{H} = \hat{H}_0 + \frac{e}{m_e}\,\boldsymbol{A}\cdot \hat{\boldsymbol{p}} = -\frac{\hbar^2}{2m_e}\nabla^2 + V(r) + \frac{e}{im_e} \,\boldsymbol{A}\cdot\nabla
\end{align}

となるね。このハミルトニアンを用いたシュレディンガー方程式

\begin{align}
i\hbar\frac{\partial}{\partial t} \psi(\boldsymbol{r}, t) = \hat{H} \psi(\boldsymbol{r}, t)
\end{align}

だね。ここから出発するよ。今回は、水素原子核に束縛された状態からとき放たれた状態も考慮するため、波動関数は波数空間で展開した

\begin{align}
\psi(\boldsymbol{r}, t) = \frac{1}{\sqrt{V}}\sum\limits_{n_x, n_y, n_z} a_{n_xn_yn_z}(t)\, e^{i\boldsymbol{k}\cdot\boldsymbol{r}}
\end{align}

で表すよ。ただし、$ \boldsymbol{k} = \frac{2\pi}{L}(n_x, n_y, n_z) $ で、$L$ は空間サイズで $V =L^3$、 $n_x, n_y, n_z$ は整数だよ。これをシュレディンガー方程式に代入して、両辺に $\frac{1}{\sqrt{V}}e^{-i\boldsymbol{k}’\cdot\boldsymbol{r}}$ を掛けて、全空間で積分すると

\begin{align}
i\hbar\frac{d a_{n_x’n_y’n_z’}(t)}{d t} = \left[ \frac{\hbar^2k’^2}{2m_e} + \frac{e}{m_e} \,\boldsymbol{A}\cdot \boldsymbol{k}’ \right] a_{n_x’n_y’n_z’}(t) + \frac{1}{V}\sum\limits_{n_x, n_y, n_z} \left[ \int e^{-i\boldsymbol{k}’\cdot\boldsymbol{r}} V(r) e^{i\boldsymbol{k}\cdot\boldsymbol{r}} dV \right] a_{n_xn_yn_z}(t)
\end{align}

となるね。さらにポテンシャル項もついでに

\begin{align}
V(r) = \frac{1}{\sqrt{V}}\sum\limits_{n”_x, n”_y, n”_z} v_{n”_xn”_yn”_z}\, e^{i\boldsymbol{k}”\cdot\boldsymbol{r}}
\end{align}

と展開して代入すると、空間積分を実行することができて、次のような形になるね。

\begin{align}
i\hbar\frac{d a_{n_x’n_y’n_z’}(t)}{d t} = \left[ \frac{\hbar^2k’^2}{2m_e} + \frac{e}{im_e} \,\boldsymbol{A}\cdot \boldsymbol{k}’ \right] a_{n_x’n_y’n_z’}(t) + \frac{1}{\sqrt{V}}\sum\limits_{n_x, n_y, n_z}\,v_{n_x’-n_x,n_y’-n_y,n_z’-n_z} a_{n_xn_yn_z}(t)
\end{align}

さらに、クーロンポテンシャル $ V(r) = -e^2/4\pi\epsilon r$のように $1/r$ の場合には、展開係数 $ v_{n”_xn”_yn”_z} $ は解析的に導出することができて、

\begin{align}
v_{n”_xn”_yn”_z} = \frac{1}{\sqrt{V}} \int V(r)e^{-i\boldsymbol{k}”\cdot\boldsymbol{r}} dV = -\frac{e^2}{4\pi\epsilon_0} \, \frac{1}{k”^2}
\end{align}

で与えられるので、最終的には

\begin{align}
i\hbar\frac{d a_{n_x’n_y’n_z’}(t)}{d t} = \left[ \frac{\hbar^2k’^2}{2m_e} + \frac{e}{im_e} \,\boldsymbol{A}\cdot \boldsymbol{k}’ \right] a_{n_x’n_y’n_z’}(t) -\frac{e^2}{4\pi\epsilon_0}\, \frac{1}{\sqrt{V}}\sum\limits_{n_x, n_y, n_z}\,\frac{1}{k_{n_x’-n_x,n_y’-n_y,n_z’-n_z}^2 }a_{n_xn_yn_z}(t)
\end{align}

となるね。これで、$a_{n_xn_yn_z}$ に関する連立常微分方程式になるね。あとは、外場を与えるベクトルポテンシャル $\boldsymbol{A}$ を計算対象の系に合わせて設定して、ルンゲ・クッタ法を用いて時間発展を計算することができるね。ちなみに、$a_{n_xn_yn_z}$ の初期状態は、

\begin{align}
a_{n_xn_yn_z}(0) = \frac{1}{\sqrt{V}}\sum\limits_{n_x, n_y, n_z} \psi(\boldsymbol{r}, 0) \, e^{-i\boldsymbol{k}\cdot\boldsymbol{r}}
\end{align}

で計算することができるので、例えば、$\psi(\boldsymbol{r}, 0) = \varphi_{100}$ と与えることで、初期状態を基底状態とした場合の計算を行うことができるよ。次回は結果を示すよ。