HANABI・4乗

 7月24日火曜日、晴れ。
 二度寝三度寝をやっていたら15時になっていた。
 だれ?らじ第110回を聞きながら(ののしりコーナーでにこにこした)MTG10「閃光☆HANABI団」をフライングゲットしにいった。あわよくばBloomingClover3巻もと期待したがさすがにまだのよう。
 結局日が沈むころになってラボ着。
 摂動論の復習がてら調和振動子+4乗ポテンシャルの数値計算をやる。実はこの数値計算そのものはやったことがなかった。

ハミルトニアン
{\begin{align}
H&=H_0+\lambda V\\
H_0&=\frac{1}{2}(p^2+x^2)=a^\dagger a+\frac{1}{2}\\
V&=\frac{1}{4!}x^4
\end{align}}
あまり意味がないけど丁寧に4!で割っておく。

摂動のない系の固有状態
{\begin{align}|n\rangle&=\frac{(a^\dagger)^n}{\sqrt{n!}}|0\rangle\\
H_0|n\rangle&=\left(n+\frac{1}{2}\right)|n\rangle
\end{align}}

x^4の行列要素
{\begin{align}
\langle n'|x^4|n\rangle=& \frac{1}{4}\sqrt{(n+4)(n+3)(n+2)(n+1)}\langle n'|n+4\rangle\\
 &-\left(n+\frac{3}{2}\right)\sqrt{(n+2)(n+1)}\langle n'|n+2\rangle\\
 &+\frac{3}{4}(2n^2+2n+1)\langle n'|n\rangle\\
 &-\left(n-\frac{1}{2}\right)\sqrt{n(n-1)}\langle n'|n-2\rangle\\
 &+\frac{1}{4}\sqrt{n(n-1)(n-2)(n-3)}\langle n'|n-4\rangle\\
\end{align}}

Rayleigh-Schrödinger展開
{\begin{align}
E_n(\lambda) = \epsilon_n+\lambda\langle n|V|n\rangle + \lambda^2\sum_{n\neq m}\frac{|\langle n|V|m\rangle|^2}{\epsilon_n-\epsilon_m}+O(\lambda^3)
\end{align}}
E_nが摂動後のエネルギーで\epsilon_nがもとの固有エネルギー。今の場合\epsilon_n=n+1/2

特に基底エネルギーは、
{\begin{align}
E_0(\lambda)=\frac{1}{2}+\frac{3}{4}\lambda-\frac{21}{8}\lambda^2
\end{align}}

Pythonプログラム。N=100番目までの固有状態で系を近似。LAPACKのdsyev(実対称行列を対角化する関数)で対角化している。

import numpy as np
import math
import scipy.linalg.lapack as la

N = 100

def H_4(m,n):
    if(abs(m-n)<=4):
        return (1.0/24.0)*(math.sqrt((n+4)*(n+3)*(n+2)*(n+1))/4.0 * int(m==n+4)\
                         -(n+3.0/2.0) * math.sqrt((n+2)*(n+1)) * int(m==n+2)\
                          +(3.0/4.0)*(2*n*n + 2*n + 1) * int(m==n)\
                           -(n-1.0/2.0) * math.sqrt(n*(n-1)) * int(m==n-2)\
                            +math.sqrt((n-3)*(n-2)*(n-1)*n)/4.0 * int(m==n-4))
    else:
        return 0

def E_1(n,l):
    return n+(1.0/2.0)+l*H_4(n,n)

def E_2(n,l):
    S = 0
    for k in(-4,-2,2,4):
        S += H_4(n+k,n)**2/(-k)
    return E_1(n,l) + l*l*S

for l in range(41):
    l = l/4.0
    H = np.zeros((N,N))
    for m in range(N):
        for n in range(N):
            H[m][n] += (n + 1.0/2.0) * int(m==n) + l * H_4(m,n)
    
    E,phi,info = la.dsyev(H)
    print(l,E[0],E_1(0,l),E_2(0,l))

結果
f:id:shironetsu:20180725012401p:plain
\lambda vs E_0(\lambda)のグラフ。"numerical"が近似ハミルトニアンの対角化による結果で"1st"が摂動の1次、"2nd"が2次。

 \lambda=0付近でよく合っているのでうまくいっているのだろうけれど、思っていたよりあっという間に外れてしまってなかなか惨めだ。ちょっとこれで遊んでみたい。