統計力学と計算機シミュレーション

川島 直輝 （都立大理）
２００３年９月２７日 久保記念シンポジウム

モンテカルロ法
1777 Buffon の針
1942～ 分子拡散の問題へのモンテカルロ法の応用(Fermi, Metropolis, Ulam, etc)
1945 ENIAC (Pennsylvania Univ.) 完成 … 弾道計算
1946 「モンテカルロ法」の命名 (Ulam による)
1952 MANIAC (Los Alamos) 真空管
1953 Metropolis, Rosenbluth2, Teller2: メトロポリス法
1955 Fermi-Pasta-Ulam （非線形問題研究の発端）
1957 Alder-Wainwright: 剛体球系における相転移
1959 Fosdick: ３Ｄ Ising model のシミュレーション
1962 Handscomb: 級数展開による量子系の数値計算
1976 Suzuki: 経路積分による量子モンテカルロ
1987 Swendsen-Wang: クラスタ更新

Buffon の針
でたらめに落とした長さ１の針が間隔１で引かれた平行線のどちらかと交わる確率は２／π．
2 （針を落とした回数）
×（針が線と交わった回数）π=

モンテカルロ法＝決定論的な問題を乱数を利用して解く方法（「モンテカルロシミュレーション」とはいえない）
3408 回の試行により π=3.1415929 (Lazzerini 1901)

中性子拡散への応用
ENIAC を用いた
モンテカルロシミュレーション
by Maria G. Mayer (1951)
d = min − λiµ log r µ)

rµ :

µµ(

一様乱数 (0 < r < 1)
λ : 材質 i 中でのチャネル µ に関する平均自由行路µi

モンテカルロシミュレーション＝確率論的問題を乱数を利用して解く方法

The criticality problem for a liquid
metal breeder reactor (reported in
U.S. Department of Commerce,
Applied Mathematics, Series 12:1920)

Metropolis
Metropolis, Nicholas Constantine (1915-1999)

1941 Ph.D. (Experimental Physica, University of Chicago)
1943 Los Alamos National Laboratory1945 Return to Chicago
1946 ``Monte Carlo method'' was named by S.Ulam1948 Return to Los Alamos
1948 The first actual Monte Carlo calculations on ENIAC
1952 MANIAC became operational
1953 The Metropolis algorithm published(Metropolis, A. Rosenbluth, M. Rosenbluth,A. Teller, and Edward Teller)
1957 MANIAC II
1957 Return to Chicago
1965 Return to Los Alamosメトロポリス法（１）Metropolis, Rosenbluth, Rosenbluth, Teller, Teller 1953

重みつきサンプリング∑W (S ) Q(S ) ∑W (S ) Q(S )P(S ) = const ⇒ Q =≈∑W (S )∑ W (S )（１）重みなし

iSiiiiSi（２）重みあり∑W (S ) Q(S ) ∑ Q(S )P(S ) ∝ W (S ) ⇒ Q =≈∑ W (S )
∑1i

Siii

Si

メトロポリス法（２）Metropolis, Rosenbluth, Rosenbluth, Teller, Teller 1953

Markov 過程と詳細つりあい
極限分布が重み W (S ) に比例するようなマルコフ過程をつくればよい．
それには遷移確率が詳細釣り合い

P (S ′ | S ) W (S ) = P (S | S ′ ) W (S ′ )を満たしていればよい．

メトロポリス法による最初の計算
Metropolis, Rosenbluth, Rosenbluth, Teller, Teller 1953

２次元剛体円盤（Ｎ＝２２４）のシミュレーション

実際にはその後Alder-Wainwright の計算により転移があるとされているが，この計算では分からなかった．

Ising Model のシミュレーション
Fosdick 1959

２元合金の実験（Cu3Au）との比較

Ｎ＝５×５×５

CF: Peierls の議論（１９３６），Onsagerの解（１９４４）

非線形振動子Fermi, Pasta, and Ulam 1955

非線形性のためにエネルギー等分配状態に近づいていくと期待されていたが，そうはならなかった．

剛体球系の相転移
Alder-Wainwright 1957

引力のない系でも液相固相転移が起きる．

相転移に関してKirkwood 1935,Born-Green 1959などの予想はあったが，近似の性質から信頼できるものではなく，実際に転移が起きるかどうかはこの数値シミュレーションによって初めて分かった．

スピングラス
Kawashima-Young 1995

典型的な「難しい」問題．
そもそも相転移があるのかないのかすらすぐには分からない．
（解析的な手法ではいまだに分かっていない．）
g: 秩序変数の分布がどの程度ガウス関数からずれているかをあらわす無次元量

モンテカルロ法による級数展開
Handscomb's method (1962, modified by Lyklema in 1982)
分配関数をスピン自由度に関する和でなく，互換演算子の様々な積 "C'' に関する和の形に書いて，C をモンテカルロ法で更新していく．

H = m J ∑ S i ⋅ S i +1 − H ∑ S izi( J > 0)i

J
∑ Ei ,i +1 + const
2 i(Eij = ``swap' ' operator)= H0 −Z =∑∞

n =0

βnn!

( ) = ∑ (± 1) ∑W (C )

Tr H

n∞

n=0nC

( C =n )n(βJ )W (C ) =Tr (e − βH E p E p L E p E p )0nn −1n!
C ≡ ( p1 , p2 , , pn ) , pk = (ik , jk )L21

モンテカルロ法による級数展開
Handscomb's methodC ≡ ( p1 , p2 ,L, p , p , L, p ) ↔ C ′ ≡ ( p , p , L, p , p , p , L, p )
k

k +1

newn

1

2k

次数を上げるとき：
新しく挿入する互換の候補と，それを挿入する場所の候補をランダムに選んで，逆向きの過程との間に，詳細釣り合いが成り立つような確率で挿入．
次数を下げるとき：削除する互換の候補を選んで詳細釣り合いが成り立つような確率でそれを削除．

k +1

長所： 統計誤差以外の誤差はない．
問題点：
１） 応用範囲が狭い．
２） フラストレーションがなくても負符号問題がある．
３） 臨界緩和がある．

n経路積分に基づく量子モンテカルロ
Suzuki 1976dτ L (S , S )∫0
Z = ∫ DS eβ−≈ ∫ DS e−&L (S k , S k )
M∑

β M

&k =1

≈ ∑∏ w(S k +1 , S k )M= ∑ W (S )
S

k =1

S

世界線のパターンS : 世界線の局所更新
W (S ') P(S ' | S ) = min1 ,W (S ) Swendsen-Wang アルゴリズム
Swendsen-Wang 1987

S

G

S’

反転しない
反転しないクラスタ
しないクラスタ
反転する
するクラスタ
反転する
クラスタ
詳細つりあい条件：

P (G | S )W (S ) = P (S | G )W (G )

ループアルゴリズム
Evertz-Lana-Marcu 1993

経路積分表示上のクラスタアルゴリズム

グラフの基本構成要素（反強磁性ハイゼンベルクモデルの場合）

グラフ割り当て確率:
P() J2 ∆τ↑↓ =連続虚数時間極限
経路積分離散化による誤差の解消

P()J↑↓ = ∆τ2ρ()J↑↓ =2大きいスピンのシミュレーション
Kawashima-Gubenatis 1995
Harada-Troyer-Kawashima 1998

Split-spin trick
2S1
Sˆ i ⇒σˆ i , µ∑2 µ =1

L

Tr ( ) = Tr ( PˆSσL)
σ 1,1 σ 1, 2 σ 2,1σ 2, 2

An SU(N) model
H=Hamiltonian:

Jαβ′J(r)J∑βα (r )N (r ,r ')[J , J ] = δ J − δ Jαβγδαδγβγβαδ

α , β , γ , δ = 1,2,L, N

Representation:
Santoro et al. (1999), present
Affleck (1985), Arovas-Auerbach (1988)
Read-Sachdev (1990)

“Magnetic” Correlation FunctionSπ ,π ≡

M π2,π−dr=L(−1)S Z (r ) S Z (0 )d∑Lr

N = 3 : ms2 = 0.085(1)(ms = 0.291(2) )N = 4 : ms2 = 0.025(2)(ms = 0.158(6) )N = 5 : ms2 = −0.0003(8)N = 6 : ms2 = 0.0004(5)Finite Neel order for N=2,3,4, but no Neel order for N>=5

確率的級数展開再び（１）
Sandvik-Kurkijarvi 1991
Z =∑∞「空箱」を挿入次数の制限を除く
ハミルトニアンを分割
スピン自由度を導入

n =0=∑Lβn( ) = ∑ n ! Tr (H )+ o(en!

Tr H


−1

 n=0 n !  n =∑{σ k }L!

=∑β (L − n ) !
nL!

CS∑

L!

L)

(σ k = 0,1)∏k LTr  H σ k 
 k =1∏L k =1

Tr  ∏ H g k 

β n (L − n ) !

= ∑W (C )G−L n

Tr  H σ k {σ k }:∑ σ k = n  k =1β n (L − n )!

= ∑∑nn=0β n  LGβnLnLL() G ≡ g1 , g 2 , , g n 
0,1,2,,=gN
 kH 

H0 ≡ 1


LS H S(C ≡ (S , G ) , S ≡ (S , S ,L, S ))S1 H iL S L S L H iL−1 S L −1212i1n1

確率的級数展開再び（２）Sandvik-Kurkijarvi 1991

Z = ∑ W (C )W C( )≡Cβ n (L − 1)!

C ≡ (S , G )L!

S1 H iL S LL
S 2 H i1 S1 L
L-1

スピン（黒丸）配置
S ≡ (S1 , S 2 ,L, S )n

ハミルトニアン配置
G ≡ ( g1 , g 2 ,3L, g )

2n

1

「空き箱」

無限次数極限「空き箱」の導入によって，級数展開モンテカルロ法の配位空間は時間が離散化された経路積分と本質的に同じものになる．
（「トロッター数」＝「展開次元の最大値」）

経路積分にもとづく方法でつかえる状態更新法が全て使える．つまりループによる更新が可能．
連続時間極限（無限次数展開極限）も同様にとることができて，結果は経路積分から出発して得られるアルゴリズムと完全に一致する．

おわりに数値シミュレーションによる「理解」とは何か？
１）実験との定量的な比較．
２）対立する理論的予想の「審判」．
３）「意外な」現象の発見．
（量子系多体問題について「意外な」現象が発見されるのはこれから？）