Pythonで無限井戸型ポテンシャルのシュレディンガー方程式を解く

はじめに

波動関数をベクトル、演算子を行列でそれぞれ表すと、シュレディンガー方程式は行列の固有値・固有ベクトル問題に帰着させることができます。この事実を活用することで、数値計算でシュレディンガー方程式を解けるようになります。本記事では、よく知られた一次元無限井戸型ポテンシャルを題材に、Pythonを用いた数値解法を紹介します。

シュレディンガー方程式を行列・ベクトルで表現する

一次元無限井戸型ポテンシャルにおける(定常状態の)シュレディンガー方程式は下記の通りです。

$$(−ℏ22md2dx2+V(x))ψ(x)=ϵψ(x) \left(-\frac{\hbar^2}{2m}\frac{\mathrm{d}^2}{\mathrm{d}x^2}+V(x)\right)\psi(x) = \epsilon \psi(x)$$

$$V(x)={0(∣x∣≤L)∞(∣x∣>L) V(x)=\begin{cases} 0 & (|x|\le L)\\ \infty & (|x|>L) \end{cases}$$

この方程式をそのままPython等プログラミング言語は困難です。そこで、波動関数をベクトル、演算子を行列で表現することで離散化します。離散化する方法は基底の選び方によって様々ですが、今回はデルタ関数を基底に選びます。具体的には、$xx$軸を一定間隔$Δx\Delta x$で$x0,...,xn−1x_0,...,x_{n-1}$のように離散化します。波動関数は、$ii$番目の要素が$ψ(xi)\psi(x_i)$になっているようなベクトル$ψ⃗\vec{\psi}$で表現します。

$$ψ⃗=(ψ(x0)⋮ψ(xn−1)) \vec{\psi}=\left(\begin{matrix} \psi(x_0) \\ \vdots \\ \psi(x_{n-1}) \end{matrix}\right)$$

上記のように波動関数をベクトルで表現すると、ハミルトニアンは以下のように行列で表現できます。

$$H^=−ℏ22m(−210⋯01−21⋯001−2⋯0⋮⋮⋮⋱10001−2)+(V(x0)00⋯00V(x1)0⋯000V(x2)⋯0⋮⋮⋮⋱00000V(xn−1)) \hat{H} = -\frac{\hbar^2}{2m} \begin{pmatrix} -2 & 1 & 0 & \cdots & 0 \\ 1 & -2 & 1 & \cdots & 0 \\ 0 & 1 & -2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 1 \\ 0 & 0 & 0 & 1 & -2 \end{pmatrix} + \begin{pmatrix} V(x_0) & 0 & 0 & \cdots & 0 \\ 0 & V(x_1) & 0 & \cdots & 0 \\ 0 & 0 & V(x_2) & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 \\ 0 & 0 & 0 & 0 & V(x_{n-1}) \end{pmatrix}$$

1項目は、2階微分を差分化したものを行列で表現しています。詳細は下記の資料が分かりやすいです。
https://www.damp.tottori-u.ac.jp/~ooshida/edu/ode/g090427.pdf
このように離散化した波動関数と演算子を用いると、シュレディンガー方程式は以下のように、行列$H^\hat{H}$の固有値$ϵ\epsilon$および固有ベクトル$ψ⃗\vec{\psi}$を求める問題に帰着します。

$$H^ψ⃗=ϵψ⃗ \hat{H}\vec{\psi}=\epsilon\vec{\psi}$$

行列の固有値・固有ベクトルを求めるプログラムは、数値計算ライブラリとして提供されており、手軽に実行できます。本記事ではPython+Numpyを用いたコード例を紹介します。

Pythonコード

import numpy as np
from scipy import constants
import matplotlib.pyplot as plt

def get_second_derivative_matrix(n: int, dx: float) -> np.ndarray:
    """2階微分の行列表現を得る"""
    diag_elements = 2 * np.diag(-np.ones(n))
    non_diag_elements_upper = np.roll(np.diag(np.ones(n)), 1, axis=1)
    non_diag_elements_upper[-1, 0] = 0
    non_diag_elements_lower = np.roll(np.diag(np.ones(n)), -1, axis=1)
    non_diag_elements_lower[0, -1] = 0
    d2_mat = diag_elements + non_diag_elements_upper + non_diag_elements_lower
    # 境界条件（無限ポテンシャル井戸）により、端点を固定
    d2_mat[0, :] = 0
    d2_mat[-1, :] = 0
    return d2_mat / dx / dx

# パラメータ設定
n = 500 # グリッド数
x = np.linspace(-1, 1, n) * 100 # 長さの単位は[nm]
dx = x[1] - x[0] # 離散化幅 [nm]
x = x[:, None] # 列ベクトルに変換

l = 50 # 井戸の端L [nm]。Pythonの慣習に倣い、小文字で表現。
inf = 1e10 # 無限を大きな数で近似
potential = np.diag(inf * (np.abs(x) > l).flatten())

# 二階微分行列の生成
second_derivative_matrix = get_second_derivative_matrix(n, dx)

hbar_nano = 6.626e-16 / (2 * np.pi) # 本プログラムでは長さの単位としてnmを用いてるので、数値を調整

# ハミルトニアンの行列表現
hamiltonian = -hbar_nano ** 2 / (2 * constants.electron_mass) * second_derivative_matrix+ potential

# ハミルトニアン行列の固有値と固有ベクトルを求めることで、シュレディンガー方程式を解く
epsilon, psi = np.linalg.eigh(hamiltonian)

# n=5までの解を可視化
n = 5
fig, ax = plt.subplots(ncols=2, figsize=(6.4 * 2, 4.8))
for i in range(n):
    r = i / (n - 1)
    b = 1 - r
    ax[0].plot(x, np.real(psi[:, i]), color=[r, 0, b], label=f'n={i+1}')
    ax[1].plot(x, np.imag(psi[:, i]), color=[r, 0, b], label=f'n={i+1}')
ax[0].set_title(r'$Re(\psi)$')
ax[1].set_title(r'$Im(\psi)$')
for i in [0, 1]:
    ax[i].set_xlabel('x [nm]')
    ax[i].legend()
plt.show()

# n=5までのエネルギー固有値を出力
epsilon_ev = epsilon * (constants.nano ** 2) / constants.electron_volt
for i in range(5):
    print(f'n={i+1}: {epsilon_ev[i]:.5f} [eV]')

n=1: 0.00004 [eV]
n=2: 0.00015 [eV]
n=3: 0.00033 [eV]
n=4: 0.00059 [eV]
n=5: 0.00093 [eV]

解析解との比較

これまでPythonを用いて数値的にシュレディンガー方程式を解いてきました。一次元無限井戸型ポテンシャルの問題は解析的にも解くことができます。

$$ψn(x)={1Lcos⁡(nπ2Lx)(n∈{2,4,...},∣x∣≤L)1Lsin⁡(nπ2Lx)(n∈{1,3,...},∣x∣≤L)0(∣x∣>L) \psi_n(x)=\begin{cases} \sqrt{\frac{1}{L}}\cos\left(n\frac{\pi}{2L}x\right)& (n \in \{2, 4, ...\} , |x|\le L) \\ \sqrt{\frac{1}{L}}\sin\left(n\frac{\pi}{2L}x\right)& (n \in \{1, 3, ...\}, |x|\le L) \\ 0 & (|x|>L) \end{cases}$$

$$ϵn=n2π2ℏ28mL2 \epsilon_n = n^2\frac{\pi^2 \hbar^2}{8mL^2}$$

比較のため、n=5まで解析解を以下に示します。

n=1: 0.00004 [eV]
n=2: 0.00015 [eV]
n=3: 0.00034 [eV]
n=4: 0.00060 [eV]
n=5: 0.00094 [eV]

波動関数は定数倍を除いて特徴がよく一致しています。また、エネルギー固有値についても、数値解では若干の過小評価があるものの、概ね一致しています。エネルギー固有値が過小評価されているのは、計算の都合上井戸の外側のポテンシャルを無限ではなく有限の大きい数で代用していることに由来するのではないかと思います。(参照: 有限井戸型ポテンシャル)

まとめ

本記事では、一次元無限井戸型ポテンシャルに対するシュレディンガー方程式を行列・ベクトルの形式で表現し、Pythonを用いて数値的に解く方法を紹介しました。有限差分法を用いた二階微分の近似や、ハミルトニアン行列の構築・固有値問題の解法を通じて、解析解と一致する数値解を得ることができました。この手法は、より複雑なポテンシャルや高次元の問題にも拡張可能であり、量子力学における数値シミュレーションの基礎として有用です。