5.2. 量子写像系

 より発展した量子力学を考える際に Classical Dynamics (離散力学系) でも用いたkicked rotator は非常に強力な模型であるので, 量子化されたkicked rotort を考えよう [1] [2]

すなわちハミルトニアンは

\[H(\hat{q},\hat{p},t) = T(\hat{p}) + V(\hat{q})\sum_n \delta(t-n), \quad\text{where}\quad [\hat{q},\hat{p}] = i\hbar\]

で定義される. ここで,\([\hat{q},\hat{p}] = \hat{q}\hat{p} - \hat{p}\hat{q}\) は正準交換関係である. kickから次のkickまでの単位時間発展演算子(Unitary Operator)は

\[\hat{U} = \exp[-\frac{i}{\hbar}T(\hat{p})]\exp[-\frac{i}{\hbar}V(\hat{q})] \quad \text{or,} \quad \hat{U} = \exp[-\frac{i}{\hbar}V(\hat{q})]\exp[-\frac{i}{\hbar}T(\hat{p})]\]

として与えられる. すなわち初期条件 \(|\psi_0\rangle\) の時間発展は \(\hat{U}^{m}|\psi_0\rangle\ (m=1,2,3\cdots)\) できまる.

 初期条件を \(\langle p_0 | \psi_0\rangle = \delta(p_0)\) として状態の時間発展 \(|\psi_1\rangle = \hat{U} |\psi_0\rangle\) を求めてみよう.

\[\begin{split}\langle p_1 | \psi_1 \rangle &= \int dp_0 \langle p_1 | \hat{U} | p_0\rangle\langle p_0 | \psi_0\rangle\end{split}\]

であり,Feyman核の

\[\begin{split}\langle p_1|\hat{U} | p_0\rangle & = \int dq_1 \langle p_1 | e^{-\frac{i}{\hbar}V(\hat{q})}| q_1 \rangle \langle q_1 | e^{-\frac{i}{\hbar}T(\hat{p})}| p_0 \rangle\\ & = \frac{1}{2\pi\hbar} \int dq_1 \exp [-\frac{i}{\hbar}\{T(p_0)+V(q_1) + (p_1 - p_0)q_1\}]\end{split}\]

の積分をフーリエ変換を用いて計算すればよいので,まず実直な方法を試す.

# -*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

twopi=2.0*np.pi

def T(x,tau):
    return x**2/2.0*tau

def V(x,k,tau):
    return k*np.cos(x)*tau
	
def qmap(ivec ,q,p,Tfunc ,Vfunc ,hbar):
    dim = len(q)
    fvec = np.zeros(dim, dtype=np.complex128)
    for i in range(dim): # loop of p’
        weight = 0.0+0.j
        for j in range(dim): # loop of q
            exp_part = Vfunc(q[j]) + Tfunc(p) + (p[i]-p)*q[j]
            weight += np.sum(np.exp(-1.j*exp_part/hbar)*ivec)
            fvec[i] = weight
    return fvec.real/dim + 1.j*fvec.imag/dim

# setting of scale information
dim = 70 # Hilbert Space dim 
qmin = 0.0
qmax = twopi
pmin = -twopi*2
pmax = twopi*2
q = np.linspace(qmin, qmax, dim, endpoint=False)
p = np.linspace(pmin, pmax, dim, endpoint=False)
h = (qmax - qmin)*(pmax - pmin)/dim # effective planck constant
hbar = h/twopi

# map parameter
tau =1.0
k=1.0
Tfunc = lambda x: T(x,tau)
Vfunc = lambda x: V(x,k,tau)

# initial state

ivec = np.array([1e-18+0.j for i in range(dim)])
ivec[dim/2] = 1.0
ivec0 = np.copy(ivec)
iter = 7
for i in range(iter):
    print("step:", i)
    abs2 = np.abs(ivec*np.conj(ivec))
    par = 1.*i/iter
    plt.plot(p,abs2,'-',label='%d step' % i,color=plt.cm.jet(par),linewidth=3)
    
    ivec = qmap(ivec ,q,p,Tfunc ,Vfunc ,hbar) # iteration

plt.ylabel(r"$|\psi_n(p)|^2$",fontsize=28)
plt.xlabel(r"$p$",fontsize=28)
plt.ylim(1e-28,1)	
plt.semilogy()
plt.legend()
plt.tight_layout()
plt.savefig("evolv_prep.png")
plt.show()
_images/evolv_prep.png

フーリエ変換はFFTアルゴリズムを用いるのが実用上(精度的,速度的)重要なので、FFTを用いて書き直す. 書き直すといっても

def qmap(ivec, q,p, Tfunc, Vfunc, hbar):
    op0 = np.exp(-1.j/hbar*Tfunc(p))		
    op1 = np.exp(-1.j/hbar*Vfunc(q))
    vec1 = op0*ivec
    vec2 = op1*np.fft.fft(vec1)
    return np.fft.ifft(vec2)

だけ,変更すれば良い.結果は期待した通り実直なフーリエ変換と同一の結果をもたらす.

ノート

本例題ではこの問題は無いが,FFTはアルゴリズム上返還後の配列の格納の順番がずれているので注意.

_images/evolv_prep_fft.png
[1]“M. V. Berry, N.L. Balazs, M. Tabor and A. Voros”, Ann. Phys., 122 (1979) 26
[2]“G. Casati, B. V. Chirikov, F. M. Izraelev and J. Ford”, Lecture Note in Physics, 93 (1979) 334