符号問題に関する最近の進展

 複素変数の方法 

藤井 宏次
東京大学 大学院総合文化研究科 広域科学専攻 相関基礎科学系
有限クォーク化学ポテンシャルの QCD では、フェルミオン行列式を含む有効作用が複素数になるた
めに、重点サンプリングの方法を直接適用できません。これを符号問題と呼びます。この問題一般
について、力学変数を複素化して解決しようとする野心的な二つのアプローチ、複素ランジュバン方
程式と thimble 上の経路積分、 があります。これらのアイデアの基本と現状を簡単に紹介します。

1

はじめに

1.1

符号問題

有限温度密度の QCD の相構造や相対論的エネルギーでの原子核衝突の物理を第一原理的な方法に基づいて解明することは、この分野の大きな目標の一つです．ハドロンの質量・相互作用や有限温度

QCD の状態方程式については、格子 QCD に基づいた数値的研究が大きな成果を挙げ、他の手法と併せて着実に研究が進展しています．ところが、有限クォーク化学ポテンシャルの状況 (µ ̸= 0) を対象にした途端に、有効作用に現れるフェルミオン行列式 detD(U ; µ) が複素数になるために、格子 QCDで有効だった重点サンプリング (importance sampling) に基づく手法の直接適用ができなくなります．
さらに実時間発展を対象にするとき、格子 QCD が虚時間のユークリッド空間で定式化されていることを思い出せば、QCD 作用に基づいた実時間発展の非摂動的な計算を行うには難しさが予想されます．そもそも、経路毎に作用 S によって決まる位相を平均するという経路積分をうまく定義しなおす
ために、経路積分を虚時間で定式化したのでした．
ヒルベルト空間の次元は粒子数や空間サイズとともに指数関数的に増大します．例えば N 格子点に1/2 スピンが存在している状態空間の次元は 2N です．この理由から、量子多体系問題をハミルトニ
アンの数値対角化に基づいて直接解くことは、系の自由度の増大につれて実質不可能になります．一方、有限温度の量子力学系の分配関数 Z を経路積分表示を用いて表すと、温度方向を虚時間軸とした

1 次元高い空間での古典統計力学系に対応付けられます．その統計平均を計算するための古典配位 (経路) の数は、やはり、系の自由度とともに指数関数的に増大するのですが、ボルツマン重みに基づく
重点サンプリング法が非常に有効にはたらくという利点があります．つまり、

dµ(x) ∼ e−H(x) dx

あるいはdµ(x) ∼ e−S(x) dx(1)

(ハミルトニアン H や作用 S ) のようにボルツマン重みにしたがうサンプルを用いてモンテカルロ (MC)√積分を行うわけです．この MC 積分は、サンプル数 N に対する統計不確かさが必ず 1/ N で抑えられるという意味で、厳密な積分方法です．(ただし別途、離散化などに伴う系統的な不確かさを考察する必要があります．) このサンプリングを実現する方法として、メトロポリス法をはじめ、遷移確率にボルツマン重み ∼ e−S を用いたマルコフ過程が用いられます (動的 MC 法)．この手法は量子モン
テカルロ法としていろいろな分野で応用されています．
ところが、興味ある物理系を経路積分表示した場合に、作用 S が複素数になる場合があります．先に挙げた、有限化学ポテンシャル (µ ̸= 0) の QCD がこの分野では最も有名な例ですが、他にも電子のフェルミオン多体系の問題や θ 項を含む理論などの複素作用を持つ物理的な系が知られています．これらの場合には e−S を確率とは解釈できません．これを符号問題と呼びます．また、実時間の量子発展では指数の肩が純虚数 ∼ eiS となるので、符号問題が最も厳しい例となっています．
これまでに、符号問題を避けつつ有限密度 QCD の情報を引き出す工夫が丹念に行われています [1]．
例えば、被積分関数の絶対値を確率重みとして採用して、複素数の位相部分を観測量として扱う再荷重 (re-weighting) 法が典型です．記号の簡単のために、積分変数を x、被積分関数を w、観測量を Oと書きます．今、w = |w| eiϕ ∈ C という状況なので、|w| 重みで観測量の期待値を評価します：∫⟨O⟩w =

dx O(x)w(x)∫=dx w(x)∫w(x)w
dx O(x) |w(x)||w(x)|/Z|w|⟩|w|⟨O |w|.
=∫w
w(x)⟨ |w| ⟩|w|
dx |w(x)| |w(x)|/Z|w|

一般に、単位体積あたりの熱力学ポテンシャル f と分配関数とは Zw =∫(2)dxw = e−βV fw のように関係します (ここで、β = 1/T は逆温度、V は系のサイズです. Z|w| に対しても同様)．すると (2) の分母に現れる因子は、重み関数 w の位相因子 eiϕ の期待値と解釈され、重み関数 w, |w| の表す二つの系の自由エネルギーの差 ∆f を用いて、⟨

w
Zw⟩|w| = ⟨eiϕ ⟩|w| == e−βV (fw −f|w| ) = e−βV ∆f
|w|
Z|w|

<1

(3)

のように振る舞うことが予想されます．これをシグナルと考えたとき、有意な結果を得るためには統計不確かさがシグナルよりも小さい必要があります：

stat.err. ∼ √1< signal,#smpl

i.e., #smpl > e+2βV ∆f .

(4)

この考察から、低温・熱力学極限に近づくにしたがって、必要な統計サンプル数 N が指数関数的に増√大することがわかります．複素数の被積分関数に対しては寄与の相殺が生じるために、1/ N の精度を保つという MC 法の利点が急激に失われてしまうのです．
この他にも、µ についてのテイラー展開法や虚数化学ポテンシャルの導入、正準アンサンブルの方法、状態密度の方法など、様々な手法で調べられています [1]．
1.2変数の複素化

符号問題を直接に解決する可能性のある野心的な試みとして、力学変数を複素数に拡張する方向性の研究が近年行われています．一つは、確率解釈によらずに動的 MC を行う手法としてランジュバン
方程式を用いる方法です．これに基づく確率過程量子化の方法が 1981 年に提案されました [2] が、その直後に実時間発展や複素作用の場合への適用 [3, 4] が国内も含めて研究されました [5]．複素作用に

ImRe
図 1: 実軸上で定義された複素作用を持つ模型に対して、複素ランジュバン方程式の生成する分布 P (x, y) (灰色領域) と、複素平面に変形された積分路 (点を繋ぐ曲線) の模式図．□は重み関数のゼロ点、○は停留点．
対応して力学変数も複素数に拡張されるので、複素ランジュバン方程式 (CLE) と呼ばれます．しかし、収束性と解の正当化の問題が早くから認識されていました．2008 年ころから改めて再考されて進展があり [6]、2014 年に有限密度 QCD の結果が報告されて注目を集めました [7]．その前後から手法の正当化についての理解が進み [8]、現在 CLE を用いて有限密度 QCD の研究が進められています．
もう一つは積分経路を変形する方法です．解析関数の積分の漸近形を得るために、複素平面内の停留点を通るように積分路を変形して評価する方法を最急降下法と呼びますが、これを多次元空間に拡張するのです [9, 10]．その際、停留点から延びる最急降下空間 J をLefschetz thimbleと呼びます．

Cristoforetti らが thimble J 上での重点サンプリングの着想を先行して発表しました [9] が、私たちも独立に、研究室文献会で紹介された Witten の論文にヒントを得て、thimble J 上の重点サンプリングを考えていました [10]．(同時期に同じことを考える人が在るものです.) そして、その後、手法の拡張が行われ [11]、実時間シミュレーションも試みられています [12]．
複素ランジュバン方程式では、実空間積分を複素空間に拡がった標本分布 P (x, y) での平均として表します。一方で、thimble の方法では、実空間から複素空間内に変形された領域での積分として表します (図 1)。以下では、これらの方法の基本的アイデアと現状を簡単に解説します．

2複素ランジュバン方程式のアプローチ

2.1

ランジュバン方程式とフォッカー・プランク方程式

まず、水に浮かぶ花粉の欠片の乱雑な運動 (ブラウン運動) を記述する確率微分方程式として知られるランジュバン方程式を、簡単に復習します．
「作用 (あるいは熱力学ポテンシャル)」S を持つ力学
変数 x に対するランジュバン方程式は、

dx
∂S
=−
+ η(t)dt∂x(5)です．ここで、第一項はドリフト項と呼ばれ作用 S(x) の停留点に向かう緩和を表し、第二項 η(t) はガウス確率分布に従うミクロな自由度による揺動を表します．十分に時間が経つと、系は作用 S(x)の極小点付近に集まり、揺動力 η(t) による一定の広がりをもって確率的に分布します．この分布の長時間平均が、ボルツマン重み ∼ e−S/T による統計平均に一致するために、揺動力の揺らぎの強さは⟨η(t)η(t′ )⟩ ∼ 2T δ(t − t′ ) と決まります．
経路積分表示は、形式的に、量子系を次元が 1 次元 (虚時間) だけ増えた古典系に対応させます．よって、ボルツマン重み ∼ e−S/ℏ に従って分布する配位で統計平均を取るためには、揺動力揺らぎの強さを⟨η(t)η(t)⟩η = 2ℏδ(t − t′ )(6)にすれば良いことが推測できます (以後 ℏ = 1)．これを確率過程量子化と呼びます [2]．この場合、この時間 t は系の物理的時間ではなく、量子揺らぎをもった配位分布を時系列として生成するために導入したものであり、ランジュバン時間と呼ばれます (以後、区別のために θ と書くことにします)1 ．
確率微分方程式であるランジュバン方程式に従う軌跡 x(θ) は、各時刻毎にある確率分布 ρ(x; θ) を持つと考えられ、それは付随するフォッカー・プランク (FP) 方程式∂∂∂
ρ(x; θ) =(S ′ (x) +)ρ(x; θ) ≡ LT
0 ρ(x; θ)
∂θ
∂x
∂x
(
)(
)
∂
∂
= −e−S/2 (−+ 12 S ′ )(+ 12 S ′ )eS/2 ρ(x; θ)
∂x
∂x

(7)

を満たします (S ′ = ∂S/∂x)．まず、ボルツマン重み ρ(x) ∝ e−S が FP 演算子 LT
0 のゼロであることがわかります．また、二行目の大括弧の中身は自己共役演算子であって、非負の固有値を持つので、T(ゼロ固有値に縮退が無い場合は) 確率分布 ρ(x; θ) = eθL0 ρ(x; 0) → ρ(x) ∝ e−S が保障されます．
2.2

複素ランジュバン方程式 (CLE)

ランジュバン方程式によって時系列を生成する際には、e−S の確率解釈は不要です．故に、複素作
用 S ∈ C をもつ系に対しても、ランジュバン方程式を設定して、問題なく解くことができます．但し、ドリフト項が − ∂S∂x ∈ C のために、θ 時間発展に伴って変数が虚数成分を持つようになるので、変数自体を複素数に拡張 (x → z = x + iy ) する必要があります．これを、複素ランジュバン方程式 (CLE)と呼びます2 ：

dx
∂S
= −Re
+ η(t) ,
dt
∂z

dy
∂S
= −Im
.
dt
∂z

(8)

CLE に対応して、複素空間 (x, y) 上の確率分布 P (x, y; θ) を考えることができ、その FP 方程式を
[∂
]
∂
∂
∂
P (x, y; θ) =(ReS ′ (z) +)+ImS ′ (z) P (x, y; θ) ≡ LT P (x, y; θ)
∂θ
∂x
∂x
∂y

(9)

と書き下すことができます．方程式には「複素」と名前がついていますが、勿論 P (x, y; θ) ≥ 0 です．

CLE の生成する標本を用いた観測量 O(x) の期待値では、O(x + iy) と解析接続したものを、平衡分布 P (x, y; θ → ∞) = P (x, y) を用いて⟨O⟩P =∫
dxdy O(x + iy)P (x, y)(10)を評価していることに相当すると考えられます．一方で、本来の期待値は、複素数値を持つ「ボルツマン重み」ρ(x) ∝ e−S(x) ∈ C を用いた⟨O⟩ρ =∫
dx O(x)ρ(x)

(11)

です．
実数である頻度分布 P (x, y) ≥ 0 を用いて、複素分布 ρ(x) ∈ C についての平均を表すことができるのでしょうか？両者を比べてみると、任意の観測量 O(x) について期待値が等しくなるためには、∫
dyP (x − iy, y) = ρ(x)
(12)
1
2

ϕ(x, t) の場の理論への応用では、ランジュバン時間 θ の時系列として場の配位 ϕ(x, t; θ) を更新していきます．
虚軸方向にも揺動力を導入する自由度が残りますが、ここでは考えません．

0.2

1
4

4
4
6
4
8
4
10

0.8

0.6

iϕ

Re <e >pq

0.15

Re <n>

4

4
4
6
4
8
4
10

0.1

0.4

0.2

0.05

0

0

0

0.25

0.5

0.75

1

0

0.25

µ

0.5

0.75

1

µ

図 2: (左) CLE で得られた粒子密度 ⟨n⟩ と化学ポテンシャル µ の関係．但し、閾値以下 (µ < µc ) の領域．(右) 再荷重で必要になる位相因子平均 ⟨eiϕ ⟩．曲線は ⟨ϕ2 ⟩ の揺らぎを取り込んだ平均場近似．格子サイズ N = 4, 6, 8, 10とともに符号問題が厳しくなり、粒子密度は正しくゼロに漸近する．(ϵ = 5 · 10−5 , m = 1, λ = 1.) [16]が必要十分であることがわかります．
実際、複素ガウス分布 S = 12 (a + ib)x2 (a, b ∈ R) の場合には、複素平面上のガウス型の確率分布

P (x, y) ∝ exp[− 12 (αx2 + 2βxy + γy 2 )](13)

(α = 2a, β = 2a2 /b, γ = 2a(1 + 2a2 /b2 )) が FP 方程式の定常解になっていて、上記の y 積分を行うことで、正しく ρ(x) が再現されることを確認できます．しかし、一般の場合に、CLE によって数値的にのみ生成される P (x, y) が上記の関係をみたすことを確認することはできません．また、ρ(x; θ) のT
満たす FP タイプの方程式：(∂/∂θ)ρ(x; θ) = LT
0 ρ(x; θ) に関しても、複素作用 S ∈ C では L0 が自己共役ではないので、収束の性質を再検討する必要があります．つまり、CLE への拡張は簡明ですが、その正当性が自明ではないのです．
2.3

CLE の有効性と問題点

CLE については、(i) 虚軸方向の無限遠に向かう古典軌道が存在して解が不安定になる場合がある
ことと (ii) 収束解が間違った期待値を与える場合があることが 80 年代当時からすでに問題になっていました．このために、しばらく研究に期間が空いたようです．しかし、2005 年に CLE を用いた実時間量子シミュレーションが試みられた [13, 14] ことを受けて、主に Aarts と Stamatescu の先導によって約 10 年前から CLE の方法の本格的な再検討が始まりました [6]．
例えば、U(1) 電荷を持つスカラー λϕ4 理論に CLE の方法を適用した場合に、閾値以下の化学ポテ
ンシャル µ < µc には系が応答せず、µ > µc ではボース凝縮が発生するという物理的な振舞いが再現されました (図 2(µ < µc ))3 ．複素共役変換は電荷の入替え操作に相当し、電荷の偏りを作り出すために µ を導入したことを考えると、µ ̸= 0 によって作用 S が複素になることを納得できます．図 2 において、符号問題が極めて厳しい (µ = 1 辺りの位相因子平均の値は ⟨eiϕ ⟩ = O(e−200 ) と見積もられています!(N = 10)) 状況でも、CLE 計算が物理的な答えを与える例として注目されました [15, 16]．
上記の懸案 (i) については、その後、ドリフト項が大きくなる領域では θ 時間発展に適応刻み ϵ を用いることで解決されました [17]．一方で、懸案 (ii) については、得られた期待値が満たすべき Dyson-

Schwinger 方程式の観点からも検討されましたが、現在では次節に述べる方向に進んでいます．
化学ポテンシャル µ < µc に 正しく 無反応を示すことを Silver Blaze 現象と名付けています．これは、名探偵シャーロックホームズにおいて「名馬『白銀号』が盗まれたときに、番犬が吠えなかったこと」が犯人の決め手になったことに由来しています．
3

2.4

CLE 正当化の条件

CLE が生成する平衡分布 P (x, y) と S(z) との解析的な関係は一般には分かっていません．CLE が正しい結果を与えるためには、少なくとも P (x, y; θ) による期待値が ρ(x; θ) による期待値に (θ → ∞で) 一致する必要があります．
P (x, y; θ) と ρ(x; θ) は、それぞれ、演算子 L(x, y) と L0 (x) で表される FP タイプの方程式に従い、その形式解は P (x, y; θ) = eθL P (x, y; 0) と ρ(x; θ) = eθL0 ρ(x; 0) と表されます．また、解析関数 O(x + iy)∂∂∂の性質： ∂xO(x + iy) = i∂yO(x + iy) = ∂z
O(z) によって、L(x, y) と L0 (x) について

)
(
∂
∂
′
O(x + iy) = L0 (z)O(z)L(x, y)O(x + iy) = −S (z) +
∂z ∂z

(14)

が成立つことに注意します．すると、初期条件 P (x, y; 0) = ρ(x; 0)δ(y) の仮定の下、二つの期待値が各時刻で一致することが次のように形式的に示されます [19, 20, 21]：∫∫T
dxdy O(x + iy)P (x, y; θ) = dxdy O(x + iy)eθL P (x, y; 0)∫= dxdy eθL O(x + iy)P (x, y; 0)∫= dx eθL0 O(x)ρ(x; 0)∫∫T= dx O(x)eθL0 ρ(x; 0) = dx O(x)ρ(x; θ) .

(15)
(16)
(17)
(18)

ここで、(16) から (17) の等式では、O(x + iy) の解析性と初期条件 P (x, y; 0) = ρ(x; 0)δ(y) を用います．また、(15) から (16) に移る際に、部分積分において表面項が無視できることを仮定しています．

O(x) ∼ eikx のような一般の観測量について表面項の無視が正当化できるためには、分布 P (x, y) が複素空間内で指数関数以上の強さで局在している必要があります [19, 20]．また、フェルミオン行列式

ln det D(x) を含む有効作用 S では、ドリフト項を複素空間に拡張すると 1/ det D(z) が極を持ちます
から、P (x, y) によってその特異点が避けられている必要があります [21, 22]．
さらに、L(x, y) に含まれるドリフト項が y 方向遠方や特異点上で発散し得る状況で、級数を用いて定義される時間推進演算子 eθL(x,y) が有限の θ に対して収束し意味を持つためには、「ドリフト項の大きさ u = |∂S/∂z| についての頻度分布 P (u) が ∝ e−c u (c > 0) 以上の強さで抑制されること」が必要であり、さらに、この条件が部分積分正当化の十分条件にもなっていると指摘されています [23]．
経験上、CLE 計算の時系列は或る一つの平衡分布 P (x, y) に収束するので、上記の条件を満たして収束すれば観測量を正しく評価していると期待されます4 ．一変数模型やカイラルランダム行列模型、

Thirring 模型などを検討し、分布 P (x, y) の振舞と結果の正誤には凡そ良い対応が見られています．
カイラル相転移とフェルミオン行列式のゼロとは密接な関係があり、相転移領域では分布 P (x, y)が行列式のゼロ近傍に分布する状況がよく生じます。そのため、カイラル転移領域では CLE が適用できなくなることが多く、注意が必要です。上記のように CLE の正当化条件の理解が整ってきたので、現在は標本分布 P (x, y) の振舞いを検討することで結果の正当性を判断しつつ、解析が進められています．
4

但し、FP 方程式に従う ρ(x; θ) が正しい「複素ボルツマン重み」ρ(x) ∼ e−S に収束しているかという点は、別に確認する必要あるように思います．

2.5

QCD のシミュレーション

2014 年に、Sexty によって、軽いクォークを含む QCD の低密度から高密度 µ/T = 0 ∼ 8 にかけての CLE の計算結果が発表されて、非常に注目されました [7]．高温相で小さい格子サイズ (43 × 4,
β = 5.7 > βc ) でしたが、符号問題のある µ ̸= 0 の条件でfull QCD の計算が実行できることには驚かされました．
これに先立って必要だったのが gauge cooling という手法です；SU(N ) ゲージ理論では、変数の複素化 (U ∈ SU(N) から SL(N, C)) に伴って、標本分布 P (U ) が、SU(N ) の空間から離れて拡って分布してしまう問題がありました．先に指摘したように CLE の正当化には、拡がった分布は好ましくありません．そこで、作用の持つ SL(N, C) 不変性に基づいて、ステップ毎に配位を SU(N ) 空間近くに移動させる変換 (gauge cooling) を挿入したのです。これによって、β が大きい (格子間隔が小さい)
場合には、分布が局在して収束した結果を得ることに成功しました [24]．
現在では、CLE の正当化条件をチェックしつつ、より低温で大きな格子規模の QCD への適用が進められています [7, 25, 26, 27, 28, 29]．例えば [28, 29] では、スタガードフェルミオン (Nf = 4) を用いて 83 × 16, β = 5.7 や、さらに大きいサイズの格子でのシミュレーションが実施されています．しかし、CLE の正当化条件を満たすために β を大きく (格子間隔を小さく) 選ぶ必要があり、低温領域に接近する上での問題になっているようです．

3 Thimble 上の経路積分と一般化
複素積分における最急降下法を振り返った後、thimble 上の積分とその一般化について説明します．
3.1

振動積分に対する最急降下法

1 次元の振動積分の例として、
Ai(x) =

1
2π

∫ ∞

s3

ei(sx+ 3 ) ds

−∞

(19)

で定義される Airy 関数を取り上げます．この被積分関数はまさに振動関数です．しかし、x を大きい√パラメタとするときに、指数部分 h(s) = i(xs + s3 /3) の停留点 s∗± = ±i x を考え、s∗+ の周りで√h(s) を二次近似して積分を評価 (鞍点法) すれば (h′′∗+ = −2 x)、
∫ ∞
2 3/2
1 ′′+ 22 3/2
1Ai(x) ≃ e− 3 x
e 2 h∗ s ds = √ 1/4 e− 3 x
(20)
2 πx
−∞
のように漸近形が得られることを、物理数学で習った方もあると思います．
以上の鞍点法は近似評価でしたが、もとの積分路を s = s∗ を通るように連続変形した経路を採って積分を行えば、積分値を変えることなく厳密な値が再現されます．その代表として停留点 (臨界点ともいう) から流れ出る最急降下経路：)(dsdh(s) ∗=−dτds(s(τ ) → s∗ as τ → −∞)(21)

を採ることができます．ここで、s = sR + isI , h = hR + ihI のように実部と虚部に分けて書きなおすと、Cauchy-Riemann 関係式から、

dsR
∂hI
∂hR
=− R =− I,
dτ
∂s
∂s

dsI
∂hI
∂hR
=− I =+ R
dτ
∂s
∂s

(22)

図 3: (左) Airy 関数の最急降下経路 (実線) と上昇経路 (破線)．♡ は停留点 (x が複素数の場合). (中)実軸上での被積分関数 (x = 1)．(右) 最急降下経路に沿う線上の被積分関数の振舞い (横軸は Re(s))．
がわかります5 ．つまり、この流れは、実部 hR が最も減少する勾配流れであると同時に、虚部 hI をハ
ミルトニアンとして保存する軌道とも解釈できます．最急降下経路では、被積分関数が振動することなく端点でゼロに収束するので、積分の収束性が極めて良いことがわかります (図 3)．実軸に沿った積分 (19) は収束性が明確ではないので、むしろ最急降下経路を含めて積分を定義するのが適当です．

(停留点 s∗− に付随する最急降下経路は、Airy 微分方程式のもう一つの独立な解に対応しています．)
3.2

Lefschetz thimble 上の経路積分

1 変数の場合の最急降下経路を、多変数に拡張したものが Lefschetz thimble です6 ．それは、N次元複素空間上の停留点 zσ から湧き出す勾配流れ (holomorphic gradient ow)
dzi
=−
dτ

(

dh(z)
dz

)∗

(23)

が覆う点の集合として定義される実 N 次元多様体 Jσ です．1 次元の場合と同様に、積分領域の変形として、実空間上の積分を thimble に沿った積分として等価に書き直すことができます：
∫
∫
∑
−S[x]
−i ImS[zσ ]dz e−ReS[z] .
dx e
=
e
Rn

(24)

Jσ

zσ

zσ は停留点．複数の thimble Jσ の集合が実空間 RN と等価になる場合には zσ で和を取ります．
これは等式変形なので、thimble 上の積分を評価すれば正解が得られます．式 (24) を眺めると、

thimble 上の積分では ImS[z] = ImS[zσ ] が一定に留まる一方、e−ReS は停留点で最大値を持ち、そこから外れるに従って単調減少するので、重点サンプリングによる積分法が効率よく適用できそうです．
この観察が thimble 上の経路積分を考える基本です．
この方法では、観測量 ⟨O⟩ の期待値を

∫
⟨O⟩ =

5

dx e−S[x] O[x]
∫
=
−S[x]
Rn dx e

Rn

∫

∑

dz
−i ImS[zσ ]zσ eJσ dξ dξe−ReS[z] O[z]∫dz−i ImS[zσ ]−ReS[z]zσ e
Jσ dξ dξ e

∑

(25)

CLE との対応を考えて h(z) = −S(z) とおいてみると、勾配流れの方程式と CLE とでは、実成分の流れの符号が反転している、という関係があります．
6
Thimble とは、裁縫のときに硬い布に針を指すときに使う指ぬきだそうです．図 3(右) を多次元に拡張して想像される、釣鐘上の形を指しての命名でしょうか．

図 4: Thirring 模型の thimble 構造 (一様な場の配位に限定 ; 左から µ = 0.6, 1.2, 1.8. β = 3, m =1, N = 4). Thimble (太線) は停留点 (黒丸) から伸びて、行列式のゼロ点 (白丸) あるいは z = ±π ± i∞に終点を持つ [32]．細線は勾配流れの流線 (図 5 の流れを定義する)．背景の濃淡は ReS の大きさを表す．
のように表します．ここで、ξ は thimble 上の点を表す N 個の実パラメータで、私たちは、停留点 zσ近傍の座標を (23) 式によって τ − τ0 だけ移すという対応で決めました [10]．|dz/dξ| ∈ C はそのヤコビ行列式です．従って、thimble の上で MC 計算を行うためには、次のような手続きが必要になります：

(i)
(ii)
(iii)
(iv)

実空間に等価な thimble の集合 {Jσ } の同定、複数の thimble の寄与にわたる平均、
一般に曲がった空間である thimble の座標 ξ の構成、ヤコビ行列式 (∈ C) の評価．

0+1 次元 Thirring 模型: thimble 構造の例として、(0+1) 次元の Thirring 模型の場合を図 4 に示

します．Thirring 模型はベクトル型 4 点結合のフェルミオン模型であり、4 点結合をほどくための補助場 z を導入した後にフェルミオン場を積分することで、分配関数を∫∫−Sb (z)Z = dz det D(z; µ) e
= dz e−Seﬀ (z)(26)という、フェルミオン行列式が現れる形に書くことができます。この模型は (ゼロ温度で) カイラル
対称性の自発的破れを持ち、化学ポテンシャル µ の増加とともにカイラル凝縮が消失してフェルミオン密度が増大します (図 6 を参照)．図 4 のパラメタでは、µ = 1.2 辺りが相転移領域に相当します．

µ = 0.6 では、原点付近から実軸上 ±π に延びる曲線が thimble です．µ = 1.2 になると、虚軸上の停
留点から左右に thimble J0 が延びて、それが行列式のゼロ点 (有効作用 Seﬀ (z) では無限大の点) に
終端点を持ちます [30]．そして、この作用無限大の点を挟んで連結する隣の thimble J±1 が ±π + i∞
まで伸びています．これらが実軸 (−π, π) に等価な積分路を構成することになります．µ = 1.8 では、
虚軸上の停留点から延びる thimble J0 が広い領域を覆っています．
上記の課題のうち、(iii), (iv) は、流れ方程式を数値的に解いて thimble を構成せざるを得ない煩雑さを反映したものです．特に行列式の評価は、数値計算ではよく知られるようにコストが大きいので、アルゴリズム実装に際して問題を引き起こします．また、thimble 上では、ヤコビ行列式が複素数になることに注意します．これを残った符号問題と呼び、再荷重の方法で考慮しています [10]．
課題 (i), (ii) は、MC アルゴリズムを作成する上で困難を引き起こします．全ての thimble を同定して、その中から実空間積分に等価な集合 {Jσ } を決定することは、ゲージ理論など一般の理論では困難

0.75

T=0.5

T=0.050.65

T=0.01Im✄ ✂ 0.70

0.60

-1.5

-1.0

-0.5

0.0

Re ✁

0.5

1.0

1.5

図 5: 主要な thimble J0 の停留点での接空間 (実軸に平行な直線) を流れ方程式の時間パラメータ T で写した像の領域 DT (T = 0.01, 0.05, 0.5). N = 8, m = 1, g 2 = 1/2 (β = 1), µ = 1.6. [11] から引用).
です．また、必要な thimble の集合が判明した場合も、ゼロ点 (作用の発散点) を跨いで繋がる thimbleがあれば、それらの寄与を平均する上で困難が生じます．無限の障壁に阻まれた両者間を遷移する MC法を構成することが自明でないのです．さらに、それぞれの thimble は異なる位相 exp(−i ImS[zσ ]) を持つのですが、特に一次相転移領域では、複数の thimble の寄与の相殺が重要になることが指摘されています [31, 32]．これを大域的符号問題と呼んでいます．
3.3

Thimble から離れる

解析関数の複素積分で重要な性質は、積分路の選択に積分値が依存しないことでした．そこで、もとの実 N 次元空間 RN を初期条件に選び、勾配流れによって時間 τ だけ流したときの像として定まる、実 N 次元領域 Dτ を考えます．この像 Dτ は RN の連続変形なので、どの τ を採っても、RN に等価な積分領域とみなすことができます．そして、τ → ∞ では、自動的に RN に等価な thimble の集合に漸近していきます．つまり、この写像は (i) 模型の thimble 構造の予備知識無しに RN に等価な積
分領域を定義し、(ii) その領域内には無限障壁がないので全体を MC 法で遷移できる、という性質がもつことがわかります．これは前の課題 (i), (ii) を解決するものです．Thirring 模型について、虚軸上停留点の接空間を初期条件に選び、流れ方程式で時間パラメータ T (= 0.01, 0.05, 0.5) だけ移動した像が図 5 に具体的に描かれています (N = 8, m = 1, g 2 = 1/2 (β = 1), µ = 1.6) [11]．
この勾配流れに沿う線上での、ReS の単調増加性と ImS の一定性から、τ の増大とともに、ImS[zσ ]と等しい虚部を持つ点 x は停留点 zσ に流れ込み、RN のそれ以外の領域は ReS = ∞ の領域に向かって引き延ばされて流れていきます．このように、ImS = 一定の領域が伸張することで符号問題が緩和されると同時に、異なる thimble に向かって漸近する領域の間には ReS による障壁がだんだん増大します．これは、RN と thimble {Jσ } 上の経路積分をうまく内挿する方法になっています [11]7 ．
図 6(左) では、Thirring 模型について、接空間から得られた領域 DT 上で MC 積分を行った場合の配位の出現頻度が、各 T 毎に示されています．横軸は、領域 DT 上の点を流れによって引き戻した対∑応する接空間での座標 ϕ̂ = n ϕn /N ∈ (−π, π) です．接空間上でのシミュレーション (つまり T = 0)
では、十分に広い範囲の配位がサンプルされています．前述のように、停留点に流れ込む 1 点を除いて、接空間の各点は勾配流れによって引き伸ばされてゼロ点に吸い込まれて行きます．従って、T が7

領域 Dτ の点は実変数 x ∈ RN で指定することができ、変数変換のヤコビ行列式は |dz(τ )/dx| に置き換わります．

図 6: (左) 接空間から時間パラメタ T = 0, 0.1, · · · , 0.5 で決まる空間上で行った MC サンプリングの頻度分布．(中) T = 0, 0.5 それぞれの空間で MC を行った場合のカイラル凝縮の評価値．曲線が厳密
解．外れている点が T = 0.5 での MC の結果．(右) tempering によってサンプルされる T = 0.5 空間
内の配位の頻度分布．(N = 16, m = 1, g 2 = 1/2, µ = 1.0) [34] から引用.
大きい程、接空間内の狭い領域のみが MC 積分に効くようになります．T = 0.2 の場合に 3 つの局在分布が見られることは、3 つの thimble (虚軸を横断するもの J0 、その両脇 J±1 ) に近い領域がサンプルされていることを表しています．しかし、T > 0.3 では J0 に対応する 1 つの局在しか見られなくな
ります．これは、積分領域が thimble に近づくに従って、ReS による領域間の障壁が高くなった結果、標本が 1 つの谷間にハマっていること (エルゴード性の破れ) を示しています．T = 0 と 0.5 で評価したカイラル凝縮の結果を図 6(中) に示しました．確かに、T = 0.5 では正しい結果を再現できません．
この状況を解決するために、tempering という統計力学の手法が適用されました [33, 34] [35]．異なる T に対応する積分領域 DT の中でそれぞれ MC を行いつつ、時折、それらの間でサンプルの交換をMC で行うという二段構えの手法です．小さい T の MC を経由することで thimble 間の遷移を加速して、短かい MC 時間でエルゴード性を保った標本の平衝化を達成させようというアイデアです．結果の配位の頻度分布が図 6(右) に示されています．(左) と異なって、隣の thimbleJ±1 の寄与が取り入れられていることを確認できます．Tempering を行うことで、各 thimble に対応する複数の谷間の配位を正しくサンプルでき、その結果、厳密解を再現できることが報告されています [33, 34]．
3.4

最適化

流れ方程式を用いて積分領域を変形する際には、パラメタ T の選び方に任意性があり、どの T の値
での領域 DT が積分の評価に最も有効なのか厳密な基準はありません．一旦、複素空間での積分領域の選択自由度を認めると、むしろ、この不定性を積極的に利用して、(thimble 構造から一旦離れて) 符号
問題を最小化するように積分路を選ぶという立場が可能になります．この問題を多次元空間内における積分領域選択の最適化問題として捉え、機械学習を応用する研究も最近進められています [36, 37]．

4

実時間発展
2008 年頃から始まる CLE 再考のきっかけは、先行した実時間シミュレーションの試みでした [13, 14]．

実時間量子系のシミュレーションが可能になると、例えば、身近には、保存量密度の長時間相関を調べることで輸送係数を直接的に評価できたり、更には非平衡現象の研究につながったりします．
準備された初期状態 |in⟩ での二点相関関数 ⟨in|ϕ(t1 )ϕ(t2 )|in⟩ = ⟨in|U (t0 , t1 ) ϕ U (t1 , t2 ) ϕ U (t2 , t0 )|in⟩のような、初期状態からの時間発展に依存する観測量の有限温度期待値は、初期時刻 t0 から始まって

t0 に戻る閉じた時間順序経路 C にわたる生成汎関数∫

Z[J] = Tr(e−βH TC (ei C J·ϕ )) =∫Dφ ei(S+J·φ)

(27)

P.B.

を使って表すことができます．作用 S が位相として現れていますが、この場合にも CLE や一般化された thimble の方法に基づいて経路積分を評価することが考えられます [13, 14] [38][11, 39]．

CLE の方法では、4 次元時空の場の配位 φ(x, t; θ) をスタート条件として準備し、ドリフト項 −i∂S/∂φを用いた CLE に従って、複素数に拡張された場を更新していきます．文献 [14] では、有限温度 λϕ4理論について、対象とする時間発展の長さが t ≲ 1/T では熱平衡解が得られるのですが、長時間発展(t ≳ 1/T ) を対象にした CLE では、物理時間 t に関する並進対称性の破れた解になってしまうことが観察されています (T は温度)．また、SU(2) ゲージ理論の実時間相関関数の評価においても不安定性が見られます (gauge cooling 無し)．
(一般化された)thimble の方法では、勾配流れ方程式を利用して、振動積分の積分領域を実 φ(x, t)の空間から変形された複素空間内の領域に移して積分を行うことによって、符号問題を和らげます．
試験的な 1+1 次元 λϕ4 理論への適用では、t ≲ 4/T まで解くことができていますが、t > 4/T では積
分領域についてエルゴード性の破れの問題が現れるようです [39]．
また、どちらの方法も物理的な結果を得るまでの収束が遅く、大きな標本を必要とするようですから、アルゴリズムの改良が不可欠になります．

5

まとめ
複素ランジュバン方程式の方法と、最急降下経路 thimble にヒントを得た一般化 thimble の方法について、これまでの流れを簡単に解説しました．最後にそれぞれの方法の特徴をまとめておきます．
5.1

CLE

変数の複素化を許してランジュバン方程式を解くのみなので、アルゴリズムは非常に簡潔です．計
算コストも、実変数の場合に比べて 2 倍になるだけです8 ．さらに模型によっては、符号問題の厳しい状況でも完全な結果を返す点で非常に優れています．しかしながら、分布 P (x, y) の振舞いによっては、適用できない状況があることが分かっているので、正当化条件に照らして、得られた結果の正しさを常に吟味する必要があります．相転移領域を調べる際には、変数変化や、模型に変形パラメタを導入するなど、問題回避に向けた工夫が試みられています。
5.2

一般化 thimble

もとの積分の等式変形なので、評価できれば正しい結果を与えます．しかし、計算コストに課題があります；MC ステップ毎に曲がった空間の数値的構成とヤコビ行列式の評価が必要になります．特に、自由度 N について行列式の評価は N 3 のコストになるので、アルゴリズムの改善が工夫されています [39]．(CLE では複素空間に拡がる分布 P (x, y) を扱い、ヤコビ行列が不要であることと対照的です．) また、エルゴード性の破れについて注意する必要があります．行列式から生じる残った符号問
題と、複数 thimble に対応して大局的符号問題が存在していることから、N を大きくした場合にそれらがどのように問題を残すのか、応用上の明確な理解がまだありません．
8

場合によっては gauge cooling が必要になります．

力学変数の複素化という点が共通しており、また、共に流れの停留点が分布や積分領域を特徴付けることから、CLE の標本分布と thimble 構造には相関も観察できます．しかし、厳密なことはわかっていません．いずれにせよ、それぞれのアプローチの利点を活かしつつ研究を進めることによって、符号問題を乗り越える緒につながるのではないか、と期待されています．

謝辞 関連する研究を進めてきた菊川さんをはじめ共同研究者の方々、いろいろな場面で議論などに協力してくださった研究者の方々に感謝します．

参考文献[1] P. de Forcrand, PoS LAT 2009, 010 (2009).
[2] G. Parisi and Y. s. Wu, Sci. Sin. 24, 483 (1981).
[3] J. R. Klauder, J. Phys. A 16, L317 (1983).
[4] G. Parisi, Phys. Lett. B 131, 393 (1983).
[5] P. H. Damgaard and H. Huel, Phys. Rept. 152, 227 (1987).
[6] G. Aarts and I. -O. Stamatescu, JHEP 0809, 018 (2008).
[7] D. Sexty, Phys. Lett. B 729, 108 (2014).
[8] E. Seiler, EPJ Web Conf. 175, 01019 (2018), and references therein.
[9] M. Cristoforetti et al., Phys. Rev. D 86, 074506 (2012); Phys. Rev. D 88, no. 5, 051501 (2013).
[10] H. Fujii, D. Honda, M. Kato, Y. Kikukawa, S. Komatsu and T. Sano, JHEP 1310, 147 (2013).
[11] A. Alexandru, G. Basar, P. F. Bedaque, G. W. Ridgway and N. C. Warrington, JHEP 1605, 053 (2016).
[12] A. Alexandru, G. Basar, P. F. Bedaque, S. Vartak and N. C. Warrington, Phys. Rev. Lett. 117, no. 8,
081602 (2016).
[13] J. Berges and I.-O. Stamatescu, Phys. Rev. Lett. 95, 202003 (2005).
[14] J. Berges, S. Borsanyi, D. Sexty and I.-O. Stamatescu, Phys. Rev. D 75, 045007 (2007).
[15] G. Aarts, Phys. Rev. Lett. 102, 131601 (2009).
[16] G. Aarts, JHEP 0905, 052 (2009).
[17] G. Aarts, F. A. James, E. Seiler and I. O. Stamatescu, Phys. Lett. B 687, 154 (2010).
[18] H. Fujii, S. Kamata and Y. Kikukawa, arXiv:1710.08524 [hep-lat].
[19] G. Aarts, E. Seiler and I. -O. Stamatescu, Phys. Rev. D 81, 054508 (2010).
[20] G. Aarts, F. A. James, E. Seiler and I. -O. Stamatescu, Eur. Phys. J. C 71, 1756 (2011).
[21] K. Nagata, J. Nishimura and S. Shimasaki, PTEP 2016, no. 1, 013B01 (2016).
[22] G. Aarts, E. Seiler, D. Sexty and I. O. Stamatescu, JHEP 1705, 044 (2017).
[23] K. Nagata, J. Nishimura and S. Shimasaki, Phys. Rev. D 94, no. 11, 114515 (2016).
[24] E. Seiler, D. Sexty and I. -O. Stamatescu, Phys. Lett. B 723, 213 (2013).
[25] G. Aarts, E. Seiler, D. Sexty and I. O. Stamatescu, Phys. Rev. D 90, no. 11, 114505 (2014).
[26] Z. Fodor, S. D. Katz, D. Sexty and C. Török, Phys. Rev. D 92, no. 9, 094516 (2015).
[27] K. Nagata, H. Matsufuru, J. Nishimura and S. Shimasaki, PoS LATTICE 2016, 067 (2016).
[28] S. Tsutsui, Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki and A. Tsuchiya, arXiv:1811.07647 [hep-lat].
[29] Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki, A. Tsuchiya and S. Tsutsui, arXiv:1811.12688 [hep-lat].
[30] T. Kanazawa and Y. Tanizaki, JHEP 1503, 044 (2015).
[31] Y. Tanizaki, Y. Hidaka and T. Hayata, New J. Phys. 18, no. 3, 033002 (2016).
[32] H. Fujii, S. Kamata and Y. Kikukawa, JHEP 1511, 078 (2015) Erratum: [JHEP 1602, 036 (2016)];JHEP 1512, 125 (2015) Erratum: [JHEP 1609, 172 (2016)].
[33] M. Fukuma and N. Umeda, PTEP 2017, no. 7, 073B01 (2017).
[34] A. Alexandru, G. Basar, P. F. Bedaque and N. C. Warrington, Phys. Rev. D 96, no. 3, 034513 (2017).
[35] R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 57, 2607 (1986); K. Hukushima and K. Nemoto,J. Phys. Soc. Jpn. 65, 1604 (1996); R. M. Neal, Stat. Comput. 6, 353-366 (1996); For review, D. J. Earland M. W. Deem, Phys. Chem. Chem. Phys. 7, 3910 (2005).
[36] Y. Mori, K. Kashiwa and A. Ohnishi, Phys. Rev. D 96, no. 11, 111501 (2017).
[37] A. Alexandru, P. F. Bedaque, H. Lamm and S. Lawrence, Phys. Rev. D 96, no. 9, 094505 (2017).
[38] Y. Tanizaki and T. Koike, Annals Phys. 351, 250 (2014).
[39] A. Alexandru, G. Basar, P. F. Bedaque and G. W. Ridgway, Phys. Rev. D 95, no. 11, 114501 (2017).