調和振動子基底展開による時間依存 Schrödinger 方程式の数値解法

2013 年 2 月
福井大学 工学部 物理工学科
平成 21 年度入学
09380337
台 亮嗣

目次
第 1 章 序論

2

第 2 章 量子力学の基礎
2.1 Schrödinger 方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 時間に依存する Schrödinger 方程式 . . . . . . . . . . . . . . . .
2.1.2 時間に依存しない Schrödinger 方程式 . . . . . . . . . . . . . . .
2.2 基底による展開 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 正規直交基底 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 完全性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 ユニタリー変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 虚時間発展法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
3
3
4
4
4
5
5
6

第 3 章 調和振動子基底3.1 一次元調和振動子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 調和振動子の行列表示 . . . . . . . . . . . . . . . . . . . . . . . . . . .

8
8
10

第4章
4.1
4.2
4.3
4.4
4.5

数値計算の精度
基底の切断（ truncation ） . . . . . . . . . . . . . . . . . . . . . . . . .
誤差の微小時間に対する依存性 . . . . . . . . . . . . . . . . . . . . . .
エネルギー期待値の保存 . . . . . . . . . . . . . . . . . . . . . . . . . .
ノルムの保存 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
冷却（ 虚時間発展）により得た固有状態の精度 . . . . . . . . . . . . . .

14
14
14
16
17
19

第 5 章 波束の時間発展
5.1 非調和振動子 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 トンネル効果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21
21
25

第 6 章 結論

29

参考文献

30

謝辞

31

付録 プログラムソースコード

32

1

第 1 章 序論
19 世紀まで，ニュートンにより始められた力学と，マクスウェルが確立した電磁気学をはじめとする，古典物理学により宇宙の森羅万象はすべて説明できると考えられていた．これを根本から覆したのが，19 世紀の最後の年に登場した量子力学である．量子力学は，非相対論的な範囲に話を限れば，理論体系は完成されており，電子工学や化学その他の様々な分野に応用され，今日ではもはや物理学の新理論ではなくなっている．
我々の日常生活からは，想像し難い状態ベクトルと呼ばれる抽象的な量を用いなけ
れば，物質を構成する原子などの科学の根幹ともいえる極微な世界は記述できないというところに，ある種の哲学的なおもしろさも，量子力学にはあるように思う．
量子力学はその解釈の難しさはあるが，ニュートン力学で登場するニュートンの運動方程式と同様に，基本的には微分方程式を扱うことで，系の状態を理解することができる．1 階の微分方程式である Schrödinger 方程式を数値的に解き，系の状態の時間発展を求めていく過程で，量子力学の基礎を学び直すことを本研究のひとつの目的としている．
本研究では，一次元ポテンシャル内の一粒子の量子力学的状態ベクトルの時間発展
を，時間に依存する Schrödinger 方程式を解くことで求める．その手法として，まず状
態ベクトルを調和振動子基底で展開する．基底で状態を展開することで時間に関わる部分は展開係数のみになるため，この展開系数の時間発展を Runge-Kutta 法で数値的に解けば状態ベクトルの時間発展も求まる．一次元量子系は実空間表示や平面波展開などの手法で扱われることが多いのだが，本研究では巨大な次元の調和振動子基底で一次元量子系を扱うことができるかど うかを論ずる．

2

第 2 章 量子力学の基礎
この章では量子力学の基本方程式である Schrödinger 方程式について述べたあと，その解を記述するためのいくつかの概念を説明する．

2.1

Schrödinger 方程式

2.1.1

時間に依存する Schrödinger 方程式

ポテンシャル V (x, t) で表される力を受けて運動する質量 m の粒子の量子力学的状態を表す状態ベクトル | ψ(x, t) i は，時間に依存する Schrödinger 方程式ih̄∂| ψ(x, t) i = Ĥ| ψ(x, t) i
∂t

(2.1)

の解として与えられる．ここで h̄ はプランク定数の 2π 分の 1 倍であり，また Ĥ は系のエネルギーを表す自己共役演算子で，ハミルトニアンと呼ばれ，

Ĥ = −

h̄2 ∂ 2
+ V (x, t)
2m ∂x2(2.2)と表される．
Schrödinger 方程式は，時間に関して 1 階の微分方程式だから，初期状態（時刻 t = 0 における状態ベクトル | ψ(x, 0) i ）を与えれば，解（ t 6= 0 における状態ベクトル | ψ(x, t) i ）が一意的に定まる．なお古典論でも量子論と同様に，系が閉じている限り（ 外界から予測不能な影響を受けることがなければ ），時刻 t = 0 における状態を与えれば，t 6= 0における状態が一意的に定まるが，古典力学の運動方程式は時間に関して 2 階の微分方程式なので，状態は力学変数の時間微分をも指定しないと定まらない．

2.1.2

時間に依存しない Schrödinger 方程式

ハミルトニアン Ĥ が時間 t を含まない，即ちポテンシャルエネルギーが V (x, t) = V (x)と表される場合には，定常状態と呼ばれる特別な状態が存在する．定常状態とは状態
ベクトルが| ψ(x, t) i = e−iωt | ϕ(x) i

3

(2.3)

のように表記できる状態であり，時間には位相因子 e−iωt を通じてのみ依存する．状態ベクトルの空間部分である | ϕ(x) i は次式のような時間に依存しない Schrödinger 方程
式に従う．
(

)

h̄2 ∂ 2−
+ V (x) | ϕ(x) i = ε| ϕ(x) i.
2m ∂x2(2.4)

ここで ε はハミルトニアン Ĥ の固有値でエネルギー固有値と呼ばれる．

2.2

基底による展開

状態ベクトルを基底で展開するとはど ういうことであるかを，以下では順を追って
説明する．

2.2.1

正規直交基底複素ベクトル空間（ C 上の線形空間）に内積が定義されている場合は，元々の基底の適当な線形結合をとることで，互いに直交するように基底を選び直すことがシュミットの直交化法などによりいつでも可能である．基底を構成する状態ベクトルを自然数
n = 1, 2, 3 · · · で番号付けして | n i と表し，| n i と | n0 i の内積を h n | n0 i と表すとき，基底を次の条件を満足するように選ぶと便利である．
{
0hn|n i =1 (n = n0 )0 (n =
6 n0 )}= δn,n0 .

(2.5)

この条件を満たす状態ベクトルの集合 {| n i} を正規直交基底（ 正規直交系）と呼ぶ．

2.2.2

完全性

複素数体 C 上の完備な内積空間 H を複素ヒルベルト空間と呼ぶ．内積空間 H が完備であるとは，H 内のどんなベクトル列 | ψ1 i, | ψ2 i, · · · でも，もしそれがlim k| ψn i − | ψm ik = 0

n,m→∞

(2.6)

を満たすベクトル列であれば，必ずその収束先 | ψ i，つまりlim k| ψn i − | ψ ik = 0

n→∞

(2.7)

を満たすベクトル | ψ i が H 内に存在することである．ここで k · · · k は状態ベクトルのノルムを表す．
どんなベクトル | ψ i ∈ H も，同じ H 上のベクトルの集合 {| n i} の線形結合として表せるとき，{| n i} は H の完全系を成すという．ひとつの自己共役演算子の固有ベクト
ルをもれなく集めてくれば，それは完全系を成すことが知られている．従って，自己4共役演算子 Â の固有ベクトルを | n i とすると，任意のベクトル | ψ i ∈ H は，n が離散
的な場合は

|ψi =

∞
∑

ψn | n i

(2.8)

n=1

のように，適当な係数 ψn ∈ C を用いて展開できる．
なお，「 完備」と「完全」は共に ”complete”の訳語であり，「 完全系」でなく「完備系」と訳す流儀もあることを申し添えておく．

2.3

ユニタリー変換

ベクトル空間からベクトル空間への 1 体 1 対応の線形写像で，ノルムを変えないものをユニタリー変換と呼ぶ．
時間に依存する Schrödinger 方程式 (2.1) 式の解は，次式のように時刻 t = 0 における状態ベクトル | ψ(0) i を時間発展演算子と呼ばれるユニタリー演算子でユニタリー変換した形をとる．

| ψ(t) i = Û (t)| ψ(0) i.

(2.9)

状態ベクトルは , 時間とともに連続的にユニタリー変換される（ 観測過程を除く）．
(2.2) 式のハミルトニアンで V (x, t) = V (x)(ポテンシャルが時間に依存しない場合) と書けるときは，(ĤÛ (t) = exp
ih̄)(2.10)と形式的に書き表すことができる．

2.4

虚時間発展法実時間 t に虚数単位 i をかけたものを虚時間と呼ぶ．
(2.1) 式の Schrödinger 方程式を考える．

ih̄∂| ψ(x, t) i = Ĥ| ψ(x, t) i.
∂t(2.11)(2.8) 式より，任意の状態ベクトルは固有ベクトルで展開できるため，初期状態を| ψ(x, 0) i =∞∑ψn | n i,(2.12)

n=1

Ĥ| n i = εn | n i (n = 1, 2 · · ·)ε1 ≤ ε2 ≤ ε3 ≤ · · · ,5(2.13)

と展開すると，(2.11) 式の解は
(

∞
∑

)

εn
ψn exp −i t | n i| ψ(x, t) i =h̄
n=1(2.14)と記述される．ここで，時間 t の虚数単位倍を τ で表し（ 即ち τ = it, t, τ は実数），(2.14) 式を τ（ 虚時間と呼ばれる）で書き換えると
(

∞
∑

)

εn
| ψ(x, τ ) i =ψn exp − τ | n ih̄
n=1(2.15)

を得る．したがって充分大きな虚時間 τ に対しては，状態ベクトルの基底 | 1 i の展開係数に対して基底 | n i (n > 0) の展開係数は
(

)

(

) =

ψn exp − εh̄n τ
ψ1 exp − εh̄1 τ

{

}

ψn
τ
exp − (εn − ε1 ) −→ 0
ψ1h̄

(2.16)(τ → ∞)のように高位の無限小となるので τ → ∞ の極限で

| ψ(x, τ ) i ∝ | n i(2.17)となる．このことは虚時間発展とともに状態ベクトルが最低エネルギーの状態に近づいていくことを意味する．本研究では，これを利用して初期状態の波動関数を準備する．

2.5

Runge-Kutta 法

本研究では，Runge-Kutta 法で時間発展演算子 Û (t) の作用を数値的に扱う．ここでは Runge-Kutta 法について説明する．
関数 y = f (x) の一階微分 y 0 が y 0 = F (x, y) の形の常微分方程式で与えられ，y0 =f (x0 ) の初期条件のもとに y = f (x) を求める初期値問題を考える．
最も簡単な方法は，y0 = f (x0 ) から始めて，近似式

f (x + h) ≈ f (x) + f 0 (x)h(2.18)

を漸化式として f (x0 + h), f (x0 + 2h), · · · を順に求めていく Euler 法（ 誤差が h の 1 乗に比例する．
）である．
計算をより高精度なものにするには，(2.18) 式のような簡単な近似式の代わりに

f (x + h) ≈ f (x) + f 0 (x)h +1
1 00
f (x)h2 + f 000 (x)h3 + · · ·2!
3!

(2.19)のような Taylor 展開を使えばよいが，Taylor 展開を使わずに，いくつかの場所で評価した F (x, y) を混ぜ合わせて同じ 結果を得る方法が Runge-Kutta 法である．本研究で6

は次式のような「 4 次の 16 公式」を使用した．

1~y = y~0 + (k~1 + 2k~2 + 2k~3 + k~4 ) + O(h5 )6~k1 = f~(x0 , y~0 )hh
k~2 = f (x0 + , y~0 + k~1 )22hh
k~3 = f~(x0 + , y~0 + k~2 )22~~
k4 = f (x0 + h, y~0 + hk~3 ).

(2.20)

これは結果の誤差が h の 4 乗に比例する．(2.20) 式を| ψ(x, t) i =∞∑ψn (t)| n i(2.21)

n=1を求めるのに利用するには，

x = t, x0 = t0

(2.22)~y = (ψ1 (t), ψ2 (t), ψ3 (t), · · ·)(2.23)のように対応させる．

7

第 3 章 調和振動子基底
調和振動子では粒子は外場（ポテンシャルエネルギー）によって空間のある領域に局在した束縛状態にある．ここでは，一次元調和振動子の振る舞いについて要点を整理して示し，重要な演算子を具体的に行列表示で表す．

3.1

一次元調和振動子

x 軸上で原点 O からの距離に比例する引力F (x) = −kx = −mω 2 x

(3.1)

1
V (x) = mω 2 x2
2

(3.2)

は，ポテンシャルから導かれる．このような力を受けている粒子のハミルトニアンはĤ = −

h̄2 ∂ 2
1
+ mω 2 x2
2
2m ∂x
2

(3.3)

と記述される．したがって，この粒子に対する時間に依存しない Schrödinger 方程式は}{

1h̄2 ∂ 2+ mω 2 x2 | ϕ(x) i = ε| ϕ(x) i
−
2
2m ∂x
2

(3.4)

と表記される．ここで，式を簡略化し，計算を容易にするために無次元の変数
√

ξ=

mω
x
= ,
h̄
b

(3.5)

2ε
h̄ω

(3.6)

λ=
√

を導入すると（ b =
えられる．

h̄
は調和振動子長と呼ばれる），(3.4) 式は，次のように書き換
mω

(

)

∂2
− 2 + ξ 2 | ϕ(ξ) i = λ| ϕ(ξ) i
∂ξ

(3.7)

ある基底 | n i における状態ベクトル | ϕ(ξ) i の波動関数 ϕn (ξ) は，

h n | ϕ(ξ) i
8

(3.8)

(

2

)

で得られる．この波動関数は，エルミート多項式と exp − ξ2 の積に比例する形をとる．
()

ξ2ϕn (ξ) ∝ Hn (ξ) exp −2
λn = 2n + 1, n = 0, 1, 2, · · · .

(3.9)

エルミート多項式 Hn (ξ) はdnexp(−x2 )ndx(3.10)

Hn+1 (ξ) = 2ξHn (ξ) − 2nHn−1 ,(3.11)

H0 = 1,(3.12)

H1 (ξ) = 2ξ,(3.13)Hn (x) = (−1)n exp(x2 )によって定義され，

(

∫ ∞
−∞

)

∂2
∂
− 2ξ
+ 2n Hn (ξ) = 0,
2
∂ξ
∂ξ
{
2

Hn (ξ)Hm (ξ) exp(−ξ )dξ =(3.14)

0n 6= m√n
2 ! π n=m(3.15)などの関係を満たす関数の一群である．
変数を x と λ にもどし，(3.15) によって規格化した関数をつくると，√ψn (x) = 

2mω/h
2n n!

1/2


(√

Hn

1
εn = (n + )h̄ω,
2

)

(

)

mω
mω 2
x exp −x ,h̄2h̄n = 0, 1, 2, · · ·(3.16)(3.17)

が得られる．
ポテンシャルエネルギーが，(3.2) で表されるとき，粒子は x = 0 のまわりに局在する（ 図 3.1 参照）．波動関数 ϕn (ξ) は状態 n で n 個のノード を持ち，エネルギー固有値は h̄ω のとびとびの値を持つ．

9図 3.1: 一次元調和振動子の定常状態での波動関数

3.2

調和振動子の行列表示

一次元調和振動子のハミルトニアン(h̄2 ∂ 2
1
h̄ω
∂2
2 2
Ĥ = −
+
mω
x
=
−
+ ξ2
2m ∂x2 22∂ξ 2)(3.18)は，

1
Ĥ| n i = h̄ω(n + )| n i2(3.19)に従う Ĥ の固有ベクトル | n i を基底とする表現で，次のような対角行列で表わされる．

√

Ĥ =




h̄ 

2mω 




1
h̄ω
2

0
0
0
..
.

0
3
h̄ω
2
0
0
..
.

0
0
5
h̄ω
2
0
..
.

0
0
0
7
h̄ω
2
..
.



···

··· 


··· 
.
··· 

.. 
.

(3.20)

上昇演算子（ 生成演算子）a† ，下降演算子（ 消滅演算子）a を(∂
1
a = √ ξ−
∂ξ
2

)

†

(

1
∂
a= √ ξ+
∂ξ
2

√

=
)

√

=

(

)

mω
h̄ ∂
x−
,
2h̄
mω ∂x
(

h̄ ∂
mω
x+
2h̄
mω ∂x

10

(3.21))(3.22)

と書くと，基底 | n i はa† | n i =√a| n i =n + 1| n + 1 i,√

(3.23)

n| n − 1 i(3.24)

により構成される．これ等の等式の両辺と | m i との内積をとると，√
h m |a† | n i = n + 1δm,n+1 ,h m |a| n i =√nδm,n−1(3.25)(3.26)となり a† | n i と a| n i の行列要素が得られる．したがって a† と a の行列表記は次のよ
うになる．




0
0
0
0
 √
 1 0
0
0

√

†

2 0
0
a = 0
√

0
3 0
 0

..
..
.. . .
.
.
.
.

···

··· 


··· 
,
··· 

.. 
.

√





0

 0


a=
 0

 0

..
.

(3.27)

1 0
0 ···
√

0
2 0 ··· √003 ··· .
000 ··· 
..
..
. . . ..
.
.
.

(3.28)∂
は(3.21)，(3.22) より ξ と ∂ξ)
1 (
ξ = √ a + a† ,
2

)
∂
1 (
= √ a − a†∂ξ2

(3.29)∂と表され，(3.5) より x と ∂xも a† と a で次のように表される．

√

x=
∂
=
∂x

)
h̄ (
a + a† ,2mω

√(3.30))mω (a − a† .
2h̄∂
このようにして，下記の x と ∂xの行列表示を得る．

√
1 0
0
0
√
 √

√
2 0
 1 √0
√
h̄ 
 0
2 0
3
x=

√
2mω 
0
3 0
 0

..
..
..
..
.
.
.
.

11

(3.31)



···

··· 


··· 
,
··· 

.. 
.

(3.32)


√
1
0
0 ···
0
√
 √

 − 1
0
2
0 ··· 
√


√
√

∂
mω 

.
0
−
2
0
3
·
·
·
=


√
∂x
2h̄ 
0
− 3 0 ··· 
 0



..
..
..
.
... .
.
.
.
.


(3.33)∑

xn の行列要素 h i |xn | j i は，完全性 kn | kn ih kn | = 1 を利用して，以下のように x1の行列表示の積として求めることができ，結果として xn を簡単な疎行列として表すことができる．
∑

h i |xn | j i =h i |x| k1 ih k1 |x| k2 ih k2 |x| · · · | kn−1 ih kn−1 |x| j i.

(3.34)

k1 ,k2 ···kn−1
2

∂
0
1
2
3
4
5
6
本研究で使用する， ∂x
2 ，x ，x ，x ，x ，x ，x ，x の行列要素を以下に示す．

√

(
)
∂2
mω  n(n − 1)1
hi| 2|j i =δi,j−2 − n +δi,j +
∂x
h̄ 
2
2

√

(n + 1)(n + 2)

δi,j+2 ,

2(3.35)
h i |1| j i = δi,j ,(h i |x| j i =(3.36)h̄mω

(h i |x2 | j i =)1/2 √ n

h i |x3 | j i =



2
√

3(n + 1)+2
h i |x4 | j i =δi,j−1 +

2




n+1
δi,j+1  ,
2
(

h̄
mω

√

)

(n + 1)(n + 2)2

δi,j−3 +√n+1
1
δi,j+1 +
2
2

3n
2

√

n
δi,j−12(

δi,j−4 + n −

(3.39))1 √n(n − 1)δi,j−2
2

)

)
3(
3 √
+ n2 + (n + 1)2 δi,j + n −(n + 1)(n + 2)δi,j+2
4√2(n + 1)(n + 2)(n + 3)(n + 4)+
δi,j+4 ,412

δi,j+2 ,2(n + 1)(n + 2)(n + 3)
δi,j+3  ,
2

4
(




(3.38))4/2  √ n(n − 1)(n − 2)(n − 3)

(3.37)

1
δi,j−2 + n +δi,j +
2)3/2  √ 1 n(n − 1)(n − 2)
h̄

mω

(

2

)2/2  √ n(n − 1)
h̄

mω
(



√

(3.40)

(h i |x5 | j i =h̄mω)5/2  √ 1 n(n − 1)(n − 2)(n − 3)(n − 4)
4√5(n − 1)+4

2n(n − 1)(n − 2)5(2n2 + 1)
δi,j−3 +2
4

5(2n2 + 4n + 1)
+
4
√

1
+
4
(

h i |x | j i =6δi,j−5

√

√

√

n+1
5(n + 2)
δi,j+1 +2
4nδi,j−1
2(n + 1)(n + 2)(n + 3)
δi,j+32(n + 1)(n + 2)(n + 3)(n + 4)(n + 5)
δi,j+5 ,
2(3.41))6/2 {h̄1√n(n − 1)(n − 2)(n − 3)(n − 4)(n − 5)δi,j−6mω86n − 9 √15(n2 − n + 1) √+n(n − 1)(n − 2)(n − 3)δi,j−4 +n(n − 1)δi,j−2885(4n3 + 6n2 + 8n + 3)
15(n2 + 3n + 3) √+(n + 1)(n + 2)δi,j+2
δi,j +
8
8
6n + 15 √+(n + 1)(n + 2)(n + 3)(n + 4)δi,j+4
8}
1√(n + 1)(n + 2)(n + 3)(n + 4)(n + 5)(n + 6)δi,j+6 .
+(3.42)
8(3.35)〜(3.42) 式が間違いないことは，(3.34) 式と数値的に一致する結果を与えることにより確認済である．

13

第 4 章 数値計算の精度
数値計算に用いたパラメータと精度の関係について論じる．第 4〜5 章では，h̄ = m =1 となる単位系を用いる．

4.1

基底の切断（ truncation ）本来ならば，基底の状態番号 n は n = 0, 1, 2 · · · ∞ と無限大まで取り得る．しかし数値計算で無限大まで考えることはできないため，本研究では n の最大値 nmax が充分大
きいが有限であるとして計算を行う．
テスト用のポテンシャルエネルギーとして(V (x) =

)

6
x
− 1.0
31.6

(4.1)

をとり，初期状態を t = 0 で x = 0 付近に局在した波束とする（ 次章で初期状態の作り方を記述してある）．これを実時間発展させ，ある時刻 t(ここでは t = 10.0) でのエネ
ルギー期待値 E の値が変化しなくなるまで（ 差が充分小さくなるまで ）nmax の値を増やしていき，計算精度上，充分大きいといえる．nmax の値を決める．図 4.1 に t = 10.0での nmax = 5000 のエネルギー期待値との差を nmax に対してプロットした．図 4.1 より nmax は 1000 で充分であると判断し，本研究ではnmax = 1000 とした．調和振動子の√h̄
パラメータとして，調和振動子長（ b = mω ）がある．n の最大値を nmax としたときの，ポテンシャルの x 座標の広がり xmax は，調和振動子長 b を用いて√
xmax = b 2nmax + 1(4.2)と表され，確保したい x 座標の広がりを得るのに，充分大きな nmax と b が必要になる．
これについては，(4.1) 式のポテンシャルに対し ，nmax = 1000 を用いる場合には，値として b = 1.0 が計算精度上充分大きいといえることを確認した．

4.2

誤差の微小時間に対する依存性

Euler 法では nmax が大きいと，dt を小さくとらなければならないことが知られている．4 次の Runge-Kutta 法においても同様の傾向があることが予想される．
dt の値が充分小さければ，計算による誤差は無視できるほど 小さい（ 表 4.1 参照）．本
研究では，微小時間 dt = 1.0 × 10−6 で，充分な精度を与えるとみなしている．

14

0.0001
E_nmax-E_5000(at t=10.0)

E_nmax-E_5000

1e-06

1e-08

1e-10

1e-12

1e-14

1e-16

0

500

1000

1500

2000

2500 3000
n_max

3500

4000

4500

5000

図 4.1: t = 10.0 における nmax のエネルギー期待値と nmax = 5000 のエネルギー期待値との差

dt
1.0 × 10−3
1.0 × 10−4
1.0 × 10−5
1.0 × 10−6

t = 10.0 でのエネルギー期待値 E
1.1374690988184457
1.1374690988184430
1.1374690988184470
1.1374690988184375

表 4.1: 微小時間 dt と t = 10.0 でのエネルギー期待値 E

15

エネルギー期待値の保存

4.3

ある物理量 A を表す演算子 Â とハミルトニアン Ĥ が交換するときには，(ただし Â
が時間に依存する場合を除く）物理量 A は保存される．
Ĥ と Ĥ の交換関係は[Ĥ, Ĥ] = 0(4.3)となり，Ĥ が時間に依存しなければエネルギーは保存される．したがって，実時間発展においてエネルギーの期待値は一定である．
図 4.2 は，ポテンシャルエネルギーを (4.1) 式のものにとり，nmax = 1000 での実時
間発展によるエネルギー期待値の変動を時間に対してプロットしたものであり，一定に保たれていることが分かる．しかし ，あくまで数値計算であるため実は完全に一定というわけではない．t = 0 での値（ E(0) = 1.1374690988184504 ）からのずれを時間の関数としてプロットしたものを図 4.3 に示す．この結果から t = 0 での値からのずれは無視できるほど 小さく，数値的実時間発展におけるエネルギー期待値は高い精度で時間に依存せず一定とみなせることがわかる．
2

E

E

1.5

1

0.5

0

0

2

4

6

8

t

図 4.2: 実時間発展におけるエネルギー期待値の時間変化

1610

2.5e-14E(t)-E(0)E(t)-E(0)2e-14

1.5e-14

1e-14

5e-15

0

0

2

4

6

8

10

t

図 4.3: ある時刻 t でのエネルギー期待値の t = 0 でのエネルギー期待値からのずれ

4.4

ノルムの保存実時間発展ではノルムは保存されるはずである．
図 4.4 は，時刻 t での規格化前の状態ベクトルのノルムを，時刻 t = 0 での状態ベク
トルのノルムで割ったものをプロットしたものである．数値計算の誤差によりノルムは完全に保存されているわけではないが，その誤差は非常に小さく，ノルムは保存されていると考えてよい．
図 4.4 のそれぞれの値の 1 からのずれを図 4.5 に示す．数値計算により規格化された
値の 1 からのずれは非常に小さく，充分規格化できていると考えてよい．なお，図 4.4，
4.5 に示した以外の計算では，Runge-Kutta 法を 1 回行うごとに，状態ベクトルの再規格化をして，この微小なノルムのずれを修正している．

17

norm

1.4

norm(t)/norm(0)

1.2
1
0.8
0.6
0.4
0.2
0

0

2

4

6

8

10

t

図 4.4: 時刻 t での状態ベクトルのノルムを時刻 t = 0 での状態ベクトルのノルムで割っ
たもの
4e-14
delta_norm
3.5e-14|norm(t)/norm(0)-1.0|

3e-14
2.5e-142e-14
1.5e-14
1e-14
5e-15
0

0

2

4

6
t

図 4.5: 図 4.4 のそれぞれの値の 1 からのずれ

18

8

10

冷却（ 虚時間発展）により得た固有状態の精度

4.5

虚時間発展ではエネルギーの期待値は時間とともに減少し ，基底状態のエネルギーに収束する．これにより虚時間発展で任意の初期状態を固有状態へと冷却することができる．
図 4.6，4.7 は，ポテンシャルエネルギーを1
V (x) = x2
2

(4.4)

として nmax = 1000 での虚時間発展におけるエネルギー期待値を時間に対してプロットしたものである．図 4.8 は，時刻 t でのエネルギー期待値と t = 100.0 でのエネル
ギー期待値との差を対数プロットしたものである．t = 40.0 付近から完全に一定値（ E = 0.249987 · · · ' h̄ω= 0.25 ）となり，e−t に比例して基底状態に収束する．本研究2では t = 10.0 で充分収束できているとみなす．
1.4

E
hbar*omega/2

1.2

E

1

0.8

0.6

0.4

0.2

0

2

4

6

8

t

図 4.6: 虚時間発展におけるエネルギー期待値の時間変化 (t = 0〜10)

19

10

1.4

E
hbar*omega/2

1.2

E

1

0.8

0.6

0.4

0.2

0

20

40

60

80

100

t

図 4.7: 虚時間発展におけるエネルギー期待値の時間変化 (t = 0〜100)
100000E(t)-E(100)exp(-x)1

E(t)-E(100)

1e-05
1e-10
1e-15
1e-20
1e-25
1e-30

0

10

20

30
t

40

50

図 4.8: t = 100 でのエネルギー期待値からのずれの時間変化

2060

第 5 章 波束の時間発展
この章では，2 つのポテンシャル V (x) を例にとり，波動関数の時間発展を求める．
まず，初期状態を準備しなければならない．本研究では，本来のポテンシャルとは異なるポテンシャル V 0 (x) 中で状態を虚時間発展させて V 0 (x) の基底状態へと冷却させることにより初期状態を作る．次に，このようにして準備した初期状態の波束を本来のポテンシャル V (x) 中で実時間発展させ，その振る舞いを調べる．

5.1

非調和振動子

まず，2 次のポテンシャル V 0 (x) 中で虚時間発展させ，x = 0 の付近に初期状態の波束を準備する．
(0

V (x) =

x
√
8

)2

(5.1)

次に，箱のような 6 次のポテンシャル V (x) を用意し，準備した初期状態の波束を実時間発展させる．
(x
V (x) =
31.6

)6

(5.2)

図 5.1 は，ポテンシャル (5.2)(5.1) およびエネルギーの期待値 E である．図 5.2 は，時
刻 t = 0, 10, 20, 40, 60, 80, 99.9 での波動関数の実部と虚部である．
古典力学では，ポテンシャルの底 (x = 0) に静止したままだが，量子力学では波束が拡散してゆくことがわかる．これは位置と運動量の不確定性関係によるものとして理
解できる．

21

5
V(x)
V’(x)
E

4

V(x)

3

2

1

0
-40

-20

0
t

20

40

図 5.1: 実時間および虚時間のポテンシャル
6
Re part
Im partt=0.0
5
t=20.04

psi(x)

t=40.03
t=60.021
0-1

t=80.0t=99.9

-40

-20

0
x

図 5.2: 波束の実時間発展

22

20

40

今度は，初期状態を x = 0 からずらしたところにつくり，同様に箱型の 6 次のポテンシャルで実時間発展させる．
実時間でのポテンシャル V (x)，初期状態の波束を準備するための虚時間でのポテン
シャル V 0 (x) はそれぞれ以下のとおりである．
(

)

x 6
V (x) =
31.6(
)2
x − 300√
V (x) =
8

(5.3)
(5.4)

図 5.3 は，ポテンシャル (5.3)(5.4) およびエネルギーの期待値 E である．図 5.4 は，時
刻 t = 0, 10, 20, 40, 60, 80, 99.9 での波動関数の実部と虚部である．
ポテンシャルエネルギーの高い地点におかれた波束は，拡散するだけでなく，古典的粒子と同様にポテンシャルの斜面をすべり落ち，ポテンシャルの平らな底をよこぎり，左側のポテンシャル壁に衝突してはね返されることが見てとれる．これは古典力
学が量子力学の近似として実現される例と考えることができる．

23

10
V(x)
V’(x)
E

8

V(x)

6

4

2

0

-2

-40

-20

0
t

20

40

図 5.3: 実時間および虚時間のポテンシャル
6
t=0.0
5
t=20.0Re part4

Im partpsi(x)

t=40.03
t=60.021
0-1

t=80.0t=99.9

-40

-20

0
x

図 5.4: 波束の実時間発展

24

20

40

5.2 トンネル効果
次に，実時間において二重の井戸型のポテンシャルを考える．初期状態をつくるた
めの虚時間のポテンシャル V 0 (x) は，x > 0 の井戸に局在した初期状態となるように，
x < 0 のポテンシャルを埋め立てるような，実時間のポテンシャルに 6 次の項を加えた形にしてある．
実時間でのポテンシャル V (x)，初期状態の波束を準備するための虚時間でのポテン
シャル V 0 (x) はそれぞれ以下のとおりである．

1 2
1 4
x +
x
10
200
)
(
x − 10 6
0V (x) = V (x) +24.4
V (x) = −

(5.5)
(5.6)

図 5.5 は，ポテンシャル (5.5)(5.6) およびエネルギーの期待値 E である．図 5.6 は，時
刻 t = 0, 10, 20, 40, 60, 80, 99.9 での波動関数の実部と虚部である．
エネルギーの期待値がポテンシャル障壁の最高点より低いので，古典的には粒子は左側のポテンシャル井戸へしみ出すことはできないはずだが，量子力学ではトンネリングによりしみ出してゆくことがわかる．

25

3

V(x)
V’(x)
E

2.5
2

V(x)

1.5
1
0.5
0
-0.5
-1
-20

-15

-10

-5

0
x

5

10

15

20

図 5.5: 実時間および虚時間のポテンシャル
6
Re partt=0.0
5
Im partt=20.04

psi(x)

t=40.03
t=60.0
2
1
0
-1
-20

t=80.0t=99.9

-15

-10

-5

0
x

図 5.6: 波束の時間発展

26

5

10

15

20

今度は，二重井戸の障壁を高くし ，同様に x > 0 の井戸に初期状態をつくり実時間発展させる．実時間でのポテンシャル V (x)，初期状態の波束を準備するための虚時間
でのポテンシャル V 0 (x) はそれぞれ以下のとおりである．

2
1
V (x) = − x2 + x4
550()
x − 10 6
V 0 (x) = V (x) +
17.8

(5.7)
(5.8)

図 5.7 は，ポテンシャル (5.7)(5.8) およびエネルギーの期待値 E である．図 5.8 は，時
刻 t = 0, 10, 20, 40, 60, 80, 99.9 での波動関数の実部と虚部である．
今度は，トンネリングが非常に遅くなり，プロットした時間範囲では，左側のポテ
ンシャル井戸へのしみ出しは見られない．また，右側の井戸中の波束の運動は，全体が一斉に位相 e−iωt で変化する定常状態によく似ている．このような右側の井戸中の波束は準定常状態と呼ぶ．

27

3
V(x)
V’(x)
E

2

V(x)

1

0

-1

-2
-20

-15

-10

-5

0
x

5

10

15

20

図 5.7: 実時間および虚時間のポテンシャル
65

psi(x)

4
3
2

Re partt=0.0

Im partt=20.0t=40.0t=60.0t=80.01
t=99.9
0
-1
-20

-15

-10

-5

0
x

5

図 5.8: 波束の実時間発展

28

10

15

20

第 6 章 結論
時間に依存する Schrödinger 方程式の調和振動子基底展開による表現を求め，それを数値的に解く C 言語プログラムを作製した (プログラムリストを付録に収録した)．また，本研究を進めるにあたって，量子力学の基礎を系統的に学び直した．
応用として，波束の実時間発展の量子力学的な側面と，古典力学と類似性のある側面を示した．これは 4 本の動画として公開しており，参考文献にその URL を記載してある．
まず，箱型のポテンシャルでの実時間発展について，x = 0 付近に初期状態の波束を準備する例では，実時間発展により波束が拡散していく様子を確認した．これは明ら
かに量子力学的な振る舞いを示している．
x = 0 からずらしたところ（ x = 30 ）に初期状態の波束を準備する例では，実時間発
展により波束が拡散しながらも，ポテンシャルにそって振動を繰り返した．量子力学
的な振る舞いと古典力学と類似性のある振る舞いが混在した様子を確認した．
実時間発展によるトンネル効果について，二重井戸のポテンシャルの障壁高が低い
場合には，波動関数のすそが x < 0 までのびているため，実時間発展ですぐに x < 0 の井戸の成分が混ざ った．ここでは量子力学特有の現象であるトンネル効果の様子が確認された．
障壁高をある程度高くすると，波動関数のすそが x < 0 までのびることができず，ト
ンネル効果は見られない．波束は x > 0 の井戸に局在し ，実時間発展で固有状態のよ
うに全体が同位相 e−iωt で振動しているように見える．ただし，x < 0 の井戸へのトンネリングが起こりにくいため，固有状態ではなく準定常状態といえる．
本研究では，巨大な次元の調和振動子基底で一次元量子系を数値的に扱うことができた．今後は，量子複雑系の例として，高次元空間の振動運動に拡張することが発展課題として考えられる．

29

参考文献
[1] 清水明：量子論の基礎 - その本質のやさしい理解のために (サイエンス社, 2003)[2] 小出昭一郎：量子力学 (I) (裳華房, 2010)[3] Leonard I. Schiﬀ：新版 量子力学 上 (吉岡書店, 2008)[4] 砂川重信：量子力学の考え方 (岩波書店, 1993)[5] 奥村晴彦：C 言語による最新アルゴ リズム事典 (技術評論社, 1992)[6] 筑波大学 課題探求実習レポート 拡散モンテカルロ法による軽い原子基底状態の数値計算
http://www.tac.tsukuba.ac.jp/ ozawa/Sato.pdf[7] 東京学芸大学 資料「昇降演算子法による調和振動子」 新田英雄
http://www.u-gakugei.ac.jp/ nitta/operator.pdf[8] 波束の時間発展
動画 1
動画 2
動画 3
動画 4

http://apphy.u-fukui.ac.jp/~nucleus/daiThesis12_movie1.m1v
http://apphy.u-fukui.ac.jp/~nucleus/daiThesis12_movie2.m1v
http://apphy.u-fukui.ac.jp/~nucleus/daiThesis12_movie3.m1v
http://apphy.u-fukui.ac.jp/~nucleus/daiThesis12_movie4.m1v

30

謝辞
本論文を作成するにあたり、田嶋直樹先生には終始丁寧なご 指導をしていただいたことに感謝し 、お礼申し上げます。また林明久先生にも本研究及び日常的なことにおいても、実に丁寧な指導、お世話をしていただきました。
本研究に対してご 意見をいただいた、多くの物理工学科の先生方をはじめ、諸先輩方や同卒研生にもお礼申し上げ、謝辞の言葉とさせていただきます。

31

付録

プログラムソースコード

以下のプログラムでは概ね以下のような処理が行われる．

• 状態ベクトルの初期値を設定• Runge-Kutta 法で虚時間発展（冷却）させ，何らかのハミルトニアンの基底状態として，初期状態の波束を準備• Runge-Kutta 法で，準備した初期状態の波束を実時間発展させる• エネルギー期待値 E ，エネルギーのゆらぎ ∆E ，x の期待値，x のゆらぎ ∆x を
出力
• 状態ベクトルおよび波動関数を出力# include <stdio.h># include <stdlib.h># include <complex.h># include <math.h># define nmax 1000
# define N (nmax+1)# define maxpower 6// 量子数 n の最大値// 基底に含める（ 量子数 n=0..nmax をもつ）基底の個数// x のべきの最大値

const double hbar = 1.0;
const int ki = 2*maxpower;
double ol_b; // b (oscillator length of the basis)
int runge(double tmax,int itmax,int dit_obs, int dit_vec,
int dit_wav,int cooling,double complex *v);
int set_hamiltonian_matrix(int option);
double deriv_xi_2 (int i, int j);
double xi_0 (int i, int j);
double xi_1 (int i, int j);
double xi_2 (int i, int j);
double xi_3 (int i, int j);
double xi_4 (int i, int j);
double xi_5 (int i, int j);
double xi_6 (int i, int j);
double complex timederivative( int i, double t, double complex *x);
double norm_square(double complex *x);
int normalize(double complex *x);
double harmonic_oscillator_wavefunction(int n, double x);
int write_vector(double time, double complex *v);

32

int write_wavefunction(double time, double complex *v);
int write_potential(double q0, double q1, double q2, double q3,
double q4, double q5, double q6);
double deriv_xi_2 (int i, int j);
double xi_2 (int i, int j);
double energy_expectation_value( double complex *x);
double energy_2_expectation_value( double complex *x);
double x_expectation_value( double complex *x);
double x_2_expectation_value( double complex *x);
int main(){// main1();
main2();// main3();
return 0;}
int main1(){double complex v[N]; // 状態ベクトル
double tmax; // 時間発展を行う時間の長さ
int itmax = 10000; // 時間の刻みの総数（ 0 番..itmax 番の刻みがある)
int dit_obs = 10; // observables を出力する時間ステップ番号の間隔
int dit_vec = 10; // 状態ベクトルを出力する時間ステップ番号の間隔
int dit_wav = 100; // 波動関数をファイル出力する時間ステップ番号の間隔// 注：dit_vec,dit_wav ともにゼロ以下の値を設定すると出力が完全に抑止される。
double omega = 0.5; // 注意：同一用途の同名変数の値を２箇所で定義している。
double beta = 0.1; // 初期状態設定用パラメータ// （ 温度の逆数、単位＝ 1/基底準位間隔）
int cooling = 0; // =0 で実時間発展、=1 で虚時間発展
int i; // 作業用変数 (カウンター)
set_hamiltonian_matrix(1);tmax = (4.0*M_PI/omega)*1;/* 初期状態の設定 */for ( i=0; i<N; i++){// v[i] = 1.0;v[i] = exp(-beta*i);}
normalize(v);// 状態の規格化/* 虚時間発展 : 何らかのハミルトニアンの基底状態として初期状態を作る */// runge(tmax, itmax, 0, 0, 0, 1, v);// runge(tmax, itmax, dit_obs, dit_vec, dit_wav, 1, v);/* 実時間発展 */runge(tmax, itmax, dit_obs, dit_vec, dit_wav, cooling, v);
return 0;}
int main2(){double complex v[N]; // 状態ベクトル33

double tmax; // 時間発展を行う時間の長さ
int itmax; // 時間の刻みの総数（ 0 番..itmax 番の刻みがある)
int dit_obs=0; // observables を出力する時間ステップ番号の間隔
int dit_vec=0; // 状態ベクトルを出力する時間ステップ番号の間隔
int dit_wav=0; // 波動関数をファイル出力する時間ステップ番号の間隔// 注：dit_vec,dit_wav ともにゼロ以下の値を設定すると出力が完全に抑止される。
double omega = 0.5; // 注意：同一用途の同名変数の値を２箇所で定義している。
double beta = 0.1; // 初期状態設定用パラメータ// （ 温度の逆数、単位＝ 1/基底準位間隔）
int cooling = 1; // =0 で実時間発展、=1 で虚時間発展
int i; // 作業用変数 (カウンター)/* 初期状態の設定 */for ( i=0; i<N; i++){// v[i] = 1.0;v[i] = exp(-beta*i);}
normalize(v);// 状態の規格化/* 虚時間発展 : 何らかのハミルトニアンの基底状態として初期状態を作る */set_hamiltonian_matrix(1);tmax = 10.0;
itmax=10000000;
dit_obs = 0;
dit_vec = 0;
dit_wav = 0;
cooling = 1;
runge(tmax, itmax, dit_obs, dit_vec, dit_wav, cooling, v);/* 実時間発展 */set_hamiltonian_matrix(2) ;tmax = 100.0;
itmax= 100000000;
dit_obs = 100000;
dit_vec = 100000;
dit_wav = 100000;
cooling = 0;
runge(tmax, itmax, dit_obs, dit_vec, dit_wav, cooling, v);
return 0;}/* ルンゲ・クッタ法による時間発展 */int runge(double tmax, int itmax, int dit_obs, int dit_vec,
int dit_wav, int cooling, double complex *v){
int nprint=6; // 出力する状態ベクトルの成分の個数（ n=0..nprint-1 を出力）
double complex vtmp[N], k0[N], k1[N], k2[N], k3[N];double t,dt; // 時刻, 時刻の刻み幅
double complex tfct; // 実時間発展では dt, 虚時間発展では -I*dt
double norm2; // 状態ベクトルのノルムの２乗
double fct; // 作業変数（ 因子）double E; // エネルギー期待値 E

34

double E_2; // E の 2 乗の期待値
double X; // x の期待値
double X_2; // x の 2 乗の期待値
double delta_E; // Δ E
double delta_x; // Δ x
int itime; // 時刻の刻みの番号=0 番..itmax 番
int i; // 作業変数（ カウンター)
dt = tmax/itmax;
if(nprint>N)nprint=N;/* t=0 */t = 0.0;
if(dit_obs > 0){ // t=0 での observables の出力
E = energy_expectation_value(v);E_2 = energy_2_expectation_value(v);X = x_expectation_value(v);
X_2 = x_2_expectation_value(v);delta_E = sqrt(fabs(E_2 - E*E));
delta_x = sqrt(fabs(X_2 - X*X));
printf("# t, E, delta_E, X, delta_x\n");printf("%f %.16f %.16f %.16f %.16f :obs\n", t, E, delta_E, X, delta_x);}
if(dit_vec > 0) write_vector(t,v); // t=0 での vector の出力
if(dit_wav > 0) write_wavefunction(t,v);/* t>0 */if(cooling==0){ tfct=dt; } else { tfct=-dt*I; }for ( itime=1; itime<=itmax; itime++){t = dt*itime;for ( i=0; i<N; i++) k0[i] = timederivative( i, t, v)*tfct;
for ( i=0; i<N; i++) vtmp[i]=v[i] + k0[i]/2.0;
for ( i=0; i<N; i++) k1[i] = timederivative( i, t + dt/2.0, vtmp)*tfct;
for ( i=0; i<N; i++) vtmp[i]=v[i] + k1[i]/2.0;
for ( i=0; i<N; i++) k2[i] = timederivative( i, t + dt/2.0, vtmp)*tfct;
for ( i=0; i<N; i++) vtmp[i]=v[i] + k2[i];for ( i=0; i<N; i++) k3[i] = timederivative( i, t + dt, vtmp)*tfct;
for ( i=0; i<N; i++) v[i] += ( k0[i] + 2.0*k1[i] + 2.0*k2[i] + k3[i])/6.0;

if(dit_obs > 0 && itime % dit_obs == 0){
E = energy_expectation_value(v);E_2 = energy_2_expectation_value(v);X = x_expectation_value(v);
X_2 = x_2_expectation_value(v);delta_E = sqrt(fabs(E_2 - E*E));
delta_x = sqrt(fabs(X_2 - X*X));
printf("%f %.16f %.16f %.16f %.16f :obs\n", t, E, delta_E, X, delta_x);35}normalize(v); // 数値計算誤差によるノルムの１からのずれを修正する
if(dit_vec > 0 && itime % dit_vec == 0) write_vector(t,v);
if(dit_wav > 0 && itime % dit_wav == 0 ) write_wavefunction(t,v);}
if(dit_wav < 0 ) write_wavefunction(t,v);
return 0;}/* hamiltonian の計算 */// double ham1[N][N];double *ham1 = NULL;
int set_hamiltonian_matrix(int option){
double mass_b = 1.0;
double omega_b = 1.0;
double mass = 1.0;
double omega = 0.5;
double x0 = 30.0;
double x1 = 30.0;
double hight1 = 1000.0;
double hight2 = 2.0;
double p, q1, q2, q3, q4, q5, q6, q0;
int i, j;double x0_pot=10.0; // location of minimum
double h_pot=0.5; // barrier height
double d_pot=0.0; // 0.25; // splitting of the depth of two minima
if(ham1 == NULL){
ham1=malloc(N*(2*maxpower + 1)*sizeof(double));if(ham1 == NULL) {fprintf(stderr,"malloc:%d\n",N);exit(1);}}ol_b = sqrt(hbar/(mass_b*omega_b));// 虚時間発展のハミルトニアンの各項の係数 (cooling により波束をつくる)
if( option == 1 ){p = -hbar*hbar/(2.0*mass_b*ol_b*ol_b); // 2 階微分の係数q1 = 0.0; // 1 次の項の係数
q2 = mass*omega*omega*ol_b*ol_b/2.0; // 2 次の項の係数
q3 = 0.0; // 3 次の項の係数
q4 = 0.0; // 4 次の項の係数
q5 = 0.0; // 5 次の項の係数
q6 = 0.0; // mass*omega_b*omega_b*ol_b*ol_b*ol_b*ol_b*ol_b*ol_b/2.0; // 6
次の項の係数q0 = 0.0; // mass*omega_b*omega_b*x0*x0*x0*x0*x0*x0/2.0; // 定数項}
else if( abs(option) == 2 ){// 実時間発展のハミルトニアン (準備した波束を振動させる)
double bxpow1=ol_b/x0_pot,bxpow=bxpow1;
p = -hbar*hbar/(2.0*mass*ol_b*ol_b);36

q0 = 0.0;
q1 = 0.0;
q2 = 0.0;
q3 = 0.0;
q4 = 0.0;
q5 = 0.0;
q6 = 1.0/995686217.814; // (1/31.6)^6;
if(option < 0){// 虚数時間発展のハミルトニアン (波束を準備する)
double vspp=h_pot/pow(2*x0_pot,6.0);
q0 += x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
q1 -= 6.0*x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
q2 += 15.0*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
q3 -= 20.0*x0_pot*x0_pot*x0_pot*vspp;
q4 += 15.0*x0_pot*x0_pot*vspp;
q5 = -6.0*x0_pot*vspp;
q6 = vspp;}}
for( i=0; i<N; i++){int js = i-maxpower;
if (js < 0) js = 0;// jstart

int je = i+maxpower; // jend
if (je > nmax) je = nmax;
for ( j=js; j<=je; j++){int k = ki*i + j;
ham1[k] = p*deriv_xi_2(i, j)+q0*xi_0(i, j)+q1*xi_1 (i, j)+q2*xi_2(i, j)+q3*xi_3(i, j)+q4*xi_4(i, j)+q5*xi_5(i, j)+q6*xi_6(i, j);}}printf("# %f =mass\n",mass);
printf("# %f =ol_b\n",ol_b);
write_potential(q0, q1, q2, q3, q4, q5, q6);
return 0;}// <i|(d^2/dxi^2)|j> 部分の計算
double deriv_xi_2 (int i, int j){if( i==(j-2) ){return sqrt(j*j-j)*0.5;}
else if( i==j ){return -j-0.5;}
else if( i==(j+2)){
return sqrt((j+1.0)*(j+2.0))*0.5;}
else {
return 0;}37}// <i|j> 部分の計算
double xi_0 (int i, int j){if( i==j ){return 1;}
else {
return 0;}}// <i|(xi^1)|j> 部分の計算
double xi_1 (int i, int j){if( i==(j-1) ){return sqrt(j*0.5);}
else if( i==(j+1)){
return sqrt((j+1.0)*0.5);}
else {
return 0;}}// <i|(xi^2)|j> 部分の計算
double xi_2 (int i, int j){if( i==(j-2) ){return sqrt(j*j-j)*0.5;}
else if( i==j ){return j+0.5;}
else if( i==(j+2)){
return sqrt((j+1.0)*(j+2.0))*0.5;}
else {
return 0;}}// <i|(xi^3)|j> 部分の計算
double xi_3 (int i, int j){if( i==(j-3) ){return sqrt(j*(j-1.0)*(j-2.0)*0.5)*0.5;}
else if( i==(j-1) ){return 3.0*j*0.5*sqrt(j*0.5);}
else if( i==(j+1) ){return 3.0*(j+1.0)*0.5*sqrt((j+1.0)*0.5);}
else if( i==(j+3) ){38

return 0.5*sqrt((j+1.0)*(j+2.0)*(j+3.0)*0.5);}
else {
return 0;}}// <i|(xi^4)|j> 部分の計算
double xi_4 (int i, int j){if( i==(j-4) ){return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0))*0.25;}
else if( i==(j-2) ){return (2.0*j-1.0)*0.5*sqrt(j*(j-1.0));}
else if( i==j ){return 3.0*0.25*(j*j+(j+1.0)*(j+1.0));}
else if( i==(j+2) ){return (2.0*j+3.0)*0.5*sqrt((j+1.0)*(j+2.0));}
else if( i==(j+4) ){return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0))*0.25;}
else {
return 0;}}// <i|(xi^5)|j> 部分の計算
double xi_5 (int i, int j){if( i==(j-5) ){return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0)*(j-4.0)/2.0)*0.25;}
else if( i==(j-3) ){return (5.0*j-5.0)*sqrt(j*(j-1.0)*(j-2.0)/2.0)*0.25;}
else if( i==(j-1) ){return (10.0*j*j+5.0)*sqrt(j/2.0)*0.25;}
else if( i==(j+1) ){return (10.0*j*j+20.0*j+15.0)*sqrt((j+1.0)/2.0)*0.25;}
else if( i==(j+3) ){return (5.0*j+10.0)*sqrt((j+1.0)*(j+2.0)*(j+3.0)/2.0)*0.25;}
else if( i==(j+5) ){return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0)*(j+5.0)/2.0)*0.25;}
else {
return 0;}}// <i|(xi^6)|j> 部分の計算

39

double xi_6 (int i, int j){if( i==(j-6) ){return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0)*(j-4.0)*(j-5.0))/8.0;}
else if( i==(j-4) ){return (6.0*j-9.0)*sqrt(j*(j-1.0)*(j-2.0)*(j-3.0))/8.0;}
else if( i==(j-2) ){return (15.0*j*j-15.0*j+15.0)*sqrt(j*(j-1.0))/8.0;}
else if( i==j ){return (20.0*j*j*j+30.0*j*j+40.0*j+15.0)/8.0;}
else if( i==(j+2) ){return (15.0*j*j+45.0*j+45.0)*sqrt((j+1.0)*(j+2.0))/8.0;}
else if( i==(j+4) ){return (6.0*j+15.0)*sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0))/8.0;}
else if( i==(j+6) ){return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0)*(j+5.0)*(j+6.0))/8.0;}
else {
return 0;}}/* 右辺 = sigma(hmn/i*hbar)*vn の計算 */double complex timederivative( int i, double t, double complex *x){
double complex s;
int j, js, je;
js = i-maxpower;
if (js < 0) js = 0;je = i+maxpower;
if (je > nmax) je = nmax;
s = 0;for ( j=js; j<=je; j++){s += ham1[ki*i+j]*x[j];}return -I/hbar*s;}
double norm_square(double complex *x){ // 状態ベクトルのノルムの２乗を計算
int i;
double s;s = 0.0;for (i=0; i<N; i++){s += creal(x[i])*creal(x[i]) + cimag(x[i])*cimag(x[i]);}return s;}

40

int normalize(double complex *x){ // 状態ベクトルを規格化
int i;
double fct;
fct = 1.0/sqrt(norm_square(x));
for ( i=0; i<N; i++){x[i] *= fct;}
return 0;}
double harmonic_oscillator_wavefunction(int n, double x){// y/m/d=13/1/17//
exp(-x^2/2) H_n(x)/sqrt(sqrt(pi) n^2 n!)// using H_n(x) = 2*x*H_{n-1}(x) -2*(n-1)*H_{n-2}(x)double a,b,c;double fnorm; // normalization factor
double y; // retuned value psi_n(x)// = value of the normalized n-node eigenwavefunction at x
int k, k1=0, k2=0;
double x2; // x*x
const double fct=1.0e50,fctreci=1.0/fct;x2=x*x;if(n>=2){b=1.0;
c=2.0*x;for(k=1;k<n;k++){
a=b;
b=c;
c=2*(x*b-k*a);if(fabs(c) > fct){k1++; a*=fctreci; b*=fctreci; c*=fctreci;}}
fnorm=1;
for(k=1;k<=n;k++){fnorm*=2*k;if(fabs(fnorm) > fct){k2++; fnorm*=fctreci;}}y=exp(-x2*0.5+(k1-k2*0.5)*log(fct))*c/sqrt(sqrt(M_PI)*fnorm);}
else if(n==1) y=exp(-x2*0.5)*2*x/sqrt(sqrt(M_PI)*2.0);
else if(n==0) y=exp(-x2*0.5)/sqrt(sqrt(M_PI));else { // if(n<0)
printf("harmonic_oscillator_wavefunction:argument error %d %f\n",n,x);
return 0.0;}
return y;}
int write_vector(double time, double complex *v){// 状態ベクトル (v) を名前に番号が付いたファイルに出力する。
//
time: 時刻. 参考情報として出力ファイルに書き加える。
//
v: 状態ベクトル（ 調和振動子基底での展開係数)41int i;static int filenumber = 0;
double s;
FILE *fo;char filename[256];
if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
snprintf(filename,256,"vec%03d.dat",filenumber++);fo=fopen(filename,"w");fprintf(fo,"# time=%.5f\n",time);
s=0.0;
for(i=0;i<=nmax;i++){fprintf(fo,"%3d %.9f %.9f\n",i,creal(v[i]),cimag(v[i]));s+=creal(v[i])*creal(v[i])+cimag(v[i])*cimag(v[i]);}fprintf(fo,"# squared norm=%.10f\n",s);fclose(fo);return 0;}
int write_wavefunction(double time, double complex *v){// 波動関数 (x) を名前に番号が付いたファイルに出力する。
//
time: 時刻. 参考情報として出力ファイルに書き加える。
//
v: 状態ベクトル（ 調和振動子基底での展開係数)
int npsi=500; // 波動関数の値を出力する x 軸上の点の個数-1 (0 番..npsi 番の点)
double xmin=-80.0; // 波動関数の値を出力する左端の点の座標
double xmax= 80.0; // 波動関数の値を出力する右端の点の座標
int i,n;static int filenumber = 0;
double complex s;
double norm2; // 状態のノルムの２乗（ x 表示で求めたもの)
double x,dx;
FILE *fo;char filename[256];
if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
snprintf(filename,256,"wf%03d.dat",filenumber++);fo=fopen(filename,"w");fprintf(fo,"# time=%.5f\n",time);
dx=(xmax-xmin)/npsi;
norm2=0.0;
for(i=0;i<=npsi;i++) {x=xmin+dx*i;
s=0.0;
for(n=0;n<=nmax;n++){s+=v[n]*harmonic_oscillator_wavefunction(n,x/ol_b);}s/=sqrt(ol_b);fprintf(fo,"%.3f %.6f %.6f\n",x,creal(s),cimag(s));
norm2+=creal(s)*creal(s)+cimag(s)*cimag(s);}fprintf(fo,"# squared norm=%.10f\n",norm2*dx);
fclose(fo);return 0;}
int write_potential(double q0, double q1, double q2, double q3,42double q4, double q5, double q6){// potential energy を名前に番号が付いたファイルに出力する。
int npsi=500; // potential の値を出力する x 軸上の点の個数-1 (0 番..npsi 番の点)
double xmin=-80.0; // potential の値を出力する左端の点の座標
double xmax= 80.0; // potential の値を出力する右端の点の座標
int i;
static int filenumber = 0;
double x,dx;
FILE *fo;char filename[256];
if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
snprintf(filename,256,"pot%03d.dat",filenumber++);fo=fopen(filename,"w");printf("# %f =q0\n",q0);
printf("# %f =q1\n",q1);
printf("# %f =q2\n",q2);
printf("# %f =q3\n",q3);
printf("# %f =q4\n",q4);
printf("# %f =q5\n",q5);
printf("# %f =q6\n",q6);
dx=(xmax-xmin)/npsi;
for(i=0;i<=npsi;i++) {x=xmin+dx*i;
fprintf(fo,"%.3f %.6f\n",x,q0+x*(q1+x*(q2+x*(q3+x*(q4+x*(q5+x*q6))))));}
fclose(fo);return 0;}// エネルギー期待値の計算
double energy_expectation_value( double complex *x){
double complex s;
int i, j, k, js, je;
s = 0;for ( i=0; i<N; i++){js = i-maxpower;
if (js < 0) js = 0;je = i+maxpower;
if (je > nmax) je = nmax;
k = ki*i+j;
for ( j=js; j<=je; j++){s += conj(x[i])*ham1[ki*i+j]*x[j];}}// s の虚部が 0 であることの確認
if (fabs(cimag(s)) > 1.0e-20){
printf("%f %f\n", creal(s), cimag(s));}
return creal(s);}43// エネルギーの 2 乗の期待値の計算
double energy_2_expectation_value( double complex *x){
double complex s, s1;
int i, j, k, js, je, ks, ke;
s = 0.0;for ( i=0; i<N; i++){js = i-maxpower*2;
if (js < 0) js = 0;je = i+maxpower*2;
if (je > nmax) je = nmax;
for ( j=js; j<=je; j++){if (i > j) ks = i-maxpower;
else ks = j-maxpower;
if (ks < 0) ks = 0;
if (i > j) ke = j+maxpower;
else ke = i+maxpower;
if (ke > nmax) ke = nmax;
s1 = 0.0;
for (k=ks; k<=ke; k++){s1 += ham1[ki*i+k]*ham1[ki*k+j];}s += s1*conj(x[i])*x[j];}}// s の虚部が 0 であることの確認
if (fabs(cimag(s)) > 1.0e-20){
printf("%f %f\n", creal(s), cimag(s));}
return creal(s);}// x の期待値の計算
double x_expectation_value( double complex *x){
double complex s;
int i, j;s = 0;for ( i=0; i<N; i++){for ( j=0; j<N; j++){s += conj(x[i])*x[j]*ol_b*xi_1 (i, j);}}// s の虚部が 0 であることの確認
if (fabs(cimag(s)) > 1.0e-20){
printf("%f %f\n", creal(s), cimag(s));}
return creal(s);}// x の 2 乗の期待値の計算

44

double x_2_expectation_value( double complex *x){
double complex s;
int i, j;s = 0;for ( i=0; i<N; i++){for ( j=0; j<N; j++){s += conj(x[i])*x[j]*ol_b*ol_b*xi_2 (i, j);}}// s の虚部が 0 であることの確認
if (fabs(cimag(s)) > 1.0e-20){
printf("%f %f\n", creal(s), cimag(s));}
return creal(s);}45