【第2話】自由粒子のスナップショット【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## 自由粒子の固有関数(E = 0.25[eV], 1.0[eV], 4.0[eV])
#################################################################
import numpy as np
import matplotlib.pyplot as plt
#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
#虚数単位
I = 0.0 + 1.0j

#################################################################
## 物理定数
#################################################################
#プランク定数
h = 6.62607015 * 1.0E-34
hbar = h / (2.0 * np.pi)
#電子の質量
me = 9.1093837015 * 1.0E-31
#電子ボルト
eV = 1.602176634 * 1.0E-19

#################################################################
## 物理系の設定
#################################################################
#電子のエネルギー
E1 = 1.0 * eV
E2 = 0.25 * eV
E3 = 4.0 * eV
#波数
k1 = np.sqrt(2.0 * me * E1 / hbar**2)
k2 = np.sqrt(2.0 * me * E2 / hbar**2)
k3 = np.sqrt(2.0 * me * E3 / hbar**2)
#空間刻み間隔
dx = 1E-9
#描画範囲
x_min = -2.0 * dx
x_max =  2.0 * dx
#描画区間数
NX = 500

#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#平面波の計算
phi1 = np.exp( I * k1 * x )
phi2 = np.exp( I * k2 * x )
phi3 = np.exp( I * k3 * x )

#グラフの描画(固有関数)
plt.title( u"自由粒子の固有関数(" + r"$ E = 0.25, 1.0, 4.0[{\rm eV}] $" + u")", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$ x\, [{\rm nm}]$ ", fontsize=30)
plt.ylabel(r"$ {\rm Re}\{\varphi(x)\} $", fontsize=30)

#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([x_min/dx, x_max/dx])
plt.ylim([-1.1, 1.1])

#グラフの描画
plt.plot(x/dx, phi1.real, linestyle='solid', linewidth = 5)
plt.plot(x/dx, phi2.real, linestyle='solid', linewidth = 5)
plt.plot(x/dx, phi3.real, linestyle='solid', linewidth = 5)

#余白の調整
plt.subplots_adjust(left = 0.1, right = 0.98, bottom=0.15, top = 0.95)

#グラフの表示
plt.show()

【第1話】プログラムの動作確認【Pythonコピペで量子力学完全攻略マニュアル】

#################################################################
## Pythonでグラフアニメーション
#################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
#図全体
fig = plt.figure(figsize=(10, 6))
#描画点数
NX = 100
#描画範囲
x_min = -1.0
x_max = 1.0
#座標点配列の生成
x = np.linspace(x_min, x_max, NX)
#アニメーション分割数
AN = 30
#アニメーション作成用
ims = []

#アニメーションの各グラフを生成
for a in range(AN):
    phi = 2.0 * np.pi * a / AN 
    #sin関数値を計算    
    y = np.sin( np.pi * x + phi )
    #各コマのグラフの描画
    img  = plt.plot(x, y, color="blue", linewidth=3.0, linestyle="solid" )
    ims.append( img )

#グラフタイトルの設定 
plt.title("sin function")
#x軸ラベル・y軸ラベルの設定 
plt.xlabel("x-axis")
plt.ylabel("y-axis")
#描画範囲を設定
plt.xlim([-1.0,1.0])
plt.ylim([-1.0,1.0])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)
#ani.save("output.mp4", writer=matplotlib.animation.AVConvWriter, dpi=300)
#グラフの表示
plt.show()

【量子コンピュータを作ろう!】(17)2つの独立した2重量子井戸に束縛された電子による制御NOT演算計算結果



今回は、前回示したハミルトニアン(クーロン相互作用+静電場+電磁波)を用いて計算した結果を示すよ。上左図は前々回に示したエネルギー準位の静電場強度依存性だけれども、$|10\rangle$ と $|11\rangle$ のエネルギー差に対応した電磁波を入射して生じる状態遷移を計算するよ。ちなみにこれは、第1量子ビットが $|1\rangle$ のときに第2量子ビットを反転させる制御・NOT演算に対応しているよ。あとで示すけれども、2量子ビットのラビ振動は量子ビット同士が絡み合った量子もつれ(エンタングル)状態を任意に生成することができるよ。

遷移状態:$|10\rangle$ と $|11\rangle$ の存在確率の時間依存性

次のグラフは初期状態 $|10\rangle$ に角振動数 $\omega = \Delta E /\hbar$ の電磁波を入射したときの、$|11\rangle$ との状態遷移の様子だよ。ラビ振動の結果、2つの状態を三角関数的に行ったり来たりするね。$|10\rangle$ 100%の初期状態に約 $ 26.1[{\rm ns}]$ 照射すると、$|10\rangle$ と $|11\rangle$ の存在確率は50%つづとなり、さらに約 $ 26.1[{\rm ns}]$ 照射すると、反対に $|11\rangle$ 100% の状態になるね。ちょうど存在確率が反転する時間の光は「$\pi$パルス」って呼ばれるよ。

先にも言ったけれども、このラビ振動は2量子ビット同士が絡み合った量子もつれ状態を生み出すことができるね。例えば、$|10\rangle$ と基底状態の $|01\rangle$ のエネルギー差 $\Delta E_{01}^{11}$ として、角振動数 $\omega = \Delta E_{01}^{11} /\hbar$ の電磁波を $\pi$パルスの半分を入射すると、

\begin{align}
|\Psi(t)\rangle = \frac{1}{\sqrt{2}} \left[ |10\rangle + |01\rangle \right]
\end{align}

という状態を作ることができるね。これは単に $|10\rangle$ と $|01\rangle $ の存在確率が 50%づつっていうだけでなくて、第1量子ビットを観測したときに、第2量子を観測しなくてもその状態が100%の確率で分かるという量子もつれ状態となっているよ。具体的には、第1量子ビットが $|0\rangle$ の場合には第2量子ビットは必ず $|1\rangle$ に存在し、また反対に第1量子ビットが $|1\rangle$ の場合には第2量子ビットは必ず $|0\rangle$ に存在するよ。ちなみに、この量子井戸の間隔を量子状態が変化しないようにゆっくりと、どこまで離していっても成り立つよ。そのため、通信に利用することができると考えられているよ。その場合は、量子井戸による量子もつれではなく、光子をを用いた量子もつれを利用するよ。

制御NOT演算に対応する波動関数の時間発展

次の図はラビ振動による制御NOT演算時の波動関数の時間経過を示したアニメーションだよ。左の第1量子井戸は電子が右側( $|1\rangle$ )、第2量子井戸は電子が左側( $|0\rangle$ )と右側( $|1\rangle$ )を行ったり来たりしているね。

これで量子コンピュータの2量子ビットマシンの動作原理シミュレーションは完成したよ。


【量子コンピュータを作ろう!】(16)2つの独立した2重量子井戸に束縛された電子による制御NOT演算のハミルトニアンと計算方法


いよいよ最終目的の2量子ビットによる制御NOT演算を行うために、時間発展に必要なハミルトニアンの導出と計算方法を示すよ。復習から入っていくね。量子ビットを表す2つの量子井戸(幅:$L=10[{\rm nm}]$、高さ$+\infty$)はそれぞれが壁(幅: $W=L/5=2[{\rm nm}]$、高さ $V_0=0.3[{\rm eV}]$)で仕切られた2部屋になっていて、前々回前回解説したとおり、2つの量子井戸の間隔 $R$ がちょうどよければ($R=20[{\rm nm}]$)、電子間のクーロン相互作用が存在しても、電子は左右のどちらかに存在する状態を作り出せるね。そして、それぞれの量子井戸で左側に電子が存在する状態を $|0\rangle$、右側に電子が存在する場合を $|1\rangle$ と表わして、2電子の状態は $|00\rangle$、$|01\rangle$、$|10\rangle$、$|11\rangle$ の4パターンで表すよ。

定常状態のハミルトニアンと固有関数

この系の固有関数は、2つの量子井戸でそれぞれ定義される関数

\begin{align}
\varphi_{n_1}(x_1) &\ = \sqrt{\frac{2}{L}}\sin\left[ k_{n_1} \left( x_1 + \frac{R}{2} + \frac{L}{2}\right)\right] \ , \
k_{n_1} = \frac{\pi(n_1+1)}{L} \\
\varphi_{n_2}(x_2) &\ = \sqrt{\frac{2}{L}}\sin\left[ k_{n_2} \left( x_2 – \frac{R}{2} + \frac{L}{2}\right)\right] \ , \
k_{n_2} = \frac{\pi(n_2+1)}{L} \\
\end{align}

を用いたその積

\begin{align}
\varphi_{n_1n_2}(x_1, x_2) = \varphi_{n_1}(x_1)\varphi_{n_2}(x_2)
\end{align}

で表される正規直交系で展開することができるね。ちなみに単純な積で表されるのは2つの電子は交わらないからだよ。この2電子系の正規直交関数を用いて、次のように展開することができるね。

\begin{align}
\psi_{l_1l_2}(x_1, x_2) = \sum\limits_{n_1n_2} a^{(l_1l_2)}_{n_1n_2} \varphi_{n_1n_2}(x_1, x_2)
\end{align}

$l_1$ と $l_2$ は $0$ または $1$ のどちらかを取り、$\psi_{00}(x_1, x_2)$ は $|00\rangle$、$\psi_{01}(x_1, x_2)$ は $|01\rangle$、$\psi_{10}(x_1, x_2)$ は $|10\rangle$、$\psi_{11}(x_1, x_2)$ は $|11\rangle$ にそれぞれ対応した固有関数だよ。$a^{(l_1l_2)}_{n_1n_2}$ のように上付きのインデックス $l_1, l_2$ をつけたのはそれぞれの状態で展開係数の値が異なるので明示的につけているよ。この固有関数に対応するハミルトニアンは次のとおりだよ。

\begin{align}
\hat{H}^{(0)} = -\frac{\hbar^2}{2m_e}\, \frac{\partial^2}{\partial x_1^2} + eE_xx_1 + V_1(x_1)-\frac{\hbar^2}{2m_e}\, \frac{\partial^2}{\partial x_2^2} + eE_xx_2 + V_2(x_2) + \frac{e^2}{4\pi\epsilon_0}\,\frac{1}{|x_1-x_2|}
\end{align}

$\hat{H}^{(0)}$ とインデックスに $(0)$ をつけているのは、前回すでに固有状態を計算できていて、次はこれに電磁波を加えることを想定しているからだよ。先の固有関数はこのハミルトニアンなので固有状態は

\begin{align}
\hat{H}^{(0)}\psi_{l_1l_2}(x_1, x_2) = E^{(0)}_{l_1l_2}\psi_{l_1l_2}(x_1, x_2)
\end{align}

と表されるよ。ここまでが前回の内容だね。ちなみに固有状態は

\begin{align}
\hat{H}^{(0)}|l_1l_2\rangle &\ = E^{(0)}_{l_1l_2}|l_1l_2\rangle \\
\langle l_1l_2| \hat{H}^{(0)} &\ = \langle l_1l_2|E^{(0)}_{l_1l_2}
\end{align}

とも表すことができるよ。これは後で使うね。

電磁波を入射したときのハミルトニアンと計算方法


次に状態遷移を起こすために電磁波を加えるよ。ベクトルポテンシャルを $\boldsymbol{A}(t)$ としてハミルトニアンは

\begin{align}
\hat{H}(t) = \hat{H}^{(0)} + \frac{e}{m_e} \boldsymbol{A}(t)\cdot \hat{\boldsymbol{p}}_1 + \frac{e}{m_e} \boldsymbol{A}(t)\cdot \hat{\boldsymbol{p}}_2
\end{align}

となるね。今回も1次元系で考えているので、ベクトルポテンシャルを$\boldsymbol{A}(t) = (A_x(t), 0, 0)$ として、

\begin{align}
A_{x_1}(t) &\ = A_0 \cos(kx_1-\omega t)\\
A_{x_2}(t) &\ = A_0 \cos(kx_2-\omega t)
\end{align}

なので、ハミルトニアンは

\begin{align}
\hat{H}(t) = \hat{H}^{(0)} + \frac{e}{m_e} A_0 \left[ \cos(kx_1-\omega t)\hat{p}_{x_1} +\cos(kx_2-\omega t)\hat{p}_{x_2} \right]
\end{align}

となるね。このハミルトニアンは時間に依存するので固有状態は存在しないので、波動関数 $\Psi(t,x_1,x_2)$ を先の固有関数 $\psi_{l_1l_2}(x_1, x_2)$ で展開するよ。

\begin{align}
\Psi(t,x_1,x_2) = \sum\limits_{l_1 l_2} b_{l_1l_2}(t) \psi_{l_1l_2}(x_1, x_2)
\end{align}

$b_{l_1l_2}(t)$ が展開係数で、この展開係数が時間とともに変化するよ。時間に依存するシュレーディンガー方程式

\begin{align}
i\hbar \frac{\partial }{\partial t} \Psi(t,x_1,x_2)= \hat{H}(t) \Psi(t,x_1,x_2)
\end{align}

に代入して、両辺に $\psi_{m_1m_2}^*(x_1, x_2)$ を掛けて全空間で積分すると、次のような連立微分方程式になるよ。

\begin{align}
i\hbar \frac{d b_{m_1m_2}(t) }{d t} = E^{(0)}_{m_1m_2} b_{m_1m_2}(t) + e A_0 \sum\limits_{l_1 l_2} b_{l_1l_2}(t) \langle m_1 m_2| \left[\cos(kx_1-\omega t)\frac{\hat{p}_{x_1}}{m_e} + \cos(kx_2-\omega t)\frac{\hat{p}_{x_2}}{m_e}\right]| l_1l_2 \rangle
\end{align}

$\hat{p}_{x_1}/m_e = [\hat{H}^{(0)}, x_1 ]/i\hbar$、 $\hat{p}_{x_2}/m_e = [\hat{H}^{(0)}, x_2 ]/i\hbar$ の恒等式と長波長近似($kx_1 = kx_2 \simeq 0$)を考慮すると次のようになるよ(特に長波長近似を課すことをしなくても数値計算自体は問題なくできるよ。でも表式が簡単になるね)。

\begin{align}
&\ \langle m_1 m_2| \left[\cos(kx_1-\omega t)\frac{\hat{p}_{x_1}}{m_e} + \cos(kx_2-\omega t)\frac{\hat{p}_{x_2}}{m_e}\right]| l_1l_2 \rangle\\
&\ =\frac{1}{i\hbar} \langle m_1 m_2| \left[\cos(kx_1-\omega t)\left(\hat{H}^{(0)}x_1 – x_1\hat{H}^{(0)}\right) + \cos(kx_2-\omega t)\left(\hat{H}^{(0)}x_2- x_2\hat{H}^{(0)}\right)\right]| l_1l_2 \rangle\\
&\ \simeq \frac{1}{i\hbar}\cos(\omega t)\langle m_1 m_2| \left[\hat{H}^{(0)}x_1 – x_1\hat{H}^{(0)} + \hat{H}^{(0)}x_2- x_2\hat{H}^{(0)}\right]| l_1l_2 \rangle\\
&\ = \frac{1}{i\hbar}\cos(\omega t) \left( E^{(0)}_{m_1m_2} – E^{(0)}_{l_1l_2} \right) \langle m_1 m_2|(x_1 + x_2)| l_1l_2
\rangle \\
&\ = \frac{1}{i\hbar}\cos(\omega t) \left( E^{(0)}_{m_1m_2} – E^{(0)}_{l_1l_2} \right) \int_{-\frac{R}{2}-\frac{L}{2}}^{-\frac{R}{2}+\frac{L}{2}} dx_1\int_{\frac{R}{2}-\frac{L}{2}}^{\frac{R}{2}+\frac{L}{2}}dx_2 \psi_{m_1m_2}^*(x_1, x_2)(x_1 + x_2)\psi_{l_1l_2}(x_1, x_2)\\
&\ \equiv \frac{1}{i\hbar}\cos(\omega t) \left( E^{(0)}_{m_1m_2} – E^{(0)}_{l_1l_2} \right) K^{m_1m_2}_{l_1l_2}
\end{align}

この $K^{m_1m_2}_{l_1l_2}$ は時間に依存しないので1度計算すれば良いね。この $K^{m_1m_2}_{l_1l_2}$ を用いると先の連立微分方程式は次のようになるよ。

\begin{align}
i\hbar \frac{d b_{m_1m_2}(t) }{d t} = E^{(0)}_{m_1m_2} b_{m_1m_2}(t) + \frac{e A_0}{i\hbar} \cos(\omega t) \sum\limits_{l_1 l_2} b_{l_1l_2}(t) \left( E^{(0)}_{m_1m_2} – E^{(0)}_{l_1l_2} \right)K^{m_1m_2}_{l_1l_2}
\end{align}

電磁波の角振動数が2準位間のエネルギー差 $\Delta E = E^{(0)}_{m_1m_2} – E^{(0)}_{l_1l_2}$ と表して $\omega = \Delta E / \hbar$ となるときに、2準位間を周期的に遷移するね。$\Delta E = E^{(0)}_{11} – E^{(0)}_{10}$ を与えることで、制御・NOT演算となることを次回シミュレーションするよ。


【量子コンピュータを作ろう!】(15)2つの独立した量子井戸に束縛された電子の基礎実験3:クーロン相互作用+静電場(2重量子井戸)


前回、電子同士のクーロン相互作用を考慮した独立した2つの2重量子井戸の電子の固有状態を計算したね。今回はさらに静電場を加えるね。壁の幅と高さは、$W = L/5= 2.0 \times 10^{-9} [{\rm m}] = 2.0[{\rm nm}]$、$V_0 = 3.0[{\rm eV}]$ として、量子井戸の間隔は $R = 20.0[{\rm nm}]$ とするよ。次の結果は、静電場の大きさを $E_x = 0.1\times10^{6}[{\rm V/m}]$ から $1.0\times10^{6}[{\rm V/m}]$ まで $ 0.05\times10^{6}[{\rm V/m}]$ づつ強くしたときの基底状態から第3励起状態までの電子存在確率の空間分布だよ。想定通りだけれども、ちょっとした静電場で第1励起状態と第2励起状態の縮退が解けて、2重量子井戸の壁の左右に分かれる形になったね。それぞれの量子井戸の左に電子がいる状態を $|0\rangle$、右にいる状態を$|1\rangle$ と表すと、それぞれの電子状態は固有エネルギーが低い順番に $|01\rangle, |00\rangle, |10\rangle, |11\rangle$ と表わすことができるね。



次の図は、上記の固有状態に対応するエネルギー準位の静電場の強度 $E_x$ 依存性だよ。想定通り、電場強度に比例して第1励起状態と第2励起状態の縮退が解けていく様子がわかるね。まさにシュタルク効果だね。今回は、エネルギー準位とは関係せずに電子が左か右かで $|0\rangle, |1\rangle$ を定義したけれども、エネルギー準位を基準として、2つの量子井戸の外側と内側で $|0\rangle, |1\rangle$ を定義すれば、固有エネルギーの低い順に $|00\rangle, |01\rangle, |10\rangle, |11\rangle$ と並ばせることもできるね。そして、下図の $\Delta E$ に対応する電磁波を外部から照射すれば、ラビ振動で $|10\rangle$ と $|11\rangle$ の状態遷移を起こすことができるので、制御・NOT演算が実現できるね。

次回は実際にラビ振動をシミュレーションして、制御・NOT演算を確かめてみよう!