5.1.2. 調和振動子の基底を用いた行列表現

 調和振動子 \(\hat{H}= \frac{\hat{p}^2}{2} + \frac{\hat{q}^2}{2}\) の固有状態 \(| n\rangle\) を用いて ハミルトニアンの行列表現を与える.

初頭的な量子力学の教科書 [1] で定義される生成,及び消滅演算子

\[\begin{split}&\hat{a} = \frac{1}{\sqrt{2\hbar}}(\hat{q} + i\hat{p})\\ &\hat{a}^{\dagger} = \frac{1}{\sqrt{2\hbar}}(\hat{q} - i\hat{p})\end{split}\]

を用いると,調和振動子のハミルトニアンは

\[\hat{H} = \hbar(\hat{a}^\dagger\hat{a} + \frac{1}{2})\]

として記述される.固有値問題は

\[\hat{H} | n\rangle = \hbar(n + \frac{1}{2})|n\rangle\]

である.いくつかの関係式

\[\begin{split}& \langle m | n \rangle = \delta_{m,n}\\ & \hat{a}^\dagger | n \rangle = \sqrt{n+1} | n+1\rangle\\ & \hat{a} | n \rangle = \sqrt{n} | n -1 \rangle\\ & \hat{a}^\dagger\hat{a} |n \rangle = n | n \rangle\end{split}\]

を思い出してもらって,調和振動子の基底 \(|n \rangle, |m \rangle\) で与えられる演算子等の行列表現を考える.

\[\begin{split}&\langle m | n\rangle = \delta_{m,n} \longleftrightarrow \left( \begin{smallmatrix} 1 & 0 & 0 & \cdots\\\\0 & 1 & 0 & \\\\0 & 0 & 1 & \\\\\vdots & & & \ddots \end{smallmatrix} \right)\\ &\langle m | \hat{a}^\dagger | n \rangle = \sqrt{n+1}\langle m | n+1\rangle =\sqrt{n+1}\delta_{m,n+1} \longleftrightarrow \hat{a}^\dagger = \left( \begin{smallmatrix} 0 & 0 & 0 & \cdots\\\\\sqrt{1} & 0 & 0 & \\\\0 & \sqrt{2} & 0 & \\\\\vdots & & & \ddots \end{smallmatrix} \right)\\ &\langle m | \hat{a} | n \rangle = \sqrt{n}\langle m | n-1\rangle =\sqrt{n}\delta_{m,n-1} \longleftrightarrow \hat{a} = \left( \begin{smallmatrix} 0 & \sqrt{1} & 0 & \cdots\\\\0 & 0 & \sqrt{2} & \\\\0 & 0 & 0 & \\\\\vdots & & & \ddots \end{smallmatrix} \right)\\ &\langle m | \hat{a}^\dagger\hat{a} | n \rangle = n\langle m | n\rangle = n\delta_{m,n} \longleftrightarrow \hat{a}^\dagger\hat{a} = \left( \begin{smallmatrix} 1 & 0 & 0 & \cdots\\\\0 & 2 & 0 & \\\\0 & 0 & 3 & \\\\\vdots & & & \ddots \end{smallmatrix} \right)\end{split}\]

故に,\(\hat{q}, \hat{p}\) 行列表現は

\[\begin{split}& \hat{q} = \sqrt{\frac{\hbar}{2}} (\hat{a}^\dagger + \hat{a}) \longleftrightarrow \sqrt{\frac{\hbar}{2}} \left( \begin{smallmatrix} 0 & \sqrt{1} & 0 & \cdots\\\\\sqrt{1} & 0 & \sqrt{2} & \\\\0 & \sqrt{2} & 0 & \\\\\vdots & & & \ddots \end{smallmatrix} \right),\\\end{split}\]\[\begin{split}\hat{p} = \sqrt{\frac{\hbar}{2}} (\hat{a}^\dagger - \hat{a}) \longleftrightarrow i\sqrt{\frac{\hbar}{2}} \left( \begin{smallmatrix} 0 & -\sqrt{1} & 0 & \cdots\\\\\sqrt{1} & 0 & -\sqrt{2} & \\\\0 & \sqrt{2} & 0 & \\\\\vdots & & & \ddots \end{smallmatrix} \right)\end{split}\]

理論的に行列の次元は無限次元だが,数値的には無理な相談なので \(N\) 次元までで切り捨てた \(N\times N\) 次元の \(\hat{q}, \hat{p}\) を用いてハミルトニアン

\[\hat{H} = \frac{1}{2}\hat{p}^2 + V(\hat{q})\]

を構成し固有値問題

\[\hat{H}|\Psi_n\rangle = E_n | \Psi_n\rangle\]

を求める. まず始めに自明な例であるが調和振動子 \(V(\hat{q}) = \hat{q}^2/2\) を解く. ([2] により詳しい解説とMATLABを用いたサンプルプログラムが紹介されている.次のプログラムはそのサンプルコードをpythonに翻訳したものである)

# harmonic_base1.py
import numpy as np
N = 50
hbar = 1.0
n = np.arange(1,N)
m = np.sqrt(n)
q = np.sqrt(hbar/2.0) * (np.diag(m,-1) + np.diag(m,1))
p = np.sqrt(hbar/2.0) * (np.diag(m,-1) - np.diag(m,1))*1.j
q,p = np.matrix(q), np.matrix(p)
H=0.5*p*p + 0.5*q*q
evals, evecs = np.linalg.eig(H)
print(evals[0:8].real)
$ python harmonic_base1.py
[ 0.5  1.5  2.5  3.5  4.5  5.5  6.5  7.5]

ノート

今考えているハミルトニアンの場合,打ち切られた行列の対称性が非常に良いため,行列の最後の成分(一番右下)が縮退してしまう. 即ちN=6としてハミルトニアンの行列要素を見ると

matrix([[ 0.5,  0. ,  0. ,  0. ,  0. ,  0. ],
    [ 0. ,  1.5,  0. ,  0. ,  0. ,  0. ],
    [ 0. ,  0. ,  2.5,  0. ,  0. ,  0. ],
    [ 0. ,  0. ,  0. ,  3.5,  0. ,  0. ],
    [ 0. ,  0. ,  0. ,  0. ,  4.5,  0. ],
    [ 0. ,  0. ,  0. ,  0. ,  0. ,  2.5]])

の様になっており,最後の打ち切りの影響がもろに出ている. 故に対角かする前に最後の 行列の最後の行と列は切り捨てて \(N-1\times N-1\) の行列として固有値問題(既に対角化されているが・・・)を解く方が良いであろう. 調和振動子以外の場合にこの問題が発生するかどうかは私はまだ理解していませんが, 理論的に縮退がないにも関わらず数値解が縮退する場合や \(N\to\infty\) をとっても変な場合はこの問題を疑った方が良いと思います.


 次に 固有関数 \(|\Psi_n\rangle\) について考える.

q-表示の波動関数は調和振動子の基底 \(| n\rangle\) を用いて

\[\langle q | \Psi_n \rangle = \sum_{n=0}^{N-1} = \langle q | n \rangle \langle n | \Psi_n\rangle\]
で与えられる.ここで \(\langle n | \Psi_n\rangle\) は調和振動子の基底で展開した固有ベクトルであり,
\(\langle q | n \rangle\)
\[\langle q | n \rangle = \left( \frac{1}{2^n n! \sqrt{\hbar\pi}} \right)^{1/2} e^{-q^2/2\hbar}\mathcal{H}_n(\frac{q}{\sqrt{\hbar}})\]

である.ここで,\(\mathcal{H}_n(x)\) はエルミート多項式

\[\begin{split}\mathcal{H}_0(x) &= 1\\ \mathcal{H}_1(x) &= 2x\\ \mathcal{H}_2(x) &= 4x^2 -2\\ \vdots\end{split}\]

である.エルミート多項式は漸化式

\[\mathcal{H}_{n+1}(x) = 2x\mathcal{H}_{n}(x) -2n \mathcal{H}_{n-1}(x)\]

を満たすので,原理的には有限回の四則演算で求めることができる [3] . しかし,幾つかの問題点があるあることが分る. 一つは \(n\) がある程度大きくなると非常に大きなもの同士の四則演算になり数値的に発散することである.

故に,\(n\) が比較的大きくなる場合はエルミート多項式の \(n\to\infty\) における漸近関係 が分ると良いのだが,まだ良く分ってないので既存のパッケージで([4] [5])数値的に求めることが出来る項まで使用しよう.

また,すごくいい加減に書いたプログラムだが,以下の様に実装し調和振動子の(自明だが)固有値問題を解いてみよう.

# harmonic_base.py

import numpy as np
from scipy import special

twopi=2.0*np.pi

class Harmonic(object):
    def T(self, x):
        return 0.5*x*x
    def V(self, x):
        return 0.5*x*x

class Operator(object):
    def __init__(self,dim, domain,Hamiltonian):
        self.dim = dim
        self.domain = domain
        self.Hamiltonian = Hamiltonian
        self.__set_domain(domain)
        self.opQ()
        self.opP()
        self.matrix()
        
    def __set_domain(self, domain):
        qmin,qmax = domain[0][0],domain[0][1]
        pmin,pmax = domain[1][0],domain[1][1]
        if (qmin > qmax) or (pmin > pmax):
            raise ValueError("Domain Error p in [%f,%f], p in [%f %f]\n" % (qmin,qmax, pmin,pmax))
        if pmax != -pmin:
            raise ValueError("p in (min,max) and max=-min\n")
        q = np.linspace(qmin, qmax,self.dim,endpoint=False)
        p = np.linspace(pmin, pmax,self.dim,endpoint=False)
        self.x = [q,p]
        self.h = (qmax - qmin)*(pmax - pmin)/self.dim

    def opQ(self):
        hbar = self.h/twopi
        m = np.sqrt(np.arange(1,self.dim))
        q = np.sqrt(hbar/2.0)*(np.diag(m,-1) + np.diag(m,1))
        self.q = np.matrix(q)

    def opP(self):
        hbar = self.h/twopi
        m = np.sqrt(np.arange(1,self.dim))
        p = np.sqrt(hbar/2.0)*(np.diag(m,-1) - np.diag(m,1))*1.j
        self.p = np.matrix(p)
        
    def matrix(self):
        self.matrix = self.Hamiltonian.T(self.p) + self.Hamiltonian.V(self.q)
    
    def eigen(self,sort=False, vector=False):
        evals, evecs = np.linalg.eig(self.matrix)
        evals = np.array(evals.tolist())
        coefs = [np.array(vec) for vec in evecs.transpose().tolist()]
        if sort:
            index = [i[0] for i in sorted(enumerate(evals.real), key=lambda x:x[1])]
        else:
            index = range(0,self.dim)
        if vector:
            vecs = [self.eigenvector(coefs[i]) for i in index]
            return (evals[index], vecs)
        else:
            return evals[index]
            
    def harmonic_eigenvector(self,n):
        x = self.x[0]/np.sqrt(self.h/twopi)
        coef = np.zeros(self.dim)
        coef[n] = 1.0
        Herm = np.polynomial.hermite.Hermite(coef,domain=self.domain[0],window=self.domain[0])
        return Herm*np.exp(-x*x/2.0) + 0.j
        
    def harmonic_eigenvectors(self,truncate):
        x = self.x[0]/np.sqrt(self.h/twopi)
        if self.dim < truncate:
            eye = np.eye(self.dim)
        else:
            eye = np.eye(truncate)
        Her_poly=[np.polynomial.hermite.Hermite(e,domain=self.domain[0],window=self.domain[0]) for e in eye]
        vecs = [Hn(x)*np.exp(-x*x/2.0) + 0.j for Hn in Her_poly]
        return [self.normalize(x) for x in vecs]
    
    def eigenvector(self,coef, truncate=100):
        base = self.harmonic_eigenvectors(truncate)
        max_index = self.dim if self.dim < truncate else truncate
        vec = np.array([coef[n]*base[n] for n in range(max_index)]).sum(axis=0)
        return self.normalize(vec)
        
    def normalize(self,x):
        norm = np.abs(np.sum(x*np.conj(x)))
        return x/np.sqrt(norm)

    def harmonic_eigenvector_asymptotic(self,nmin=50):
        pass

とでも書いておいて,

# harmonic.py
from harmonic_base import *
import numpy as np
import matplotlib.pyplot as plt
twopi = 2.0*np.pi
dim=100
qmin,qmax = -5.0,5.0
pmin,pmax = -5.0,5.0
domain = [[qmin,qmax],[pmin,pmax]]
op = Operator(dim,domain, Harmonic())
hbar = op.h/twopi
evals,evecs = op.eigen(sort=False, vector=True)
#evals = op.eigen(sort=False, vector=False)

fig=plt.figure(figsize=(6,6))

for n in range(dim):
    if n > 50: break
    plt.plot(op.x[0],np.abs(evecs[n]*np.conj(evecs[n]))+evals[n].real, '-',linewidth=3)

y=Harmonic().V(op.x[0])
plt.plot(op.x[0],y,'-k',linewidth=2)
plt.ylabel(r"$V(q)$,$|\psi_n(q)|^2$",fontsize=20)
plt.xlabel(r"$q$",fontsize=20)
plt.ylim(0.0,2.0)
plt.xlim(-3.0,3.0)
plt.savefig("eigen_state_HB.png")
#plt.show()
plt.close()
fig=plt.figure(figsize=(6,6))

x = np.arange(0,dim-1)
plt.plot(x,evals[0:dim-1].real/hbar,'-o')
plt.xlabel(r"$n$ (quantum number)")
plt.ylabel(r"$E_n/\hbar$")
plt.legend(loc=2,numpoints=1)
plt.savefig("eigen_valus_HB.png")
#plt.show()

を実行すると,次の結果が得られる.

fig1 fig2

Footnotes

[1]例えば J. J. Sakurai (櫻井純), 現代の量子力学 (1989) 吉岡書店
[2]H. J. Korsch, H. J. and M. Glück, Euro. J. Phys, 23 (2002) 413
[3]エルミート多項式の直交性を確かめる
[4]Numpy:A Hermite series class.
[5]Scipy:Special functions