「14日で作る量子コンピュータ」を読んでみる。4日目 量子井戸の再現

はじめに

　今回は量子井戸内部の波動関数と量子井戸に電場を加える方法についてまとめる。

4 量子井戸

4.1 井戸型ポテンシャル内部の波動関数

　ここでは井戸型ポテンシャルと言われる、以下のようなポテンシャル領域を考える。

V(x)=\left\{\begin{array}{ll}
0 & (x < |\frac{L}{2}|)\\
\infty & (x > |\frac{L}{2}|)
\end{array}
\right.
\quad (4.1)

$\frac{L}{2}$より遠い領域ではポテンシャルが無限なので、量子がポテンシャルの壁を越えられないため、この領域では量子の存在確率はゼロ、つまり、

$$
\varphi(x)=0(x > \frac{L}{2})
$$

　したがって境界条件は、

$$
\varphi(\frac{L}{2})=0, \varphi(-\frac{L}{2})=0
\quad (4.2)
$$

　である。これを整理すると、

$$
(4.3: 量子井戸中の時間依存因子)
$$

\varphi_n(x)=\left\{\begin{array}{ll}
Asin(k_nx) & (n=1,3,5,7,...)\\
Acos(k_nx) & (n=0,2,4,8,...)
\end{array}
\right.,
k_n=\frac{\pi(n+1)}{L}

　条件より量子は$x<|\frac{L}{2}|$の範囲に存在しているので、規格化条件は

$$
1=\int_{-\frac{L}{2}}^{\frac{L}{2}}|\varphi(x)|^2dx
$$

　これに(4.3)を代入すると$A=\sqrt{\frac{2}{L}}$とわかるので、これに時間依存因子をかけることで波動関数を求めることができる。

$$
\psi(x,t)=\sqrt{\frac{2}{L}}sin\left[ k_n(x+\frac{L}{2}) \right]e^{-i\omega_n t}, k_n=\frac{\pi(n+1)}{L}
$$

また

$$
E_0 = \frac{\hbar^2k_0^2}{2m}=0.376[eV]
$$

$$
T_0=\frac{2\pi}{\omega}=\frac{2\pi \hbar}{E_0}
$$

knはnによって決まるので、

$$
E_n = E_0 \times (n-1)^2, T_n=\frac{T_0}{(n-1)^2}
$$

4.2 量子井戸に電場を加える

　量子井戸内のポテンシャルがゼロである際は上のように解析的に波動関数を求めることができるが、解析的に求めることが難しい場合もある。

$$
V=eEx
$$

　したがって、シュレディンガー方程式に与える。ポテンシャルエネルギーは

$$
(4.4: 電場を加えた際のポテンシャルエネルギー)
$$

V(x)=\left\{\begin{array}{ll}
eE_xx & (|x|\leq \frac{L}{2})\\
+\infty & (|x| > \frac{L}{2})
\end{array}
\right.

　この場合のシュレディンガー方程式は先ほど(ポテンシャルが存在しない場合)のように解析的に解けることは少ないのが、数値的に計算する方法がある。
　シュレディンガー方程式を数値的に解く常套手段は、「計算したい未知の解を基地の解(正規直交系)の和で展開して、展開係数を数値的に求める」と言う流れである。フーリエ展開と同じような感じらしい。
1.まず解きたいハミルトニアンを$\hat{H}$、固有関数を$\varphi(x)$、エネルギー固有値を$E$とすると、固有方程式は

$$
\hat{H}\varphi(x)=E\varphi(x)
\quad (4.5)
$$

2.未知状況のハミルトニアンを、既知のハミルトニアン$\hat{H}_0$とそれ以外の項$V$の和として表す。

$$
\hat{H}=\hat{H}_0+V
\quad (4.6)
$$

• この時$\hat{H}_0$についての固有関数と固有エネルギーを$\varphi_n(x),E_n$とおくと次の式を満たす

\hat{H}_0\varphi_n(x)=E_n\varphi_n(x)
\quad (4.7)

3.未知の固有関数$\varphi(x)$が既知の固有関数$\varphi_n(x)$の和で次のように展開できると考える。

$$
\varphi(x)=\sum_{n}a_{n}\varphi_n(x)
\quad (4.8)
$$

4.上の式の展開係数($a_n$)を求める。(4.8)の両辺に$\varphi_m^*(x)$をかけて全空間で$x$について積分すると。

\int_{-\infty}^{\infty}\varphi_m^*(x)\varphi(x)dx=
\sum_{n} a_n \int_{-\infty}^{\infty}\varphi_m^*(x)\varphi_n(x)dx=
\sum_{n}a_{n}\delta_{mn}=am

  まとめると

a_n=\int_{-\infty}^{\infty}\varphi_{n}^{*}(x)\varphi(x)dx
\quad (4.9)

5.ここで(4.5)に(4.6)と(4.8)を代入した式を出しておく。

\sum_{n}a_{n}(\hat{H}_0+V)\varphi_n(x)=E\sum_{n}a_n\varphi_n(x)
\quad (4.10)

6.ハミルトニアン$\hat{H}_0$は(4.7)より既知の固有関数が存在しているので次のように変換できる。

$$
\sum_{n}a_n(E_n+V)\varphi(x)=E\sum_n a_n \varphi_n(x)
\quad (4.11)
$$

7.(4.11)から$a_n$を抜き出すために、両辺に$\varphi_m^*(x)$をかけて全空間で積分する。

\sum_n a_n \int_{-\infty}^{\infty}\varphi_m^*(x)(E_n+V)\varphi_n(x)dx=E\sum_n a_n \int_{-\infty}^{\infty}\varphi_m^*(x)\varphi_n(x)
\quad (4.12)

  右辺の積分内部は正規直交関係式であり$\varphi_n$は正規直交であり、$E$は定数なので

$$
(4.12の右辺)=E\sum_n a_n \delta_{nm}=Ea_m
$$

  左辺第１項目の$E_n$は空間に依存しない定数なので、これも同じように正規直交交換式が利用でき

$$
(4.12の左辺第１項目)=\sum_n a_n E_n \int_{-\infty}^{\infty}\varphi_m^{*}(x)\varphi_n(x)dx=\sum_n a_n E_n \delta_{nm}=a_m E_m
$$

  左辺の第２項目の$V$は空間に依存するので、この積分はそのまま置いておいて

$$
(4.12の左辺第２項目)=\sum_n a_n \int_{-\infty}^{\infty}\varphi_m^{*}(x)V\varphi_n(x)=\left< m|V|n \right>
$$

  これらをまとめて(6.9)は以下のようになる

$$
a_m E_m + \sum_n a_n =Ea_m
\quad (6.10)
$$

8.(6.10)の左辺について以下のようにおく　

M\equiv
\left(\begin{matrix}
E_0+<0|V|0> & <0|V|1> & <0|V|2> & \ldots\\
<1|V|0> & E_1+<1|V|1> & <1|V|2> & \ldots\\
<2|V|0> & <2|V|1> & E_2+<2|V|2> & \ldots\\
\vdots & \vdots & \vdots & \ddots\\
\end{matrix} \right)

$$
\boldsymbol{a}=\left(
\begin{matrix}
a_0\
a_1\
a_2\
\vdots
\end{matrix}
\right)
$$

  これらを用いると(6.10)は

$$
M\boldsymbol{a}=E\boldsymbol{a}
\quad (6.11)
$$

　この式のうち$\boldsymbol{a},E$は未知$M$の内部は数値計算によって近似できるので、行列による固有方程式として解く「固有値問題」であると言うことがわかる。次はこれを実装してみる。

4.3 量子井戸の実装

　書籍の方では、これらの計算を関数としてそれぞれ行っていたが、量子井戸の動作に必要なものだけをまとめたQWellクラスを実装してみようと思う。

4.3.1 動作内容

　量子井戸を作るに当たって、以下のクラスメソッドを実装する

• phi(n, x): ポテンシャルが無い場合の空間依存因子を出力
• V(x): ポテンシャル項を出力
• Energy(n): 第n励起状態のエネルギー
• integral_eigenMatrixElement(x, n1, n2):$\boldsymbol{M}$の要素の数値計算を行う際に積分対象となる関数
• eigenMatrix(): $\boldsymbol{M}$を求める関数
• solve_eigenValue(): 固有エネルギー、固有ベクトルを求める関数
• get_proximatePhi(n): 固有ベクトルを使って、第n励起状態の波動関数を返す関数
• set_electricField(Ex): 電場の強さを再設定しその際の固有エネルギー、固有ベクトルを求める関数。計算する際にはそれより前にこの関数を実行する

4.3.2 実装

　最初にクラスの外形を作っていく

quantumWell

class QWell:
    def __init__(self):
        pass

    def phi(self, n, x):
        pass

    def V(self, x):
        pass

    def Energy(self, n):
        pass

    def integral_eigenMatrixElement(self, x, n1, n2):
        pass

    def eigenMatrix(self):
        pass

    def solve_eigenValue(self):
        pass

    def get_proximate_phi(self, n):
        pass

    def set_electrixField(self, Ex):
        pass

　
　次にクラス変数を設定していく。まず初期化

quantumWell

def __init__(self):
    # 量子井戸の幅
    self.L = 1.0e-9
    # 計算区間
    self.x_min = -self.L / 2.0
    self.x_max = self.L / 2.0
    # 状態数
    self.n_max = 10
    self.dims = self.n_max + 1
    # 電場の強さ
    self.Ex = 1.0e9
    # 電場に応じて
    self.set_electricField(self.Ex)

次に必要な定数を求める関数を実装していく

quantumWell

#　ポテンシャルが存在しない場合の空間依存因子
def phi(self, n, x):
    kn = math.pi * (n + 1)/ self.L
    phi_value = math.sqrt(2.0/self.L) * math.sin(kn * (x + self.L/2.0))
    # return phi_value * math.sqrt(self.L/ 2.0) # Scaling
    return phi_value 
#　ポテンシャル項
def V(self, x):
    return (e * self.Ex * x)
#　(第n励起状態のエネルギー)
def Energy(self, n):
    kn = math.pi * (n + 1)/self.L
    En = (hbar ** 2) * (kn ** 2) / (2 * me)
    return En

　次はこのクラスの核となる行列計算の関数を実装する

quantumWell

# 固有対を求めるための行列の要素を求める
def integral_eigenMatrixElement(self, x, n1, n2):
    return self.phi(n1, x) * self.V(x) * self.phi(n2, x) / eV

# 固有値問題の対象となる行列
def eigenMatrix(self):
    matrix = []
    for n1 in range(self.dims):
        col = []
        for n2 in range(self.dims):
            result = integrate.quad(
                self.integral_eigenMatrixElement,
                self.x_min,
                self.x_max,
                args=(n1, n2)
            )
            real = result[0]
            imag = 0.0j

            En = self.Energy(n1)/eV if (n1 == n2) else 0
            col.append(En + real)
        
        matrix.append(col)
    
    return np.array(matrix)

def solve_eigenValue(self):
    matrix = self.eigenMatrix()
    # print("m",matrix)
    result = LA.eig(matrix)
    # 固有エネルギー
    eig = result[0]
    # 固有ベクトル
    vec = result[1]

    index = np.argsort(eig)
    eigenvalues = eig[index]

    vec = vec.T
    vectors = vec[index]

    return eigenvalues, vectors

def get_proximatePhi(self, n):
    # 固有エネルギーの数に対応する
    vector = self.eigenvectors[n]
    m = list(range(len(vector)))
    return lambda x:sum(map(lambda m,an,x: an * self.phi(m, x), 
        m, vector, x * np.ones_like(vector)))

# 電場の強さを設定する
def set_electricField(self, Ex):
    self.Ex = Ex
    self.eigenvalues, self.eigenvectors = self.solve_eigenValue()

　クラスの実装は以上である。次はこれを用いて量子井戸内部の様子を観察してみる。

4.3.3 量子井戸内部の様子の観察

　次は量子井戸内部の観察を行っていく。上のクラスを用いることで簡単に実装できる。

StarkEffect.py

import matplotlib.pyplot as plt
import numpy as np
from quantumWell import QWell

# 今回は基底状態に対して電場を加える。
n_state = 0
# ----- x軸の設定 -----
# 実験を行うxの最小値と最大値
x_min = - 0.5e-9
x_max = -x_min
# 分布をプロットするx座標
xl_real = np.linspace(x_min, x_max, 500)
# プロットする際のx軸の座標(単位をnmに変換している)
xl_show = xl_real * 1.0e9

# ----- 量子井戸の定義 -----
qwell = QWell()

# ----- 電場の設定 -----
# ここでは1.0e9ずつ0から1.0e10まで増やしていく
de = 1.0e9
E_min = 0 * de
E_max = 10 * de
NE = 10
el = np.linspace(E_min, E_max, NE + 1)

# ----- 分布の描画 -----
# グラフの設定
fig = plt.figure(figsize=(10, 6))
plt.title(f"Existence probability at position (n={dim})")
plt.xlabel("Position[nm]")
plt.ylabel("|phi|^2")
plt.xlim([-0.5, 0.5])

# 数値計算
for e in el:
    qwell.set_electricField(e)
    # 近似した波動関数を取得する
    phi_fn = qwell.get_proximatePhi(dim)
    phi_fn_abs = lambda x: (abs(phi_fn(x)) ** 2 ) 
    phi_value_col = list(map(phi_fn_abs, xl_real))

plot.show()

4.3.4 実行結果

これを実行することで以下のようなグラフが得られる

電場が強くなっていくに連れて左側に分布が偏っていくのが確認できた。