様々なポテンシャルにおける量子系の Metropolis アルゴリズムに基づくモンテカルロ計算
相木悠一, 近藤千尋, 堀敬一朗, 松尾一輝

2018 年 4 月 1 日

1目次
イントロダクション

3

理論物理と数値計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

経路積分法による量子力学へのアプローチとモンテカルロ計算

42.1

経路積分法の導出と虚時間形式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42.2

経路積分で基底状態における期待値を計算する方法. . . . . . . . . . . . . . . . . . . . . . .

7

2.3

Metropolis アルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8調和振動子ポテンシャル

10

3.1

1 次元調和振動子ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

3.2

3 次元調和振動子ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

3.3

調和振動子のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

Woods-Saxon ポテンシャル

22

4.1

殻模型と Woods-Saxon ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.2

ハイパー核と Woods-Saxon ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.3

1 次元 Woods-Saxon ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

4.4

3 次元 Woods-Saxon ポテンシャル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

クーロンポテンシャル

33

5.1

クーロンポテンシャルの特異性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

5.2

極座標表示の経路積分 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

5.3

物理量の計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

5.4

相関関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

5.5

l ≥ 2 に対するシミュレーション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

6

フーリエ展開による運動エネルギーの計算

41

7

まとめ

42

7.1

結論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

7.2

今後の展望・課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

1
1.1
2

3

4

5

2

1 イントロダクション
1.1 理論物理と数値計算
昨今の技術の発展によりコンピューターの計算速度は飛躍的に向上し、比較的安価なものでも従来に比べて格段に高性能なものとなった。それに伴い理論物理においてもコンピューターを用いて数値的（近似）解を求めるという手法はより一般的なものとなり、理論物理にとって強力なツールとして定着している。そこで本課
題研究では数値計算手法の理解と実践を目指すものとし、より具体的（標語的）には、『様々なポテンシャルの量子系で虚時間形式の経路積分を利用し、Metropolis アルゴリズムに基づいたモンテ
カルロ計算を実行する。』とする。
次章ではまず背景知識として虚時間形式の経路積分と Metropolis アルゴリズムに基づくモンテカルロ計算について説明をする。

3

2 経路積分法による量子力学へのアプローチとモンテカルロ計算
この章では虚時間形式の経路積分と Metropolis アルゴリズムについて述べる。

2.1 経路積分法の導出と虚時間形式
まず初めにこれからの議論における仮定を述べておく。

1. この導出においては表記を簡単にするため 1 次元 1 粒子系を扱う。（この節の終わりで 3 次元 1 粒子の場合についても結果と変更点を述べる。
）

2. 位置と運動量の固有状態の完全性の積分区間は ±∞ であるとする。
すなわち、

∫ ∞
−∞

∫ ∞

dq |q⟩ ⟨q| = 1,−∞

dp |p⟩ ⟨p| = 1.

(1)

3. ハミルトニアンは時間に陽に依存せず、以下の形をしているとする
Ĥ =

p̂2
+ V (q̂).
2m(2)

以上の 3 つを踏まえて、ある時刻 ti で状態 ψi にあった系が時刻 tf (> ti ) で状態 ψf に遷移する遷移振幅は

Heisenberg picture でH ⟨ψf | e−iĤ(tf −ti )/ℏ|ψi ⟩H(3)と書ける。これを Schrödinger picture で書き直すと、S ⟨ψf , −tf |ψi , −ti ⟩S .

(4)となる。（ブラやケットの添字 H, S はそれぞれ Heisenberg picture, Schrödinger picture であることをあらわに書いたものである。以降は必要に応じて適宜書いたり省略したりする。
）ここで、|q, t⟩H = eiĤt/ℏ |q⟩S

(5)について少し調べておくと、(∫∫
dq |q, t⟩ ⟨q, t| = eiĤt/ℏ)

dq |q⟩ ⟨q| e−iĤt/ℏ = 1(6)よりこの状態は完全で、また、H ⟨q, t|ψ, −t⟩S = ⟨q| e−iĤt/ℏ iĤt/ℏe|ψ⟩ = ⟨q|ψ⟩ = ψ(q)(7)

が成り立つことがわかる。この結果を利用すると、∫S ⟨ψf , −tf |ψi , −ti ⟩S ==∫dqf dqi ⟨ψf , −tf |qf , tf ⟩ ⟨qf , tf |qi , ti ⟩ ⟨qi , ti |ψi , −ti ⟩(8)dqf dqi ψf∗ (qf ) ⟨qf , tf |qi , ti ⟩ ψi (qi ).

(9)

上式で ψf∗ (qf ) と ψi (qi ) は時間に依存せず、系の時間発展の情報は遷移要素 ⟨qf , tf |qi , ti ⟩ に集約されておりこ
れを調べれば遷移振幅が計算できることになる。

4

そこで遷移要素H ⟨qf , tf |qi , ti ⟩H ≡ Z(qi , ti , qf , tf )(10)

= S ⟨qf | e−iĤ(tf −ti )/ℏ |qi ⟩S

(11)

2

p̂
を見ていく。仮定 2 の Ĥ = 2m
+ V (q̂) であることと、トロッター積の公式

)
(
)]N
[
(
]
t
t
Â expB̂
exp t(Â + B̂) = lim exp
N →∞
N
N
[

(12)

を利用すると

[ ( 2
)
]
p̂
Z = ⟨qf | exp −i+ V (q̂) (tf − ti )/ℏ |qi ⟩
2m
[
(
)
(
)]N
it p̂2it= lim ⟨qf | exp −exp −V (q̂)
|qi ⟩
N →∞
ℏN 2m
ℏN
[
(
)
(
)]N
iϵ p̂2iϵ= lim ⟨qf | exp −exp − V (q̂)|qi ⟩ .
N →∞
ℏ 2m
ℏ

(13)
(14)
(15)

ここで t = tf − ti , ϵ = t/N とした。
（後の便利のために tn = ti + nϵ としておく。また、t0 = ti , tN = tf .）さらに N − 1 個の完全性の条件∫
Z(q0 , t0 , qN , tN ) = lim

N →∞

∫

dqn |qn ⟩ ⟨qn | = 1(n = 1, · · · , N − 1) を間に挿入すると、
2
iϵ p̂

iϵ p̂

2

dq1 · · · dqN −1 ⟨qN | e− ℏ 2m e− ℏ V (q̂) |qN −1 ⟩ · · · ⟨q1 | e− ℏ 2m e− ℏ V (q̂) |q0 ⟩ . (16)iϵiϵ

なお上式において qf = qN , qi = q0 と表記を変えた。

(16) の積になっている被積分関数のひとつを取り出すと、
[
]
[
]
[
]
[
]
iϵ p̂2iϵiϵiϵ p̂2⟨qj+1 | exp −exp − V (q̂) |qj ⟩ = exp − V (qj ) ⟨qj+1 | exp −
|qj ⟩
ℏ 2m
ℏ
ℏ
ℏ 2m
∫
p2
iϵ
iϵ j
= e− ℏ V (qj ) dpj ⟨qj+1 |pj ⟩ e− ℏ 2m ⟨pj |qj ⟩ .

(17)
(18)

上の 2 番目の等号で運動量の固有状態の完全性（仮定 2）を用いた。
また、運動量の固有状態 |p⟩ について ⟨p|q⟩ = √ 1 e−ip·q/ℏ （平面波解）より (18) は、
2πℏ

e

− iϵ
ℏ V (qj )∫dpj − iϵ p2j ipj ·(qj+1 −qj )/ℏ
e ℏ 2m e
=
2πℏ

∫

[
]
dpj
iϵ
qj+1 − qjexp
pj ·− H(pj , qj ) .
2πℏ
ℏ
ϵ

(19)

あるいは上式において pj 積分を実行して− iϵℏ V (qj )e

ここで形式的に

∫

[
]
2
dpj
iϵ
m
imϵ qj+1 − qj 2
exp −pj − (qj+1 − qj ) +
(
)
2πℏ
2mℏ
ϵ
2ℏ
ϵ
(
√
(
)2 )
1
2πℏ
imϵ
q
−
q
j+1
j
− iϵ
V
(q
)
j
=e ℏ
·
·
exp
2πℏ
iϵ
2ℏ
ϵ
}]
[ { (
√
)2
m
iϵ m qj+1 − qj
− V (qj )
.
exp
=
2iπℏϵ
ℏ
2
ϵ
qj+1 − qj
≡ q̇j
ϵ→0
ϵ
lim

5

(20)
(21)

(22)

(23)

と定義すると ϵ → 0 の極限（N → ∞ の極限）において、(22) は

√

[ {
[
]
}] √ miϵ m 2iϵmexpq̇j − V (qj ) =exp
L(q̇j , qj )
2iπℏϵ
ℏ 2
2iπℏϵ
ℏ

(24)

と書くことができる。L(q̇, q) はラグランジアンである。
以上から、



N
−1
( m )N/2∑iϵZ(q0 , t0 , qN , tN ) = lim
dq1 · · · dqN −1
exp L(q̇j , qj )
N →∞
2iπℏϵ
ℏ j=0
[ ∫ tN]∫i= Dq exp
dtL(q̇, q)
ℏ t0
∫

(25)

(26)

m
ここで Dq = limN →∞ ( 2iπℏϵ)N/2 dq1 · · · dqN −1 .

これが Z(q0 , t0 , qN , tN ) =H ⟨qN , tN |q0 , t0 ⟩H を経路積分で書いた表式である。ただし q0 = qi , qN = qf という境界条件が課されていることに注意する。
以上の導出のプロセスからこの表式に登場する積分∫

Dq は境界条件 q0 = qi , qN = qf を満たす（区分的に滑らかな）全ての経路についての和を意味しており、これを一般化すると、この積分は粒子の取りうる全ての軌跡についての和と考えることができる。そして、それぞれの軌跡を古典作用の重みで足し合わせたものが遷移要素である。あえてこのことを強調した表式で書き表すと、

Z=

∑

eiS/ℏ

(27)

paths

とった具合になる。
次に (22) においてウィック回転と呼ばれる変換 t → −iτ (ϵ → −iϵτ ) を行うと、

√

[ {
[
{ (
}] √
}]
(
)2
)2
ϵτ
m qj+1 − qjϵτ m qj+1 − qj
m
m
exp
−
exp −
− V (qj )=+ V (qj )
2πℏϵτ
ℏ
2
ϵτ
2πℏϵτ
ℏ
2
ϵτ
(28)

より



)N/2 ∫
N−1∑m
ϵτZ(q0 , t0 , qN , tN ) = limLE dq1 · · · dqN −1 exp −
N →∞ 2πℏϵτ
ℏ j=0
[
]
∫
1
= Dq exp − SEℏ∫
SE = dτ LE(

LE (q̇, q) =

1
m
2

(

dq
dτ

(29)

(30)
(31)

)2
+ V (q)(32)と虚数の現れない形で書ける。この操作のことをユークリッド化と呼び、SE をユークリッド作用と呼ぶ。また、経路積分のこの表現を虚時間形式と呼ぶ。
最後に 3 次元の場合について考えてみる。最終的な表式を得るまでの導出の流れは上でやった 1 次
元の場合とまったく同じである。注意点としては (18) 式に相当する式の処理において、平面波解が、

6

1−ip·q/ℏ⟨p|q⟩ = (2πℏ)と変更される点である。結果を書いておくと、
3/2 e



N
−1
( m )3N/2∑iϵZ(q0 , t0 , qN , tN ) = lim
dq1 · · · dqN −1
exp L(q̇j , qj )
N →∞
2πℏϵ
ℏ j=0
[ ∫ tN]∫i= Dq exp
dtL(q̇, q)
ℏ t0
( )2
1
dq
L(q̇, q) = m+ V (q).
2
dt
∫

(33)

(34)
(35)

2.2 経路積分で基底状態における期待値を計算する方法
ここでは量子系における物理量として基底状態による期待値をとりあげ、これが経路積分形式の遷移要素から得られることを見る。
つ ま り 、|0⟩ を 基 底 状 態 と し た と き 任 意 の 物 理 量 Â の 期 待 値 ⟨0| Â |0⟩ を 経 路 積 分 形 式 を 利 用 し て（Schrödinger 方程式をあらわに解くことなく）計算するにはどうするかを考える。
遷移要素 ⟨qN , tN |q0 , t0 ⟩ = ⟨qN | e−iĤt/ℏ |q0 ⟩ で q0 = qN = q（周期境界条件）とおいて、虚時間 t → −iτ を導入すると、Z̄ = ⟨q| e− ℏ Ĥ |q⟩ =τ∑e−SE /ℏ .

(36)

paths

ここで β ≡ τℏ とし、いま求めた周期境界条件を満たす遷移要素を全 q について足しあげると、∫∫

dq ⟨q| e−β Ĥ |q⟩ = Tre−β Ĥ .

dq Z̄ =
上の 2 つめの等号では完全性∫

dq |q⟩ ⟨q| = 1 を用いた。

∑

続いて、エネルギー固有状態 |n⟩ もまた完全系をなしている (

∑

(
⟨n| e

−β Ĥ

|n⟩ =

n

よって、

∑

n |n⟩ ⟨n| = 1) とすると、

)
e

−βEn

= Tre−β Ĥ

(Ĥ |n⟩ = En |n⟩).

∫∑⟨n| e−β Ĥ |n⟩ =n∑e−βEn .

dq ⟨q| Âe−β Ĥ |q⟩ =∑⟨n| Âe−β Ĥ |n⟩ =

n

以上から、

∑

(39)

n

次に任意の演算子 Â に対して虚時間で ⟨q| Âe−β Ĥ |q⟩ を考えると、上と同じく

∫

(38)

n

dq Z̄ =から、

(37)

Â |n⟩ e−βEn
n ⟨n|
∑ −βE
n
ne

∑

∫

dq |q⟩ ⟨q| =⟨n| Â |n⟩ e−βEn .

∑n |n⟩ ⟨n| = 1

(40)

n

∫
=

dq ⟨q| Âe−β Ĥ |q⟩
∫
.
dq Z̄

(41)

上式の左辺について分母、分子両方に eβE0 （E0 は基底エネルギー）をかけると∑Â |n⟩ eβ(E0 −En ) β→∞n ⟨n|∑ β(E −E )−−−−→ ⟨0| Â |0⟩ .
0
n
ne
7

(42)

なお、E0 < En

(n = 1, 2, · · · ) を用いた。

よって∫dq ⟨q| Âe−β Ĥ |q⟩ β→∞∫−−−−→ ⟨0| Â |0⟩
dq Z̄

(43)

と左辺の量において β → ∞ の極限（β = τ /ℏ であり ℏ は定数なのでこの極限は τ → ∞ の極限を考えることに等しく以降では τ → ∞ と書くこともある）をとると Â の基底状態における期待値が得られることがわかる。

Â が q̂ のみの関数、すなわち Â = A(q̂) であるとすると上式の左辺は、]∫[∑A(q) paths e−SE /ℏ dq]∫ [∑−SE /ℏ dq
e
paths

(44)

と書くことができる。あるいはこれは実際的な用途を見据えての表式だが、

∫

dq

∑
⟨A(q(τ ))⟩ e−SE /ℏ∫ paths∑
dq paths e−SE /ℏ(45)と書くこともできる。これは重み e−SE /ℏ で分布する ⟨A(q(τ ))⟩ の期待値という格好をしている。ここで⟨A(q(τ ))⟩ はある経路（path）における虚時間平均を表す。
以上から、経路積分を用いて ⟨0| Â |0⟩ を計算する手順を（計算機を用いることを見据えて）構成すると、

1. まず、ある周期境界条件（q0 = qN = q ）を満たす path を生成する。
2. 上で生成された path の虚時間方向の A(τ ) の平均（これを ⟨A⟩ とする）をとる。
3. ⟨A⟩ を q0 = qN = q を満たすあらゆる path について求め、またあらゆる q についても同様の計算を行う。

4. 求まった各 path（あるいは配位と呼ぶこともある）における ⟨A⟩ に対し、e−SE /ℏ （重み）をかけて足しあげる。
となる。以上の step を毎度完全に random に（ただし周期境界条件は満たしつつ）配位を生成して実行してやれば（β → ∞ の極限で）原理的には ⟨0| Â |0⟩ が求まる。しかし、いくら計算機といえどこれでは大変な手間である。そこで Metropolis アルゴリズムに基づいたモンテカルロ計算なる手法を導入する。

2.3 Metropolis アルゴリズム
Metropolis アルゴリズムとは重み e−S に従う分布を生成することができるアルゴリズムである。
2.3.1 Metropolis アルゴリズムの仕組み
まず初めにこの Metropolis アルゴリズムで上記の重み e−S に従う分布が生成される仕組みを説明する。
まず、ある状態 r と s（例えば位置 r と s）についてそれぞれの状態への遷移が等価な状況（これを詳細釣り合い（detailed balance）と呼ぶ）を用意する。式で書くと

P (r → s) = P (s → r).
上式で P は確率を意味している。

8

(46)

このとき、ある時点でのそれぞれの状態の状態数を nr , ns とし、また各状態に作用（正確には各状態間の作用の差）が定義されていてその値がそれぞれ Sr , Ss (Sr > Ss ) であるとする。そして Metropolis ステップとして上の詳細釣り合いを満たしている環境に以下の条件を追加する。

(i) r から s（作用が大きい方から小さい方）に移動する数 ：nr P (r → s)(ii) s から r（作用が小さい方から大きい方）に移動する数 ：ns P (s → r) × e−(Sr −Ss )つまり作用が大きい状態から小さい状態への遷移（あるいは移動）は元の詳細釣り合いに従い、作用の小さい状態から大きい状態へは元の詳細釣り合いとそれに加えて e−(Sr −Ss ) で抑制された確率で遷移するという状況である。
以上の状況で平衡状態（見かけの移動なし）というものを考えると、

nr P (r → s) = ns P (s → r) × e−(Sr −Ss ) .

(47)

詳細釣り合い（P (r → s) = P (s → r)）から

nr eSr = ns eSs = const.

∴ ns ∝ e−Ss .

(48)

よって状態数は e−S に比例する。つまり、この Metropolis ステップを繰り返すことで平衡状態で重み e−S の分布を生成することができる。

2.3.2 Metropolis アルゴリズムの手順
前の小節で Metropolis アルゴリズムで重み e−S の分布を生成できることを見た。ここでは、より具体的に（わかりやすい形で）Metropolis アルゴリズムの手順をまとめておく。ここで重みは e−S[x] と設定する。

1. ある状態 xnew を random に生成する（ただし詳細釣り合いを満たすように生成しなければならないことに注意する）。

2. 作用の差 S[xnew ] − S[xold ] を評価する。これが負ならばこの xnew を採用（accept）しこれを改めて
xold と定義し始めのステップに戻る。
3. 前のステップで評価した作用の差が正ならば、0 < r < 1 を満たす r を random に生成し続いてe−(S[xnew ]−S[xold ]) − r を評価する。これが正ならばこの xnew を採用しこれを xold と定義し直して再び最初のステップに戻る。一方、負の場合には今回生成した xnew は採用せず（reject）、xold は変えずに最初のステップからやり直す。
このループを十分な回数繰り返すことで重み e−S[x] の分布を生成することができる。
一番最初のステップで詳細釣り合いを満たさなければならないことに注意すると書いたが本課題研究では

xnew = xold + hstep · (2r − 1) (0 ≤ r ≤ 1, r : random, hstep : f ixed)とすることでこれを満足するようにした。

9

(49)

3 調和振動子ポテンシャル
この章から先では具体的なポテンシャルの下で実際に数値計算を行い、結果を考察していく。この章では解析的に解くことができ、結果のよく知られている調和振動子ポテンシャルのもとで基底状態の粒子存在確率および、エネルギー準位を計算し、Schrödinger 方程式の解と比較する。

3.1 1 次元調和振動子ポテンシャル
以下では特段の記述のない限り、Nτ = 100、Nsweep = 10000、hstep = 1.0 のもと、数値計算を行った。

1 次元調和振動子のポテンシャルは角振動数を ω とすると
V (x) =

1mω 2 x22(50)と書けて、ここでは各パラメーターを下記のように設定し、シミュレーションを行った。

ℏ = 1,

m = 1,

ω=1(51)a) 熱平衡化（thermalization）までの sweep 数
まず、ここで初めて出てきた熱平衡化という言葉について説明を加えておく。前章の理論は無限回の sweepを考えていたので初期条件（初期配位）をどのようにとっても正しい値を得ることができるというものだったのだが、実際の計算機による計算では有限回の sweep になるため、計算結果は初期値依存性を少なからず持ってしまう。この初期値依存性は本来の物理を反映したものではないので、できるだけ抑えたいもので、ここで熱平衡化というプロセスを利用する。どのような初期条件に設定しても配位生成を繰り返していくことで、すべて似た配位に収束していくということが経験的に知られていてこのプロセスを熱平衡化と呼ぶ。

実際の計算結果（図 1）を見てみると、（この初期条件の範囲内では）どの初期条件から出発しても多めに見積もっても 50 sweep までには各配位における x の期待値（虚時間平均）は同じような値に収束しているのが見てとれる。
そこで初めの 50 sweep までのデータを無視してやれば初期値の寄与を大きく減らせると考え、今後は、50

sweep 以降のデータのみを対象として各物理量を計算することとする。

10

図 1 1 次元調和振動子における sweep 数と各 sweep での位置 x の虚時間平均の関係性 : 横軸は sweep
数、縦軸は各 sweep における位置 x の虚時間平均を表す。初期配位として各虚時間 τi における xτi を

xτi = 1, 3, 5, 10, −5 として、そこから配位生成をスタートしたものをそれぞれ色分けして描いてある。初期条件を変えたときに、各配位における x の期待値が配位生成とともにどのように変化するか確認できる。

b) ヒストグラムと波動関数の対応
経路積分で [x, x + ∆x] に状態 (path) を見出す確率は、

1
P (x; β) =
β

∫ β
dτ

′

0

∫ x+∆xdx′

xZ(xF , β; x′ , τ ′ )Z(x′ , τ ′ ; xI , 0)Z(xF , β; xI , 0)(52)

∆x が十分小さいとき、∫ x+∆xdx′ = ∆xδ(x − x′ )

(53)

x

であり、遷移確率の一般式は、

Z(xF , τF ; xI , τI ) =∞∑e−En (τF −τI ) ϕn (xF )ϕ∗n (xI )

(54)

0

このとき、τ 積分は実行可能になり、計算を進めると、P (x; β)= |ϕ0 (x)|2 + O((β∆E)−1 )
∆x

(55)

確率密度は規格化されたヒストグラムとみなすこともでき、
tot
∑
P (x)
1
θ(∆x − |xi − x|)=∆x∆xNtot

N

{i}

11

(56)

よって確率密度は、
tot∑
1|ϕ0 (x)| =θ(∆x − |xi − x|)∆xNtot

N

2

(57)

{i}

と表すことができる。
一次元調和振動子の基底波動関数は、

1ϕ = N0 exp(− α2 x2 )
2√
mω
α=
ℏ
( mω )1/4
N0 =
πℏ

(58)
(59)
(60)

とあらわせる。各位置での存在確率は、
2|ϕ| = N02 exp(−α2 x2 ) = N02 exp(−

mω 2
x )
ℏ

(61)

ここで、現在の m, ω, ℏ の値を代入すると、|ϕ| = π −1/2 exp(−x2 )2であり、ガウシアンの形になることがわかる。
図 2 には、数値計算したときの確率密度を示した。

図 2 1 次元調和振動子における確率密度とフィッティングしたガウシアン：横軸は位置 x、縦軸は確率密度を表す。また、f(x) は式（62）に基づく確率密度、g(x) はフィッティング関数を表す。 g(x) =2

a exp(− (x−b))。a = 0.585(±0.000), b = 0.004(±0.000), c = 0.929(±0.001).
c

式（62）とガウシアンフィッティングの値が、対応しているのがわかる。

12

(62)

c) エネルギー固有値
Schrödinger 方程式から、調和振動子のエネルギー固有値は、
En =1(n + 1)ℏω (n = 0, 1, 2...).
2(63)

今回 ℏ = 1, ω = 1 としているので、基底状態のエネルギー固有値は、

E0 =

1
1
ℏω =
2
2

(64)

となることが期待される。
シミュレーションの結果、E0 = 0.464(±0.077) であった。ただし、運動エネルギー K の計算は、ビリアル
の定理

⟨K⟩ =

1
⟨xV ′ ⟩
2

(65)を用いて評価した。
次に、E0 の ℏ 依存性を探る。図 3 のように、E0 と ℏ が比例関係を示しているのがわかる。比例係数は、

a = 0.448(7) となっており、想定値の 1/2 におおむね近い値になっている。

図 3 1 次元調和振動子ポテンシャルにおける基底エネルギーと ℏ の関係性：横軸はプランク定数 ℏ、縦軸
は基底エネルギー E0 。点は ℏ を 0.1 と 0.5 から 10 まで 0.5 きざみでの値で計算したものである。直線はフィッティング関数で f (ℏ) = ax + b である。a = 0.448(7), b = -0.00068(24)

第一励起状態のエネルギーについても計算することができる。基底状態に対しての二点間の相関は次のよう
になる。

13

Γ(τ ) = ⟨0|x(0)x(τ )|0⟩ − ⟨0|x|0⟩2
2= ⟨0|x(0) exp(Ĥτ ℏ)x(0) exp(−Ĥτ /ℏ)|0⟩ − ⟨0|x|0⟩∑2=⟨0|x(0) exp(−Ĥτ /ℏ)|n⟩ ⟨n|x(0) exp(Ĥτ /ℏ)|0⟩ − ⟨0|x|0⟩n=∑exp(−(En − E0 )τ /ℏ) | ⟨0|x|n⟩ |2 − ⟨0|x|0⟩

2

n

=

∑

exp(−(En − E0 )τ /ℏ) | ⟨0|x|n⟩ |2n(̸=0)→ exp(−(E1 − E0 )τ /ℏ) | ⟨0|x|1⟩ |2(τ → ∞)

よって少しずれた時間間隔の相関関数の比をみることで第一励起状態のエネルギーについての情報を得ることができる。

Γ(τ + ∆τ )= exp(−(E1 − E0 )∆τ /ℏ)
T →∞
Γ(τ )
(
)
Γ(τ + ∆τ )ℏ
E1 = E0 −
log
∆τ
Γ(τ )
lim

(66)
(67)

図 4 が一次元調和振動子における位置の相関関数である。τ が大きくなると相関が小さくなっている
ことが分かる。周期境界条件を用いているため τ = 5 のあたりで最も小さくなり、また大きくなっている。y = a exp(bx) で fitting した時の b の大きさがだいたい 1 に近い値になっているためエネルギー差ℏω = 1 (ℏ = 1, ω = 1) を正しく出せたと言える。式 (67) のように、得られた時間間隔の相関関数と少しずれた時間間隔 (∆τ = dτ ) の相関関数の比の log をとると図 5 のようになる。τ = 0 周辺もしくは τ = 10 周辺で正しい値がでている。

図 4 1 次元調和振動子ポテンシャルにおける相関関数：横軸は τ 、縦軸は Γ(τ ) で log スケールになっている。y = a exp(bx) で fitting した。a = 0.485(2), b = -1.017(6), WSSR/ndf(reduced chisquare) =0.2538514図 5 1 次元調和振動子ポテンシャルにおける E1 − E0 ： 横軸は τ 、縦軸は log

Γ(τ +∆τ ).
Γ(τ )d) 不確定性関係
調和振動子における n 番目の励起状態の不確定性関係は

1∆x∆p = ℏ(n + )2(68)である。したがって基底状態の不確定性は 1/2ℏ となる。色々な ℏ における不確定性を計算すると図 6 となった。おおむね傾き 1/2 の直線にのっており、正しく基底状態にアクセスできているといえる。

図6

1 次元調和振動子ポテンシャルにおける基底状態の不確定性と ℏ の関係：横軸はプランク定数 ℏ、

縦軸は不確定性 ∆x∆p。ℏ を 0.1 と 0.5 から 10 まで 0.5 きざみの値で計算したものである。直線はフィッ
ティング関数で y = ax + b である。a = 0.4911(7),b = -0.00063(24)153.2 3 次元調和振動子ポテンシャル
以下では特段の記述のない限り、Nτ = 100、Nsweep = 10000、hstep = 1.0 のもと、数値計算を行った。
調和振動子のポテンシャルは動径距離を r として

V (r) =

1mω 2 r22(69)と書けて、ここでは各パラメーターを下記のように設定し、シミュレーションを行った。

ℏ = 1,

m = 1,

16

ω=1

(70)

a) 熱平衡化までの sweep 数
図 7 より、（この初期条件の範囲では）50 回目の sweep には熱平衡に達しているとみなし、今後は 50 sweep
以降のデータのみを対象として各物理量を計算することとした。

図 7 3 次元調和振動子における sweep 数と各 sweep での動径距離 r の虚時間平均の関係性: 横軸はsweep 数、縦軸は各 sweep における動径距離 r の虚時間平均を表す。初期条件として各虚時間 τi における xτi , yτi , zτi を 1, 3, 5, -5 にそろえてから配位生成をスタートしたものをそれぞれ色分けして描いてある。初期条件を変えたときに、各配位における動径距離 r の期待値が配位生成とともにどのように変化す
るか確認できる。

17b) ヒストグラムと波動関数の対応
まず、3 次元空間での各格子点付近の count 数を図 8 に示した。
図 8 から原点に向かうにつれ count 数が増加することが定性的にわかる。

図 8 3 次元調和振動子ポテンシャルにおける 3 次元ヒストグラム：図の x, y, z 軸は実空間の x, y, z（座標）に対応する。count 数は点の色の濃淡で表現してある。この図から原点に近づくほど、count 数が増えていることが定性的にわかる。

18

次に、r（動径）方向に等間隔に区切ったヒストグラムを図 9 に示した。

図 9 3 次元調和振動子ポテンシャルにおける r 方向ヒストグラム：横軸は動径距離、縦軸は count 数を表す。動径方向に等間隔に区切った各球殻上で count したものが描画してある。

図 9 は図 8 の原点に向かうにつれ count 数が増加するという結果と一見相反しているように見えるが、こ
れは、

∫

∫
dx3 ρ = 4π

drr2 ρ(r)(71)

における r2 ρ(r) に対応したもので ρ(r) ではないことに注意する。（動径方向に等間隔に区切ってヒストグラ
ムを作成していた。）
そこで、各 count 数を r2 で割った、ρ(r) に対応するヒストグラムを改めて図 10 に示す。
図 10 からは原点に向かうにつれ count 数が増加するのが正確に見てとれる。

19

図 10 3 次元調和振動子ポテンシャルにおける動径方向の確率密度 ρ(r) に対応する r 方向ヒストグラム:
図 9 の縦軸を r 2 で割ったもの。ガウシアンのような形になっていることがわかる。

c) エネルギー固有値
Schrödinger 方程式から、3 次元調和振動子基底状態のエネルギー固有値は、
E0 =

3
ℏω.
2

(72)

今回 ℏ = 1, ω = 1 としている。
シミュレーションの結果、E0 = 1.510(±0.137) であり 3/2 は、誤差範囲の中にあり想定どおりの値がでていることがわかる。
次に、E0 の ℏ 依存性をを探る。図 11 のように、E0 と ℏ が比例関係を示しているのがわかる。比例係数
は、a = 1.293 (± 0.006) となっており、想定値の 3/2 におおむね近い値になっている。

20

図 11 3 次元調和振動子ポテンシャルにおける基底エネルギーと ℏ の関係性：横軸はプランク定数 ℏ、縦軸は基底エネルギー E0 を表す。フィッティング関数は y = ax + b で ω = 1, 2, 3, 4, 5 のデータでフィッ
ティングを行った。a = 1.293 (± 0.006) 、b = 0.215 (± 0.012) , (reduced chisquare) = WSSR/ndf :0.00369936

3.3 調和振動子のまとめ調 和 振 動 子 ポ テ ン シ ャ ル に つ い て 、Metropolis ア ル ゴ リ ズ ム に 基 づ く モ ン テ カ ル ロ 計 算 に よ っ て 、Schrödinger 方程式による解析解と近い値を算出することができた。このことは、今回用いた手法の信頼性を支持する。
一方で誤差を含めても解析解と僅かに一致しない部分もあり、これは計算の精度を上げる（虚時間の分割数を増やす、虚時間方向をより遠くまでとる、sweep 数をより多くとる、など）ことで改善することが期待されるがこの事実は解析解が分かっていない問題にこの手法を適用する際には慎重にパラメータを選ばなければならないことを意味する。実際には徐々に精度を上げていきどのような値に漸近しそうかをチェックすることになるであろう。また、計算機そのものが持つ性質による問題（丸め誤差の問題や疑似乱数の 乱数らしさ の問題など）があることも常に忘れてはならない。

21

4 Woods-Saxon ポテンシャル
この章では、Woods-Saxon ポテンシャルと呼ばれる原子核物理に登場するポテンシャルに束縛された 1 粒子系を扱う。
まず初めに Woods-Saxon ポテンシャルについて原子核物理の観点から解説し、これを取り扱う意義と今回の手法が利用できることを確認したのち、実際の計算結果に移る。

4.1 殻模型と Woods-Saxon ポテンシャル
原子核内における核子の性質を記述する方法として、平均場を用い、ある核子（1 粒子）とそれ以外の核子
によるポテンシャルという形で記述する殻模型（shell model）が知られている。
このモデルによって本来多くの核子による多体問題である原子核の問題を 1 核子とポテンシャルという 1 体
問題に帰着させることができ、その点において非常に有意である。
原子核の動径方向の電荷密度分布は、2 つのパラメーターをもつ Fermi 関数、

ρ(r) =

ρ(0)1 + e(r−c)/a(73)でよく近似できることが分かっており、この密度分布を再現するポテンシャルとして Woods-Saxon ポテン
シャル

Vcenter (r) = −V0
1 + e(r−R)/a(74)が知られている（図 12 参照）。
そして、この中心力ポテンシャルに LS 項を付け加えた

V (r) = Vcenter (r) + Vls (r)⟨l · s⟩ℏ2(75)

によって、原子核内の核子のエネルギー準位（特に縮退度）をよく説明することができる。
この Wood-Saxon ポテンシャル（と LS 項）によって決定されるエネルギー準位に各々の核子（通常陽子と中性子）が Pauli の排他律に従って、下から詰まっていくという描像である。

4.2 ハイパー核と Woods-Saxon ポテンシャル
ハイパー核とは核子として陽子と中性子の他にハイペロン（s クォークを含むバリオン）が存在している原子核のことで、実験的にはこのハイペロンが陽子や中性子と 区別できることから原子核内の個々の核子のエネルギー準位を調べる探子（probe）として利用することもあり、原子核の構造を知るにあたってこのハイ
ペロン（あるいはハイパー核）は重要な存在である。
ここで、ハイペロンを 1 つだけ含む原子核を殻模型の描像（Woods-Saxon ポテンシャル）で考えてみる。
ハイペロンは陽子や中性子と異なる粒子なので、ハイペロンが 1 個だけの状況ではこのハイペロンが Pauliの排他律を感じる ことはなく、この場合ハイペロンは最も低いエネルギー準位（基底状態）に落ち着く。

また、基底状態においては軌道角運動量 l = 0 なので LS 力が働くことはなくその意味でこのハイペロンの感じるポテンシャルはV (r) = Vcenter (r) = −22

V0
1 + e(r−R)/a(76)

ということになる。
以上の話のうち、• ハイペロンは基底状態に配置される。
• ハイペロンにはその基底状態においては LS 力は働かない。
の二点により今回の手法でこの物理を扱えることがわかる。
そしてこの章ではこのハイペロンを 1 つ含む原子核という系を背景に Woods-Saxon ポテンシャルのもとでの 1 粒子系の基底状態における振舞いを Metropolis アルゴリズムに基づいたモンテカルロ計算で見ていく。
これから実際の計算結果に移るが、以下では特段の記述のない限り、Nτ = 100、Nsweep = 10000、hstep = 1.0のもと、数値計算を行った。
改めて Woods-Saxon ポテンシャルの表式を書くと、V (r) = −V0
1 + e(r−R)/a(77)とこれは 3 つのパラメーター R, a, V0 を持つ。この 3 つのパラメーターについて今回は以下のようにとった。

R = 4.0,

a = 0.5,

V0 = 4.0(78)

注意として本来これらのパラメーターは実験値等から決定するものであるのだが、ここでは 勝手に このよ
うに設定した。
このときポテンシャルは図 12 のような形になっている。

図 12 Woods-Saxon ポテンシャル：横軸は原点を通る任意の直線 x の座標、縦軸はポテンシャルの数値を表す。式（77）において R = 4.0, a = 0.5, V0 = 4.0 としたもの。

さらにその他のパラメーターは下記のように設定し、シミュレーションを行った。

ℏ = 1,

m=1
23

(79)

4.3 1 次元 Woods-Saxon ポテンシャル
本来 Woods-Saxon ポテンシャルは 3 次元空間上で考えるものであるが、簡単のために 1 次元で行ってみた結果を載せておく。

a) 熱平衡化までの sweep 数
図 13 を見ると、初期条件のとり方が ±5 までの範囲のものについては、200 回目の sweep までの間に熱平衡に至っているとみなしても良さそうである。そこで、今後の計算では 200 sweep 以降のデータのみを対象として各物理量を計算することとした。

図 13 1 次元 Woods-Saxon ポテンシャルにおける sweep 数と各 sweep での位置 x の虚時間平均の関係
性：横軸は sweep 数、縦軸は各 sweep における位置 x の虚時間平均を表す。初期配位として各虚時間 τiにおける xτi を xτi = 1, 3, 5, 10, 15, -5 として、そこから配位生成をスタートしたものをそれぞれ色分けして描いてある。初期条件を変えたときに、各配位における x の期待値が配位生成とともにどのように変化するか確認できる。

24b) ヒストグラムと波動関数の対応
図 14 にヒストグラムとフィッティングしたガウシアンを示した。
ここでこのヒストグラムをガウシアンでフィッティングした理由は、量子力学の一般論、すなわち束縛状態の基底波動関数はゼロ点（node）を持たないこと、エネルギー固有値がポテンシャルより小さい領域では波動関数は exp で減少すること、そして波動関数の規格化可能性と（規格化された）ヒストグラムが波動関数の振幅の二乗であることを考慮にいれ、ガウシアンでフィッティングできるのではないかと考えたからである。
実際ヒストグラムの形からもガウシアンが想起される。

図 14 1 次元 Woods-Saxon ポテンシャルにおけるヒストグラムとヒストグラムにフィッティングしたガ
ウシアン：横軸は位置 x、縦軸は count 数を表す。f (x) = a exp(−(x−b)2) にてフィッティングしている。
2c2

a = 107.3(±1.2), b = 0.014(±0.002), c = 2.945(±0.027), reducedchisquare = 19.9569

25

図 15 にヒストグラムから算出した確率密度とフィッティングしたガウシアンを示した。

図 15 1 次元 Woods-Saxon ポテンシャルにおける確率密度とフィッティングしたガウシアン： 横
軸は位置 x、縦軸は確率密度を表す。f (x) = a exp(−(x−b)2) にて フィッティングしている。ac=0.296(±0.000), b = 0.011(±0.002), c = 3.75(±0.01), reducedchisquare = 3.6848e − 005

これを見ると、ガウシアンでフィッティングできるのではないかという期待に沿った結果が得られている。
そこで基底状態の確率密度がガウシアンになることが分かっている調和振動子ポテンシャルを参考に確率密度の表式（式（61））を図 15 のフィッティング関数に対応させて換算の ω を計算すると、

ω = 1/c = 0.266(80)と求まる。ただし上式における c は図 15 のフィッティングパラメーター c である。
さらにこの換算 ω を角振動数とする調和振動子ポテンシャルと今回採用した Woods-Saxon ポテンシャルとを比較したものを図 16 に示した。ただし、図 16 では形の比較がしやすいように原点における調和振動子のポテンシャルの値を Woods-Saxon ポテンシャルの値に一致させて（y 軸に沿って平行移動させて）ある。図16 から、この換算 ω による調和振動子ポテンシャルは Woods-Saxon ポテンシャルの底の形を反映したものになっていることがわかる。
また、図 16 のポテンシャルのもとで両者の基底エネルギー EW （Woods-Saxon）, Eo （調和振動子）をそれぞれ計算すると、

EW = −3.961 ± 0.056,

Eo = −3.857 ± 0.043

(81)

と求まる。

EW の誤差内の最大値 EW max = −3.905 と Eo の誤差内の最小値 Eomin = −3.900 とを比較するとわかる
ように、誤差幅込みでも一致はしていない。

26

図 16 Woods-Saxon ポテンシャル（図 12）と角振動数 ω の値が式（80）の調和振動子ポテンシャル：横軸は位置 x、縦軸はポテンシャルの数値を表す。形の比較がしやすいように原点における調和振動子のポ
テンシャルの値を Woods-Saxon ポテンシャルの値に一致させて（y 軸に沿って平行移動させて）表示してある。

c) エネルギー固有値
基底エネルギー E0 の ℏ 依存性を図 17 に示した。
この図からは一見比例関係にあるように見えるが、傾きが a=0.080 であることと y 切片として-4 よりわず
かに小さな値が出ていることから必ずしも比例関係とはこの結果だけでは結論付けられない。よって ℏ 依存性をより正確に確認するのであればさらに多くのデータが必要となる。

27

図 17 1 次元 Woods-Saxon ポテンシャルにおける基底エネルギーと ℏ の関係性：横軸はプランク定数ℏ、縦軸は基底エネルギー E0 を表す。フィッティング関数は y = ax + b で ω = 1, 2, 3, 4, 5 のデー
タでフィッティングを行った。a = 0.080 (± 0.008) 、b = -4.04 (± 0.06) , (reduced chisquare) =WSSR/ndf : 0.0365799

4.4 3 次元 Woods-Saxon ポテンシャルa) 熱平衡化までの sweep 数
図 18 を見ると、初期条件のとり方が ±5 までの範囲のものについては、200 回目の sweep までの間に熱平衡に至っているとみなしても良さそうである。そこで、今後の計算では 200 sweep 以降のデータのみを対象として各物理量を計算することとした。

28

図 18 3 次元 Woods-Saxon ポテンシャルにおける sweep 数と各 sweep での動径距離の虚時間平均の関係性：横軸は sweep 数、縦軸は各 sweep における動径距離 r の虚時間平均を表す。初期条件として各虚時
間 τi における xτi , yτi , zτi をそれぞれ 1, 3, 5, -5 にそろえてから配位生成をスタートしたものをそれぞれ色分けして描いてある。初期条件を変えたときに、各配位における動径距離 r の期待値が配位生成とともにどのように変化するか確認できる。

b) ヒストグラムと波動関数の対応
まず、3 次元空間での各格子点付近の count 数を図 19 に示した。
図 19 から原点に向かうにつれ count 数が増加することが定性的にわかる。

29

図 19 3 次元 Woods-Saxon ポテンシャルにおける 3 次元ヒストグラム：図の x, y, z 軸は実空間の x, y,z（座標）に対応する。count 数は点の色の濃淡で表現してある。この図から原点に近づくほど、count 数が増えていることが定性的にわかる。

次に、r（動径）方向に等間隔に区切ったヒストグラムを図 20 に示した。
調和振動子ポテンシャルの時と同様、これは、∫∫
3dx ρ = 4π

drr2 ρ(r)(82)

における r2 ρ(r) に対応したものである。
そこで、各 count 数を r2 で割った、ρ(r) に対応するヒストグラムを改めて図 21 に示す。
図 21 において大きな r については単調減少になっているが r の小さい領域では調和振動子ポテンシャルの時ほどきれいな形になっていないのは必ずしもこれが物理を反映したものというわけではなく、単に小さな rで割ると図 20 のヒストグラムの僅かな差が強調されてしまうことに起因するものなので寧ろ小さな r における図 21 のヒストグラムはあまり信用できないものといえる。
上の事情をふまえて図 21 を見ると、原点に向かうにつれ count 数が増加する傾向は概ね見えているとい

30

える。

図 20 3 次元 Woods-Saxon ポテンシャルにおける r 方向ヒストグラム：横軸は動径距離、縦軸は count数を表す。動径方向に等間隔に区切った各球殻上で count したものが描画してある。

図 21 3 次元 Woods-Saxon ポテンシャルにおける動径方向の確率密度 ρ(r) に対応する r 方向ヒストグ
ラム：図 20 の縦軸を r 2 で割ったもの。ガウシアンのような形になっていることがわかる。

31c) エネルギー固有値
基底エネルギー E0 の ℏ 依存性を図 22 に示した。
この図は 1 次元 Woods-Saxon ポテンシャルの結果（図 17）よりは比例関係になっているように見えるが、やはり y 切片として-4 より小さな値が出ているので ℏ 依存性を判断するにはデータが足りないと言わざるを得ない結果と言える。

図 22 3 次元 Woods-Saxon ポテンシャルにおける基底エネルギーと ℏ の関係性：横軸はプランク定数ℏ、縦軸は基底エネルギー E0 を表す。フィッティング関数は y = ax + b で ω = 1, 2, 3, 4, 5 のデー
タでフィッティングを行った。a = 0.692 (± 0.009) 、b = -4.21 (± 0.02) , (reduced chisquare) =WSSR/ndf : 0.0192774

32

5 クーロンポテンシャル
次のような 3 次元のクーロンポテンシャルに対して諸物理量を求めてみる。

V (r) = −αr(r = |r|)

5.1 クーロンポテンシャルの特異性
以下では次のパラメータでシミュレーションを行った。

ℏ = 1,

m = 1,

α=1(83)

様々な初期配位に対して、sweep を繰り返していくと平均半径 r がどのように変わっていくか見るために、
各 sweep 数における r を計算すると図 23 のようになった。

図 23

ク ー ロ ン ポ テ ン シ ャ ル に お け る sweep 数 と 各 sweep で の 平 均 半 径 x : 初 期 配 位 を rτi =0.001, 0.005, 0.015, 0.1, 1.3 として、sweep を重ねていった。

初期配位を原点に近づけていくと平均半径はその分小さくなり、どの初期配位に対しても共通の一定の値への収束が見られない。これはクーロンポテンシャルが中心 r = 0 での特異性を持っているからであり、中心に近づくほどエネルギーが下がり安定になる、古典的描像を見せる。この手順では水素原子の基底状態を達成することができない、もしくは基底状態の実現には sweep 数が到底足りていないと言える。このように下に有界でないポテンシャルに対して今用いている Metropolis 法は有効でない。

5.2 極座標表示の経路積分
原点での特異性を除くために遠心力ポテンシャルをラグランジアンに加える。この項が現れるのは運動方程式を極座標表示で表した場合であり、そのために今まではデカルト座標で表示していた経路積分を極座標で考33

えてみる。
デカルト座標のときと同様に微小遷移振幅の計算を行う。極座標表示で中心力ポテンシャル V (r) がある系のハミルトニアンは

pˆr 2
L̂2
+
+ V (r̂)2m 2mr̂2
ℏ1 ∂
r
pˆr ≡
i r ∂rĤ =となる。演算子 pˆr の固有状態 |pr ⟩ に対する固有値を pr とする。

pˆr |pr ⟩ = pr |pr ⟩演算子 pˆr は次の交換関係を満たすため、確かに動径方向の運動量であるといえる。

[r̂, pˆr ] = r

ℏ1 ∂
ℏ1 ∂
r−
r r = iℏi r ∂ri r ∂r

極座標表示における位置の動径方向と運動量の動径方向のブラケットは次の微分方程式を満たす。

pr ⟨r|pr ⟩ = ⟨r|pˆr |r⟩ =ℏ⟨r|pr ⟩iよって次のような球面波で表すことができる。

⟨r|pr ⟩ =1i
exp( rpr )rℏ

遷移振幅に対してリートロッター分解を適用すると[]
1
ĤT |xi ⟩
iℏ
{
[
]
[
]
[
]}N
∆t pˆr 2∆t∆t L̂2exp
exp
V (r̂)|xi ⟩= lim ⟨xf | exp
N →∞iℏ 2mr̂2
iℏ 2m
iℏ
]
[
[
]
[
]
∫
∆t pˆr 2∆t∆t L̂2
exp
exp
V
(r̂)
|xN −1 ⟩ ⟨xN −1 | · · · |x1 ⟩= lim
dxN −1 · · · dx1 ⟨xf | exp
N →∞iℏ 2mr̂2
iℏ 2m
iℏ
[
]
[
]
[
]
∆t L̂2∆t pˆr 2∆t· ⟨x1 | exp
expexp
V (r̂) |xi ⟩iℏ 2mr̂2iℏ 2miℏ⟨xf , tf |xi , ti ⟩ = ⟨xf | exp

微小遷移振幅を抜き出して計算する。

]
[
]
[
]
∆t pˆr 2∆t∆t L̂2exp⟨xn+1 | exp
exp
V (r̂) |xn ⟩iℏ 2mr̂2
iℏ 2m
iℏ
]
[
[
]
[
]
∫
∆t pˆr 2∆t∆t L̂2
|x⟩
⟨x|
exp
exp
V
(r̂)
|xn ⟩
= dx ⟨xn+1 | expiℏ 2mr̂2
iℏ 2m
iℏ
[

B

A

34

(84)

]
∆t L̂2(A) =⟨xn+1 | exp|l, m⟩ ⟨l, m|x⟩iℏ 2mr̂2l,m[]∑∆t L̂2=⟨xn+1 | exp|l, m⟩ ⟨l, m|x⟩iℏ 2mr̂2
l,m][∑∆t ℏ2 l(l + 1)
|l, m⟩ ⟨l, m|x⟩=⟨xn+1 | exp
2iℏ 2mrn+1
l,m[]∑∆t ℏ2 l(l + 1)=
exp
⟨rn+1 |r⟩ ⟨θn+1 , ϕn+1 |l, m⟩ ⟨l, m|θ, ϕ⟩2iℏ 2mrn+1
l,m[]∑∆t ℏ2 l(l + 1) 1∗=
exp
δ(rn+1 − r)Ylm (θn+1 , ϕn+1 )Ylm
(θ, ϕ)
2iℏ 2mrn+1
r2
∑

[

l,m

[
]
[
]
∆t pˆr 2
∆t
dp
⟨x|
exp
|p⟩
⟨p|
exp
V
(r̂)
|xn ⟩
(2πℏ)3
iℏ 2m
iℏ
[
]
[
]
∫
dpr p2r∆t pr 2
∆t
=
dΩ
⟨r,
θ,
ϕ|
exp
|p
,
θ
,
ϕ
⟩
⟨p
,
θ
,
ϕ
|
exp
V
(r
)
|r, θn , ϕn ⟩
p
r p
p
r p
p
n
(2πℏ)3
iℏ 2m
iℏ
[
]
∫
dpr p2r∆t pr 2
∆t
=
dΩ
exp
+
V
(r
)
⟨r, θ, ϕ|pr , θp , ϕp ⟩ ⟨pr , θp , ϕp |rn , θn , ϕn ⟩pn(2πℏ)3
iℏ 2m
iℏ
[
]
∫
∆t pr 2∆t
dpr p2r
exp
+
V
(r
)
⟨r|pr ⟩ ⟨pr |rn ⟩ ⟨θ, ϕ|θn , ϕn ⟩=n(2πℏ)3
iℏ 2m
iℏ
[
]
∫
dpr p2r∆t pr 2
∆t
1
i
1
i
1
=
exp
+
V
(r
)
exp( pr r) exp(− pr rn )δ(θ − θn )δ(ϕ − ϕn )n(2πℏ)3
iℏ 2m
iℏ
r
ℏ
rn
ℏ
sin θ
[
]
[
]
∫ ∞
∆t
1
∆t pr 2
1
dpr 2
i
=
exp
V (rn )
δ(θ − θn )δ(ϕ − ϕn )
p
exp
+
p
(r
−
r
)
r
n
rrn
iℏ
sin θ
(2πℏ)3 riℏ 2m ℏ
0
}]
[
[
{
]
∫ ∞
1
1
1
dpr 2
∆t
r − rn 2i∆t m r − rn 2=(
) − V (rn )δ(θ − θn )δ(ϕ − ϕn )
p
exp
−i
(p
−
m
)
exp
r
rrn
ℏ
2
∆t
sin θ
2 −∞ (2πℏ)3 r
2mℏ
∆t
[
{
}]
3
(
)
1
m
1
i∆t m r − rn 2
1
2
=
exp
(
) − V (rn )δ(θ − θn )δ(ϕ − ϕn )4π 2πi∆tℏ
rrn
ℏ
2
∆t
sin θ
∫

(B) =

]
∆t ℏ2 l(l + 1) 1∗
(84) = drr dΩδ(rn+1 − r)Ylm (θn+1 , ϕn+1 )Ylm
exp
(θ, ϕ)
2iℏ 2mrn+1
r2
l,m
[
{
}]
i∆t m r − rn 2
1
1 ( m ) 32 1
exp(
) − V (rn )δ(θ − θn )δ(ϕ − ϕn )·4π 2πi∆tℏ
rrn
ℏ
2
∆t
sin θ
}]
[
{
∑ ( m ) 32
1i∆t m rn+1 − rn 2ℏ2 l(l + 1)∗=
Ylm (θn+1 , ϕn+1 )Ylm
(θn , ϕn )exp(
) − V (rn ) −22πi∆tℏ
4πrn+1 rnℏ
2∆t
2mrn+1
∫

2

∑

[

l,m

求めた微小遷移振幅を合成して有限区間の遷移振幅を計算する。遷移振幅の dΩ 積分を実行すると球面調和関数の性質より、微小遷移振幅中の l, m は等しくなる。遷移振幅中の積分は dr だけになり

35

∫ N∏
−1 {

(

1

m ) 32 2
rj drj2πi∆tℏ}

1 ( m ) 32
N →∞4πrj+1 rj
4πr1 ri 2πi∆tℏ
j=1

{
}
N
−1
2
∑
∑
i
m
r
−
r
ℏ
l(l
+
1)
j+1
j 2
∗
 Ylm (θf , ϕf )Ylm
·
exp 
∆t
(
) − V (rj ) −(θi , ϕi )
2
ℏ j=0
2
∆t
2mrj+1
l,m
[ ∫
]
∫
∑
1
i
∗
=
Dr
exp
dtL(r, ṙ, l) Ylm (θf , ϕf )Ylm
(θi , ϕi )4πrf riℏ⟨xf , tf |xi , ti ⟩ = liml,m

1 ( m ) 32
drj
N →∞4π 2πi∆tℏ
j=1

Dr ≡ limN−1∏

これを見ると中心力のもとで l, m は保存していることがわかる。周期境界条件を課して初期位置と終位置に対して全領域積分をすると

[ ∫
]
∫
∫
∑
1
i
dxf Drexp
dtL(r, ṙ, l)
4π
ℏ
l,m
[
]
∫
∫
∑
i
= rf2 drf Dr
exp
dtL(r, ṙ, l)
ℏ
l,m
]
[
∫
∫
∑
i
2
dtL(r, ṙ, l)
=r
Dr
exp
ℏ

∫

dxi dxf δ(xi − xf ) ⟨xf , tf |xi , ti ⟩ =l,m

最後は dxf を Dr に含めた。この作用に対してユークリッド化を行うと次のようになる。

SE =

m
2

(

dr
dτ

)2
+ V (r) +ℏ2 l(l + 1)
2mr2

このパラメータ l は方位量子数であり非負整数をとる。問題点は原点でポテンシャルが発散することだったので以下では l = 1 として遠心力ポテンシャルを設ける。α = 5 としてシミュレーションを行うと各虚時間での粒子の位置 r は図 24 のようになる。

36

図 24 クーロンポテンシャルにおける半径 r の軌跡 : 緑の線はクーロンポテンシャルと遠心力ポテンシャ
ルの和である。比べやすいよう、ポテンシャルは底上げして表示している。

合計のポテンシャルは下に有界であり、半径 r はポテンシャルの底を中心に移動している。大きい r の範囲
ではポテンシャルがなだらかであり、拘束が弱くなっている。

5.3 物理量の計算
α = 1 でエネルギーを求めると E = −0.11(2) となった。Schrödinger 方程式を解いた時の厳密解は式 (85)であり、方位量子数を l = 1 ととっているため、許されるのは n ≥ 2 でありエネルギーの最低値は n = 2 の−0.125 である。シミュレーション結果と誤差の範囲内で一致している。
E=−

α2
mα2
=
−
2n2 ℏ
2n2

(n ≥ l + 1)

(85)

様々な α に対してエネルギーを計算すると図 25 のようになった。α の変更に対して式 (85) のような二乗の依存性を見せていることが分かる。

37

図 25 クーロンポテンシャルにおけるエネルギーの α 依存性 : 各 α = 0.1, 0.5, 1, 5, 10, 20, 30, 40 に対してエネルギーを求めた。曲線は y = axb である。a = −0.25, b = 2

l = 1 の場合、エネルギーが最低になるのは量子数が n = 2, m = +1, 0, −1 の時で、波動関数は次のようになる。

(
Ψnlm (x) = Ψ21m (x) =

1
2a0

) 32

√

r
r
exp(−
)Y1m (θ, ϕ)
2a
3a0
0

このときの半径 r の期待値を計算すると

⟨r⟩ =

5ℏ2
mα

(86)

となる。エネルギーと同様に α を変えて平均半径をシミューレーションすると図 26 のようになった。厳密
解は式 86 にシミュレーション用パラメータを代入して ⟨r⟩ = 5/α となり、図 26 である。α が小さい場合は
平均半径は解析解 (y = f (x)) に近い値になっているが、α が大きくなると解析解からずれ始めポテンシャルが極小になる場所 (y = g(x)) に移っていく。これは α が大きくなると束縛が強くなり、粒子がポテンシャルの底を停留しているからである。つまり、α によって量子から古典への遷移が見られる。

38

図 26 クーロンポテンシャルにおける平均半径の α 依存性 : 各 α = 0.1, 0.5, 1, 5, 10, 20, 30, 40 に対して
平均半径を求めた。緑の曲線 f (α) = 5/α は平均半径 ⟨r⟩ の解析解, 紫の曲線 g(α) = 2/α はポテンシャ
ルが極小になる半径 r を表している。

実際に α = 40 に対して量子的振る舞いを左右するパラメータ ℏ を変えながら ⟨r⟩ /ℏ2 を計算すると図 27 のようになる。解析的には ⟨r⟩ /ℏ2 は式 (86) より

⟨r⟩
5ℏ2 1
5
=
=
= 0.125
2
2
ℏ
mα ℏ
mα

(m = 1, α = 40)となり、ℏ によらない値をとる。図 27 を見ると ℏ を大きくすると、つまり量子性を高めると ⟨r⟩ /ℏ2 は解析
解の 0.125 に近づいていることが分かる。

図 27 クーロンポテンシャルにおける平均半径の ℏ 依存性 : ℏ = 0.1, 1, 2, 3, 4, 5, 6 に設定して r/ℏ2を求めた。緑の直線は解析解 ⟨r⟩ /ℏ2 = 0.125 を表す。

39

5.4 相関関数
ここでは、粒子の位置 r(τ ) の虚時間軸上での相関 Γ を求める。調和振動子においては基底状態に対して相関を求めたことで第一励起状態のエネルギーを探ったが、クーロンポテンシャルでは人為的に l = 1 としているため、系の最低エネルギー状態は |nlm⟩ = |21m⟩ となり、求められる一つ上の励起状態のエネルギーは第二励起状態 n = 2 のものとなる。相関関数 Γ(τ ) は次のように求める。

2

Γ(τ ) = ⟨21m|r(0)r(τ )|21m⟩ − ⟨21m|r|21m⟩

2= ⟨21m|r(0) exp(−Ĥτ /ℏ)r(0) exp(Ĥτ /ℏ)|21m⟩ − ⟨21m|r|21m⟩∑
2=⟨21m|r(0) exp(Ĥτ /ℏ)|n1m′ ⟩ ⟨n1m′ |r(0) exp(−Ĥτ /ℏ)|21m⟩ − ⟨21m|r|21m⟩n(≥2)m′=∑exp(−(En − E2 )τ /ℏ) | ⟨21m|r|n1m′ ⟩ |2 − ⟨21m|r|21m⟩2n(≥2)m′=∑exp(−(En − E2 )τ /ℏ) | ⟨21m|r|n1m⟩ |2

n(≥3)→ exp(−(E3 − E2 )τ /ℏ) | ⟨21m|r|31m⟩ |2(τ → ∞)

相関関数は図 28 のようになった。

図 28

クーロンポテンシャルにおける位置の相関関数 Γ(τ ) と τ

5.5 l ≥ 2 に対するシミュレーション40

6 フーリエ展開による運動エネルギーの計算

41

7 まとめ7.1 結論
調和振動子ポテンシャルについては、Metropolis 法による数値計算によって、シュレディンガー方程式による解析解と近い値を算出することができた。このことは、今回用いた手法が信頼できることを示唆している。
クーロンポテンシャルに対して行ったシミュレーションも同様に確からしい物理量を返した。Woods-Saxon
ポテンシャルについては、存在確率をガウシアンでよく近似できるということがわかった。

7.2 今後の展望・課題
Woods-Saxon ポテンシャルをシュレディンガー方程式より解いた結果と今回のシミュレーションで得られた結果を比較したい。
また、実際のハイパー核におけるパラメーターを用いて計算を行い、実験値との比較を行い、どこまで再現可能かを検証しこの手法の有効性を考えたい。また、中性子数や陽子数が異なる様々な原子核について、基底状
態がどうふるまうのかを調べたい。
今後は、調和振動子ポテンシャルとの対応関係をより深く追求していきたい。

42

謝辞
課題研究を進めることができたのは色々な方々の協力によるところが大きいと思います。
菅沼秀夫さんには研究テーマ全体に関するアドバイスや理論的側面における指導、時には個人的な悩みの相談などを多くのサポートをして頂きました。
前澤祐さんにはシミュレーション手法の知恵や実際の計算機上での計算の実装法などを頂きました。

TA の杉浦巧さんには事あるごとに私達の質問に親切に対応して頂きました。
多くのご指導に大変感謝いたします。

43

参考文献
[1] R. P. Feynman and A. R. Hibbs 著, D. Styer 校訂, 北原和夫 訳 「量子力学と経路積分」 新版, みすず
書房[2] B. Povh, K. Rith, C. Sholtz and F. Zetsche 著, 柴田利明 訳 「 素粒子・原子核物理入門 : Particles andNuclei 」 丸善出版, 改訂新版
[3] 前澤祐 『簡易講義：量子力学のモンテカルロシミュレーション』[4] Mark S. Swanson 著, 青山秀明, 川村浩之, 和田信也 訳 「経路積分法 量子力学から場の理論へ」 吉岡
書店[5] M. Creutz and B. Freedman ”A Statistical Approach to Quantum Mechanics” Ann. Phys. 132, 427462(1981)44