21世紀の現代では、素粒子原子核物理の研究にコンピュータの存在は欠かせないものになっています。このノートは、理化学研究所のスーパーコンピュータ「富岳」など、世界中のスーパーコンピュータで行われている最先端の素粒子原子核研究を紹介するためのノートです。実際にスマホやパソコンを使って、研究を追体験してみましょう!
このノートで取り扱う内容は以下の3つです。
ここで出て来た2つのキーワード「量子場の理論」と「マルコフ連鎖モンテカルロ法」が、このノートのメイントピックです。最後にはこれらを駆使して実際にシミュレーションをしてみます。とても難しい内容ですが、順を追って少しずつ進めていきましょう!!
このノートにある説明はとても大雑把なものになっています。その代わり、このノートの最後に参考文献をまとめてあるので、より詳細に立ち入りたい方はそれらを参考に調べてみてください。また、文中にキーワードを残してありますので、Googleなどで調べてみるだけでも、多くのことが知れると思います。
より詳しく知りたい方は、理学部物理学科への進学を検討してみてくださいね。
このノートでは、いくつかの数値シミュレーション用のプログラムを用意しています。プログラム言語は、pythonとC++です。実行方法については、「Google colaboratory python(C++)」などで検索してみてください。
この章では、素粒子原子核の物理を記述する理論である「量子場の理論(Quantum Field Theory, QFT)について、簡単に紹介します。はじめに断っておくと、2022年現在東北大学のカリキュラムでは、この量子場の理論は大学院講義の講義で扱われる代物です。なので、この章に書かれていることがさっぱり分からなくでも「自分は向いていないんだろうか。。」と悲観しないでください。とはいえ、何も分からないままでは辛いので、感覚だけでも感じられる様に努力しました。一緒に頑張って行きましょう!
キーワード:スケール、物質の階層構造、場の理論、作用関数、経路積分
科学を学ぶ上で非常に大切な概念として「スケール(scale)」というものがあります。物理にとどまらず、あらゆる自然現象を考察する上で欠かせない概念です。ここでは具体的に、「大きさのスケール」に着目してみましょう。
何らかの自然現象を研究しようと考えた時、その対象はある典型的な「大きさ」を持っています。例えば
の様になっています。重要なのは、自然現象をそれぞれが持つ典型的な「大きさ」で分類することができて、そのそれぞれが科学の分野を作っている、つまりその大きさの世界を支配するルールがあるという点です。ここでは「大きさ」を例に取りましたが、エネルギーや温度など様々な物差し(スケール)で現象を分類し、その世界のルールは何かを考えることが物理の第一歩になります。
もちろん物理学を学んでいくと、その様な典型的な大きさを持たない「普遍性」を示す現象にも遭遇します。これは原子の物理を星の物理に適用できるということで、非常に興味深いですが、ここでは割愛します。
ところで、上に列挙した様に、ミクロな世界では素粒子・原子核・原子・分子と何桁ものスケールにわたる階層構造があります。しかし、自然は大きさによってルールをいちいち変えたりするでしょうか?そうでないとすれば、どの様にして階層構造が現れているのでしょうか?近年では、異なる階層の間を橋渡しするための研究も盛んに行われています。
突然ですが、ドラマや小説などはお好きですか?様々な人物が登場し、恋愛や事件を巡って物語が展開して行きます。そのドラマの中で人物たちは各々の役割を担い、様々な関係性を持って関わり合い、1つの脚本の上でドラマを演出してくれます。ミクロの世界(素粒子の世界)でも、事情は同じです。様々な性質を持った素粒子が登場し、ある物理法則に従いながら多種多様の現象を見せてくれます。
量子場の理論は、素粒子などミクロの世界を記述する理論です。その中で、素粒子たちはどの様にして記述されているのでしょうか。
「場」とは何でしょうか。物理学の業界では、空間的・時間的に広がっていれば物理的なものは何でも「場」と呼ばれます。例えば、電磁波や結晶の振動、電子などの素粒子も「場」と呼ばれますが、もっと身近な例はニュースの天気予報を見ればわかりやすいです。
天気予報で見る場の例
大雑把にいうと、空間(+時間)の各点に特定の数字を与えることができればOKです。
$10^{-15}\ \mathrm{m}$以下のとても小さな(大きさのない)素粒子が空間的・時間的に広がっている場であるというのは不思議な感じがするかもしれません。ここでのポイントは、「場の空間的・時間的広がり」と「粒子の大きさ」は一致しないという点です。つまり、場としての広がりはそれ自体が粒子の大きさに対応するのではなく、「その何処かに粒子がいる範囲」としての広がりに対応しています(量子性)。
重要なことは、現代物理学では素粒子は場として記述される、という点です。
キーワード:スカラー場とベクトル場、波動関数、量子性、素粒子標準模型
現代物理学において、様々な物理法則は作用関数(action)と呼ばれる1つの式で表現されることが多いです。作用関数が何者かは割愛しますが、詳しく知りたい方は、解析力学・ラグランジュ形式などで調べてみてください。ここではその代わり、代表的な物理法則の作用関数を紹介します。
量子電磁力学(Quantum Electrodynamics, QED):電子と光の理論。
$$J=\int d^4x \left[-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+\overline{\psi}(i\gamma_{\mu}D^{\mu}-m)\psi\right]$$
量子色力学(Quantum Chromodynamics, QCD):強い相互作用の理論。
$$J=\int d^4x \left[-\frac{1}{2}\mathrm{tr}\ G^{\mu\nu}G_{\mu\nu}+\overline{q}(i\gamma_{\mu}D^{\mu}-m)q\right]$$
超弦理論 タイプIIB行列模型:初期宇宙の理論。 $$J=-\frac{1}{g^2}\mathrm{tr}\left(\frac{1}{4} [A_{\mu},A_{\nu}]^2 + \frac{1}{2}\overline{\Psi}\Gamma^{\mu}[A_\mu, \Psi] \right)$$
一般相対性理論:重力の理論。 $$J=\frac{1}{2\kappa}\int d^4x R\sqrt{-g}$$
さっぱり分からなくても何の問題もありません。現代物理学の代表的な物理法則が作用関数という1つの式で表現されている、ということを感じてもらえればOKです。
キーワード:解析力学、力学系、素粒子標準模型
素粒子(登場人物)と物理法則(脚本)を決めると1つの作用関数(企画)が決まります。後はこの作用関数(企画)に従って計算(実演)するだけです。この計算(実演)の仕方を取りまとめたルールが量子場の理論です。
量子場の理論とは、場の量子的振る舞いを記述する理論です。その始まりは1950年代まで遡り、数多くの成功を収めてきた一方で、その数学的枠組みは未成熟で今もなお発展しています。その枠組みの中で様々な計算手法や定式化がなされていますが、ここでは「経路積分量子化」という手法で、どの様に計算するかを簡単に紹介します。
量子論においてエネルギー$E$などの物理量を計算することは、量子化の枠組みの中でその期待値$\langle E \rangle$を計算することに相当します。登場する素粒子を一纏めに$\Phi$、作用関数を$J^{E}[\Phi]$と書くと、ある手続きのもとで経路積分量子化では物理量$O$の期待値を
$$\langle O \rangle = \frac{\int \mathcal{D}[\Phi]O\mathrm{e}^{-J^E[\Phi]}}{\int \mathcal{D}[\Phi]\mathrm{e}^{-J^E[\Phi]}}$$と計算します。この$\int \mathcal{D}[\Phi]$は少し特殊な積分で経路積分と呼ばれますが、この式の意味するところは今は気にしなくても大丈夫です。重要なことは、1つの作用関数が決まれば、量子場の理論のルールに従ってエネルギーなど実験と比較できる物理量が計算できるという点です。
しかし、残念ながらこの計算は人間には非常に複雑で難しい計算になっています。この困難を打開するために生まれたのが「格子上の量子場の理論」という手法です。これは1974年にケネス・ウィルソンによって提案され、1980年にマイケル・クロイツによってモンテカルロ法を用いた数値計算として成功を収め、現代に至るまで有効活用されています。次の章では、この数値計算に欠かせない技術である「マルコフ連鎖モンテカルロ法」についてみて行きましょう。
キーワード:量子論、解析力学、経路積分法、格子上の量子場の理論
この章では、素粒子原子核の物理をはじめとして様々な分野に応用が効く数値計算技術「マルコフ連鎖モンテカルロ(Markov Chain Monte Carlo, MCMC)法」について、簡単に紹介します。このMCMC法が使われるのは、例えば
という時です。なのでMCMC法は、素粒子原子核物理に限らず、複雑な積分が要求される物性物理や統計的手法が重要な機械学習・金融など、幅広い分野で力を発揮するスゴイ技術です!そんなスゴイ技術ですが、もし天才にしか扱えない技術であれば、あまり嬉しくはありません。しかし、このMCMC法の一番スゴイところは、その発想はシンプルであるにも関わらず、
天才がどう足掻いても解けない問題でも、「根気」さえあれば凡人がコツコツとシミュレーションすることで数値的に解ける
ところです。そんなスゴイ技術であるMCMC法を眺めていきましょう。
キーワード:確率と期待値、区分求積法、モンテカルロ法、確率分布の生成
MCMC法は複雑な積分の計算もできる数値計算技術です。しかし、コンピュータで積分の計算とは何か。まずはここを見て行きましょう。
学校のテストで出てくるような(定)積分の問題を考えてみましょう。
問題:関数$f(x)=x$を区間$a\le x \le b$で積分せよ。 $$\int^a_b x dx = \left[\frac{1}{2}x^2\right]^a_b = \frac{1}{2}(a^2-b^2)$$
一次関数の積分公式を使って、このように積分を評価することができました。しかし、この積分公式が意味するところはなんでしょうか?
実は、積分の計算は関数と区間に囲まれた面積を求めていることに相当しています(区分求積と言います)。今回の場合では、台形の面積を求めることで確認できますが、この考え方は複雑な関数の積分を対象にしても有効です。つまり、刻み幅がとても小さな長方形をたくさん足し合わせることで、積分を表現することができるようになります。
$$\int^a_b f(x) dx = f(a) + f(a+0.0000\dots1)+\cdots +f(b) $$この様に積分も頑張った足し算とみなせるというのが1つ目のポイントです。
人生逆転の夢を見て、宝くじを1口だけ買ったとしましょう。この時もらえる金額の期待値はいくらでしょうか。
宝くじが1等・2等・・・である事象を$A_1,A_2\cdots$、それぞれの当選確率を$P(A_1),P(A_2),\cdots$と書きましょう。この時、ある事象$A$を与えると、その確率を返す関数を「確率分布」と呼び$P(\{A\})$の様に書きます。
また、宝くじは1等からハズレまでどれかには該当するはずです。これは確率の和が1であることを意味しており、確率の定義の1つです。
$$P(A_1) + P(A_2) + \cdots = 1$$また、それぞれの事象でもらえる当選金額を$f(A_1),f(A_2),\cdots$と書くと、期待値は以下の様に計算されます。
$$\langle f(\{A\}) \rangle =\sum_{i=1,2,\cdots}f(A_i)P(A_i) = f(A_1)P(A_1)+f(A_2)P(A_2)+\cdots$$この様に、期待値もその定義からある関数$f(A)P(A)$の足し算とみなせるところが2つ目のポイントです。
今までで確認した通り、積分と期待値はどちらもなんらかの関数の足し算でかけることがわかりました。この2つをつなげることはできないでしょうか?
手始めに区間$a\le x\le b$で定義された確率分布$P(x)$を考えてみましょう。この分布は確率なので次の2つの性質を満たしているはずです。
この時、関数$f(x)$の期待値は
$$\langle f(x) \rangle = f(a)P(a) + f(a+0.0\dots1)P(a+0.0\dots1)+\cdots +f(b)P(b) = \int^a_b f(x)P(x) dx $$とかけます。故に、確率分布について区間$a\le x\le b$で一様であるもの、例えば$P(x)=\frac{1}{b-a}$を用意すれば
$$\langle f(x) \rangle = \frac{1}{b-a}\int^a_b f(x) dx $$と積分と期待値をつなぐことができます。もちろん、一様な確率分布でなくとも、$f(x)/P(x)$の期待値を計算すれば
$$\langle \frac{f(x)}{P(x)} \rangle = \int^a_b \frac{f(x)}{P(x)}\times P(x) dx = \int^a_b f(x) dx $$の様に、積分値を評価することができます。
この様にして積分と期待値の関係を見ることができましたが、実際に期待値はどの様に計算すればいいのでしょうか。答えを言ってしまうと、この期待値はコンピュータを使って求めることができ、この方法は「モンテカルロ法」と呼ばれます。
モンテカルロ法とは、乱数を用いた数値計算手法の総称です。ここではモンテカルロ法を用いて、積分(期待値)をどの様に評価するのかを見て行きましょう。
キーワード:乱数、擬似乱数、確率分布、ヒストグラム、モンテカルロ法
乱数とは、なんらかの確率分布$P(x)$に従ってランダムに生成される数字列のことを言います。具体的にサイコロを使って乱数について考えてみます。
今、サイコロを振って1〜6の内、どの目が出たかを記録していくとします。
例:4→6→1→2→4→3→2→••••••
実はサイコロが一様分布の乱数生成機であることが、次の2点からわかります。
一様分布に従っている。
ここで出たサイコロの目は一様な確率分布に従って生成されていることが知られています。これを確認するためには、どの値が何回出たかを数え上げてグラフにしてみる(ヒストグラムと呼びます)と分かりやすいです。たくさんサイコロを振った時にどの目もだいたい同じくらいの回数が出ていれば、どの値も等確率で生成される一様分布であることが確認できるはずです。
ランダムに生成されている
サイコロは過去にどの値を出したかによらず、次にどの値を生成するかは予測できません。この過去の履歴に依存しないという性質を「ランダム」と呼びます。
このサイコロに限らず、全てのコンピュータには一様分布の乱数を生成する機能が備わっています。具体的に、pythonというプログラム言語でコンピュータに一様分布の乱数を生成させるコードを下に載せました。何回か実行してみてバラバラな値が出力されることを確認してみてください。
#0.0以上以上1.0以下の一様分布から乱数を生成するプログラム
import random
print(random.uniform(0.0,1.0))
0.5439912915686401
(※)ここで確認した乱数を含め、コンピュータで生成される乱数は正しくは「擬似乱数」と呼ばれるもので、厳密な乱数ではありません。非常に大規模な計算ではこの違いが原因となって、思わぬバグを生んでしまう場合があります。気になる方は、擬似乱数・線形合同法・メルセンヌツイスタなどで検索しみてください。
それでは続いて、この「乱数」を使って数値積分を評価してみましょう。まずは一様分布を例にとって、乱数による数値積分を実行してみましょう。評価する数値積分も簡単のために、区間が$0\le x \le 1$の一次関数$f(x)=x$積分を考えます。
この時、数値積分を評価するには次の順序で計算します。
これをプログラムとして実装すると以下の様になります。
#一次関数関数f(x)=xを 0<=x<=1 で積分するプログラム
import random
# N:繰り返しの回数
N = 100
sum = 0
ave = 0
# 積分する関数
def function(x):
return x
for i in range(N):
s = random.uniform(0.0,1.0)
sum += function(s)
#期待値の計算と出力
ave = sum / N
print(ave)
0.5247123710670295
出力した値は公式を使って計算した値から少しずれていることが確認できると思います。このズレは計算が間違っていることを意味していません。実際に、繰り返しの回数を大きくすることでズレを小さくする、つまり精度を上げることができます。モンテカルロ法では、計算に乱数を導入し確率的な手法で積分を評価するため、計算は常にある大きさの誤差(統計誤差)を持つことになります。この統計誤差の評価方法については詳しく紹介しませんが、統計学の技術を使って評価することができます。
キーワード:統計学、統計誤差、母集団分布と標本分布、リサンプリング
この節の最後に、3次元球の体積をモンテカルロ法で計算するc++のプログラムコードを載せておきます。興味のある方は実際にプログラムを書いてみて、好きな積分を評価してみてください。
%%writefile MC_sph.cpp
//=====================//
// モンテカルロ法による半径1の球の体積を求めるプログラム
// 乱数の範囲は 0 <= x,y <= 1
//
// S = {x>0, y>0, x^2+y^2<1}
// P(x,y) = 4 / \pi (x^2+y^2<1)
// 0 (x^2+y^2>=1)
// \int_S dxdy f(x,y) P(x,y) = 1 / (\pi/4) x \int_S dxdy f(x,y)
//=====================//
#include <stdlib.h>
#include <iostream>
#include <random>
int main() {
//=====================//
// In c++, the <random> header is necessary.
// the random-distribution can be obtained by two steps.
// ENGINE : genetates random values
// DISTRIBUTION: molds them as a target distribution
//
// random_device: integer random value used as a seed
// mt19937 : pseudo-random value
//=====================//
int loop_N = 1000; // 何回体積を計算するか
int niter = 10000; // 1回の計算につき何回乱数を生成するか(期待値の和Nに対応)
int rand_max = 10000000;
double ave = 0.;
double sav = 0.;
double err = 0.;
for (int k = 0; k < loop_N; k++) {
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_int_distribution<int> rand_dist(0, rand_max);
int n_in = 0;
double sum_z = 0.0;
for (int iter = 0; iter < niter; iter++) {
double x = (double)rand_dist(mt) / rand_max;
double y = (double)rand_dist(mt) / rand_max;
if (x * x + y * y < 1.) {
n_in += 1;
sum_z += sqrt(1. - x * x - y * y);
}
}
ave += sum_z / n_in / loop_N * 2.0 * M_PI;
sav += sum_z * sum_z / n_in / n_in / loop_N * 2.0 * M_PI * 2.0 * M_PI;
}
err = sqrt(sav - ave * ave);
std::cout << "公式の結果:" << 4. * M_PI / 3. << std::endl;
std::cout << "モンテカルロ法による評価:(中心値、誤差)=(" << ave << "、 " << err << ")" << std::endl;
return 0;
}
Overwriting MC_sph.cpp
!g++ MC_sph.cpp -o RUN.x; ./RUN.x
公式の結果:4.18879 モンテカルロ法による評価:(中心値、誤差)=(4.18867、 0.00524855)
前置きが長くなってしまいましたが、本題「マルコフ連鎖モンテカルロ法とは何か?」に入りましょう。答えを言ってしまうと、マルコフ連鎖モンテカルロ法はお望みの確率分布の乱数を生成するための技術です。
詳細は高度なので立ち入りませんが、MCMC法の最も有名な例であるメトロポリス法について具体例を紹介します。
キーワード:マルコフ連鎖、既約性、非周期性、詳細釣り合い条件、不変分布
ここでは、「ある定められたステップをきちんと行えばお望みの確率分布を生成できる」ことを、実際に分布を生成して確認して行きたいと思います。
まずは目標となる確率分布を決めます。ここでは具体的に、
$$p(x) = \frac{\mathrm{e}^{-J(x)}}{Z},\ \mathrm{where}\ J(x) = \frac{x^2}{2}\ \mathrm{and}\ Z=\sqrt{2\pi}$$の場合を考えます。物理学の用語で$J(x)$は作用(action)と呼び、$Z$は確率の和$\int^{\infty}_{-\infty} P(x)dx=1$を満たす様に決めた係数です。特に今の様な形の分布はガウス分布と名前がついています。
MCMC法の中でもメトロポリス法では、ある初期値$x_0$から次の手続きを踏みながら、次々に$x_{2}, x_{3},\dots$を生成します。具体的に$x=x_{i}$という値を持っているとして、次にどの様な値$x_{i+1}$を出力すべきか考えてみましょう。
メトロポリス法では具体的に$x_{i}$という値を持っている時、次に出力する$x_{i+1}$の候補$y$を
$$y=x_{i}+\xi$$と定めます。ただし、$\xi$は左右対称な確率分布の乱数(例:区間$-1\le x \le 1$の一様分布)です。この$y$を$x_{i+1}$に採用するか否かは次の様に場合わけして考えます。
$p(y)\le p(x)$の時
候補$y$を採択し、$p(y)/p(x)$の確率で$x_{i+1}=y$とする。
$p(y)>p(x)$の時
無条件に候補$y$を採択する。
両者を合わせて、メトロポリス法での採択確率は
$$\mathrm{min}\left\{1, \frac{p(y)}{p(x)}\right\}$$とかけます。つまり、メトロポリス法では
を繰り返しながら、数列$x_0, x_1, x_2,\dots$を生成して行きます。
大雑把ではありますが、以上がMCMC法の例:メトロポリス法の流れです。実際にガウス分布を生成するc++コードを動かして、ヒストグラムを作成してみましょう。一様分布の乱数を使っているにもかかわらず、生成される乱数のヒストグラムはガウス分布を描いています。
%%writefile Metropolis_gauss.cpp
//=====================//
// Metropolis sampling
//
// target: P(x) = exp(-0.5*x^2) / Z
// action = 0.5*x^2
//=====================//
#include <stdlib.h>
#include <iostream>
#include <map>
#include <random>
#include <vector>
int main() {
//=====================//
// In c++, the <random> header is necessary.
// the random-distribution can be obtained by two steps.
// ENGINE : genetates random values
// DISTRIBUTION: molds them as a target distribution
//
// random_device: integer random value used as a seed
// mt19937 : pseudo-random value
//=====================//
int niter = 100000;
int rand_max = 10000000;
double step_size = 0.5;
//=====================//
// Initial input
// x=100 is far from distribution
// -> Need termalization
//=====================//
double x = 0.;
int naccept = 0;
std::vector<double> list(niter);
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_int_distribution<int> rand_dist(0, rand_max);
for (int iter = 1; iter < niter + 1; iter++) {
double backup_x = x;
double action_init = 0.5 * x * x;
//=====================//
// Uniform distribution
// [1] [-0.5:0.5] : satisfy the detailed balance (c=0.5)
// [2] [-0.5:1.0] : break
//=====================//
//[1]
double dx = (double)rand_dist(mt) / rand_max;
dx = (dx - 0.5) * step_size * 2.;
x += dx;
//[2]
// double dx = (double)rand_dist(mt) / rand_max * 0.75;
// dx = (dx - 0.25) * 2.;
// x += dx;
double action_fin = 0.5 * x * x;
//=====================//
// Metropolis test
//=====================//
double metropolis = (double)rand_dist(mt) / rand_max;
if (std::exp(action_init - action_fin) > metropolis) {
naccept += 1;
} else {
x = backup_x;
}
list[iter - 1] = x;
}
//=====================//
// Histogram
//=====================//
double range = 5.25;
double step = 0.1;
std::map<double, int> histogram;
for (double i = -range; i <= range; i += step) {
int number = 0;
for (auto& pl : list) {
if (i - step * 0.5 < pl && pl <= i + step * 0.5) {
++number;
}
}
histogram.insert(std::make_pair(i, number));
}
// for (auto itr = histogram.begin(); itr != histogram.end(); ++itr) {
// std::cout << itr->first << " " << itr->second << std::endl;
//}
//=====================//
// Thermalization
//=====================//
// for (int iter = 1; iter < niter + 1; iter++) {
// std::cout << iter << " " << list[iter - 1] << std::endl;
// }
return 0;
}
Writing Metropolis_gauss.cpp
!g++ Metropolis_gauss.cpp -o RUN.x; ./RUN.x
ここで紹介したメトロポリス法はMCMC法の中でも最も単純なものですが、そのエッセンスは全て詰まっています。より発展したMCMC法を知りたい方は、ギブス法・熱浴法・メトロポリスヘイスティング法・ハイブリッドモンテカルロ法などを調べてみてください。また今回の紹介では、なぜメトロポリス法で目的の確率分布が得られるのかに触れませんでした。この点を知りたい方は、MCMC法の一般論や詳細釣り合いなどをキーワードに調べてみてください。
以上がMCMC法の概要です。今まで紹介してきたことをまとめると、MCMC法で複雑な積分や確率を計算できるというのは
の2点に集約されます。ここだけみるとあまりご利益がないかもしれませんが、解きたい問題が複雑になればなるほど、力を発揮します。例えば素粒子原子核分野の1つである格子量子色力学では、非常に多次元の積分が要求され、単純な問題でもMCMC法を使わなければ宇宙の年齢より長い計算時間をかけても答えを出すことができません。また近年人工知能や機械学習の文脈で目にすることも多いベイズ推定の問題にも応用が可能です。この様に分野の詳細によらず、さらに複雑な問題にも応用可能ということでMCMC法は重宝されています。
キーワード:MCMC法、格子量子色力学、ベイズ推定
今までに紹介してきた知識を使って、実際にシミュレーションをしてみましょう。
今回考察する作用関数は$J=N\mathrm{tr}W\left(\phi\right)$です。$\phi$はNxNのエルミート行列というもので、数字を縦と横に並べた「行列」の中でも特殊な性質を持つものです。また、$W(\phi)$とは$\phi$の多項式です。
具体的に、行列とは(2x2行列の場合):
$$ \left( \begin{matrix} 1 & 5 \\ 7 & 4 \end{matrix} \right) $$の様に書かれ、記号$\mathrm{tr}$はその行列の対角成分を足し合わせる記号です。
$$ \mathrm{tr} \left( \begin{matrix} 1 & 5 \\ 7 & 4 \end{matrix} \right) =1+4=5 $$キーワード:行列、線形代数
この形の作用関数は行列模型( 0次元の場の理論)と呼ばれ、これに似た積分は様々な非常に面白い問題に頻繁に現れます。
今回は特に
$$J=N\mathrm{tr}\left(\frac{1}{2}\phi^2 + \frac{1}{4}\phi^4\right)$$の形に限定してシミュレーションします。
プログラムの簡単な説明はソースコード内に書き入れました。
!g++ --version
g++ (Ubuntu 7.5.0-3ubuntu1~18.04) 7.5.0 Copyright (C) 2017 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%%writefile main.cpp
//=====================//
// Hybrid Monte-Carlo -Multi variables
//
// target: P(\phi) = exp(-J[\phi]) / Z where \phi is NxN hermitian matrix
// action = N Tr { 0.5 x \phi^2 + 0.25 x \phi^4 }
//
// Leapfrog method
// [1] intro n = 0
// phi_ij( 0.5 * dt ) = phi_ij( 0 ) + pi_ij( 0 ) * 0.5 * dt
//
// [2] main n = 1 ~ Nt-1
// pi_ij( n * dt )
// = pi_ij( (n-1) * dt ) - \pdif H / \pdif phi_ij( (n-0.5) * dt ) * dt
//
// phi_ij( (n+0.5) * dt )
// = phi_ij( (n-0.5) * dt ) + pi_ij( n * dt ) * dt
//
// [3] outro n= Nt
// pi_ij( Nt * dt )
// = pi_ij( (Nt-1) * dt ) - \pdif H / \pdif phi_ij( (Nt-0.5) * dt ) * dt
//
// phi_ij( Nt * dt )
// = phi_ij( (Nt-0.5) * dt ) + pi_ij( Nt * dt ) * dt * 0.5
//
//=====================//
#include <stdlib.h>
#include <algorithm>
#include <cmath>
#include <complex>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <vector>
//=====================//
// rand_max: uniform_int_distribution # 変更しないでください
// niter: sample number
// ntau : step number of leapfrog # 変更しないでください
// dtau : step size of leapfrog # 変更しないでください
// ndim : number of variables
// ninit: 1 -> new config., 0 -> old config.
//=====================//
const int rand_max = 1000000;
const int niter = 50000;
const int ntau = 20;
const double dtau = 0.05;
const int ndim = 2;
const int ninit = 1;
int BoxMuller(double& p, double& q, int r_max) {
//=====================//
// In c++, the <random> header is necessary.
// the random-distribution can be obtained by two steps.
// ENGINE : genetates random values
// DISTRIBUTION: molds them as a target distribution
//
// random_device: integer random value used as a seed
// mt19937 : pseudo-random value
//=====================//
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_int_distribution<int> rand_dist(0, r_max);
double pi = 2.0 * asin(1.0);
double x = (double)rand_dist(mt) / r_max;
double y = (double)rand_dist(mt) / r_max;
p = sqrt(-2. * log(x)) * sin(2. * pi * y);
q = sqrt(-2. * log(x)) * cos(2. * pi * y);
return 0;
}
// 作用関数の計算
double calc_action(const std::complex<double> phi[ndim][ndim]) {
double action = 0.;
// \phi^2 パート
std::complex<double> phi2[ndim][ndim];
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
phi2[idim][jdim] = std::complex<double>(0., 0.);
}
}
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
for (int kdim = 0; kdim != ndim; kdim++) {
phi2[idim][jdim] += phi[idim][kdim] * phi[kdim][jdim];
}
}
}
for (int idim = 0; idim != ndim; idim++) {
action += 0.5 * phi2[idim][idim].real();
}
// \phi^4 パート
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
action += 0.25 * (phi2[idim][jdim] * phi2[jdim][idim]).real();
}
}
// 全体に掛ける因子
action = action * ndim;
return action;
}
// ハイブリッドモンテカルロ法に必要な計算(ハミルトニアンの微分)
int calc_delH(const std::complex<double> phi[ndim][ndim],
std::complex<double> (&delH)[ndim][ndim]) {
// \phi^2 と \phi^3 の初期化
std::complex<double> phi2[ndim][ndim];
std::complex<double> phi3[ndim][ndim];
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
phi2[idim][jdim] = std::complex<double>(0., 0.);
phi3[idim][jdim] = std::complex<double>(0., 0.);
}
}
// \phi^2 パート
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
for (int kdim = 0; kdim != ndim; kdim++) {
phi2[idim][jdim] += phi[idim][kdim] * phi[kdim][jdim];
}
}
}
// \phi^3 パート
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
for (int kdim = 0; kdim != ndim; kdim++) {
phi3[idim][jdim] += phi2[idim][kdim] * phi[kdim][jdim];
}
}
}
// ハミルトニアンの微分(delH)の初期化
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
delH[idim][jdim] = std::complex<double>(0., 0.);
}
}
// delH パート
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
delH[idim][jdim] = (double)ndim * (phi[idim][jdim] + phi3[idim][jdim]);
}
}
return 0;
}
// ハイブリッドモンテカルロ法のための計算(ハミルトニアンの計算)
double calc_H(const std::complex<double> phi[ndim][ndim],
const std::complex<double> pi[ndim][ndim]) {
// ポテンシャルエネルギー
double H = calc_action(phi);
// 運動エネルギー
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
H += 0.5 * (pi[idim][jdim] * pi[jdim][idim]).real();
}
}
return H;
}
// ハイブリッドモンテカルロ法のための計算(分子動力学)
int Molecular_Dynamics(std::complex<double> (&phi)[ndim][ndim], double& H_ini,
double& H_fin) {
// initialize pi
std::complex<double> pi[ndim][ndim];
for (int idim = 0; idim != (ndim - 1); idim++) {
for (int jdim = idim; jdim != ndim; jdim++) {
double r1, r2;
BoxMuller(r1, r2, rand_max);
pi[idim][jdim] =
r1 / sqrt(2.) + r2 / sqrt(2.) * std::complex<double>(0., 1.);
pi[jdim][idim] =
r1 / sqrt(2.) - r2 / sqrt(2.) * std::complex<double>(0., 1.);
}
}
for (int idim = 0; idim != ndim; idim++) {
double r1, r2;
BoxMuller(r1, r2, rand_max);
pi[idim][idim] = r1;
}
// initialize delH
std::complex<double> delH[ndim][ndim];
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
delH[idim][jdim] = std::complex<double>(0., 0.);
}
}
H_ini = calc_H(phi, pi);
// First step
// CAUTION: 0.5
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * 0.5 * dtau;
}
}
// 2 ~ Ntau step
for (int step = 1; step != ntau; step++) {
calc_delH(phi, delH);
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
pi[idim][jdim] = pi[idim][jdim] - delH[idim][jdim] * dtau;
phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * dtau;
}
}
}
// Final step
// CAUTION: 0.5
calc_delH(phi, delH);
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
pi[idim][jdim] = pi[idim][jdim] - delH[idim][jdim] * dtau;
phi[idim][jdim] = phi[idim][jdim] + pi[idim][jdim] * 0.5 * dtau;
}
}
H_fin = calc_H(phi, pi);
return 0;
}
// ジャックナイフ用のテンプレートクラス
template <typename T>
class Jackknife {
public:
Jackknife(int n, T m) : N((double)n), M(m){};
T operator()(const T& t) const {
if (N == 1) return t;
return (M * N - t) * (1. / (N - 1.));
};
private:
double N;
T M;
};
int main() {
//=====================//
// In c++, the <random> header is necessary.
// the random-distribution can be obtained by two steps.
// ENGINE : genetates random values
// DISTRIBUTION: molds them as a target distribution
//
// random_device: integer random value used as a seed
// mt19937 : pseudo-random value
//=====================//
std::random_device rd;
std::mt19937 mt(rd());
std::uniform_int_distribution<int> rand_dist(0, rand_max);
//=====================//
// 配位(configuration)の初期化
// MCMC法では、生成した乱数の組みを「配位」と呼びます。
// この配位生成は通常非常に時間がかかるため、時短のために既に生成した乱数を初期値とする場合が
// 多いです。
// このプログラムでは、ninit が1の時は全て0で初期化、0の時は読み込んで初期化します。
// デフォルトは1です。
//=====================//
std::complex<double> phi[ndim][ndim];
if (ninit == 1) {
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
phi[idim][jdim] = std::complex<double>(0., 0.);
}
}
} else {
std::ifstream configufile("configuration.dat");
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
configufile >> phi[idim][jdim];
}
}
configufile.close();
}
int naccept = 0;
int nconf = 0;
double sum_action = 0.;
// confirmation of Thermalization
std::vector<double> list(niter);
// allocate the measured action
std::vector<double> action(niter);
for (int iter = 0; iter != niter; iter++) {
std::complex<double> backup_phi[ndim][ndim];
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
backup_phi[idim][jdim] = phi[idim][jdim];
}
}
double H_ini, H_fin;
Molecular_Dynamics(phi, H_ini, H_fin);
double metropolis = (double)rand_dist(mt) / rand_max;
if (std::exp(H_ini - H_fin) > metropolis) {
naccept += 1;
} else {
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
phi[idim][jdim] = backup_phi[idim][jdim];
}
}
}
list[iter] = calc_action(phi);
if (iter >= 20000 && iter % 100 == 0) {
nconf += 1;
sum_action += calc_action(phi);
action.push_back(calc_action(phi));
}
}
//=====================//
// Thermalization(熱化)
// 生成した乱数の組みが初期値によらないこと(熱化)を確認する必要があります。
// このプログラムでは、乱数の値をテキストファイルに書き出しています。
// 初期値の近くで値が大きくぶれていれば、その値は初期値に強く影響されており、大きな誤差のもと
// になります。実際の計算では、その様な箇所は使いません。
//=====================//
std::ofstream action_monitor("action_monitor.dat");
for (int iter = 0; iter != niter; iter++) {
action_monitor << iter << " " << std::setprecision(16) << list[iter]
<< std::endl;
}
action_monitor.close();
//=====================//
// Save final configuration
//=====================//
std::ofstream configfile("configuration.dat");
for (int idim = 0; idim != ndim; idim++) {
for (int jdim = 0; jdim != ndim; jdim++) {
configfile << phi[idim][jdim] << std::endl;
}
}
configfile.close();
//=====================//
// 平均値と誤差の計算
// MCMC法では、自分自身の情報を使って乱数の生成を行っているので、データ間に相関を持ってしまっ
// ています。そのため、相関を加味した誤差評価(ここではジャックナイフ法)を行うことになります。
//=====================//
double ave_action = sum_action / action.size();
transform(action.begin(), action.end(), action.begin(),
Jackknife<double>(action.size(), ave_action));
double ave = 0.;
double sav = 0.;
for (auto paction : action) ave += paction / action.size();
for (auto paction : action)
sav += (paction - ave) * (paction - ave) * (1. - 1. / action.size());
std::cout << ave << " " << sqrt(sav) << std::endl;
return 0;
}
Writing main.cpp
!g++ main.cpp -o MatrixModel.x; ./MatrixModel.x
0.00865388 0.000632883
このノートでは、「スマホで始める素粒子原子核の物理」と題して
の2つに焦点を当てて、その概要を簡単に紹介しました。最後のシミュレーションでは、それら2つを駆使して、コンピュータを使って物理の問題を解いてみました。これらを通して、最先端の物理研究を感じることができれば、このノートの目的は達成です。
説明の都合上、詳細を省いていたり、厳密な議論をしていない箇所が多く見受けられます。しかし、それら1つ1つの内容が数多くの研究によって、初めて成り立っていることを考えれば、このノート一遍では到底全てを賄えません。アイザック・ニュートンは、論敵であるロバート・フックに宛てた手紙の中で次の様に引用しています。
"If I have seen further it is by standing on the shoulders of Giants."
ー私がかなたを見渡せたのだとしたら、それはひとえに巨人の肩に乗っていたからです。
いつの時代でも、最先端の研究は須く先人達の膨大な研究の基で成り立っています。大学での4年間は、それら1つ1つを腰を据えて考えることのできる絶好の機会です。一歩一歩着実に、知の巨人の肩によじ登っていきましょう。
このノートを作成するにあたって参考にした、または内容を補足するになるべく平易なものをここにまとめました。いくつかの出版物を除き、基本的にはウェブ上で読めるものばかりですので、参考にしてみてください。
量子クラスターで読み解く物質の階層構造:http://be.nucl.ap.titech.ac.jp/cluster/kenkyu.html
菅野礼司、物理教育 第45巻 第4号 (1997):https://www.jstage.jst.go.jp/article/pesj/45/4/45_KJ00005896867/_pdf
EMANの物理学:https://eman-physics.net
【場の理論入門】物理学では相互作用と粒子は場(Filed)という1つの言葉でまとめて記述されます~物理体系の視点の理論の解説〜:https://youtu.be/cXLKwgdgLNI
橋本幸士、講談社、宇宙のすべてを支配する数式」をパパに習ってみた 天才物理学者・浪速阪教授の70分講義:https://bookclub.kodansha.co.jp/product?item=0000147764
R.P.ファインマン 著, 釜江常好 訳, 大貫昌子 訳, 岩波書店、光と物質のふしぎな理論 私の量子電磁力学:https://www.iwanami.co.jp/book/b255808.html
中村純、2013 年度原子核三者若手夏の学校 素粒子パート 場の理論 講義録 格子 QCD シミュレーション入門 —パソコンの上で格子場を楽しむことを目標に—:http://www2.yukawa.kyoto-u.ac.jp/~soken.editorial/sokendenshi/vol18/yamada.pdf
大川正典、石川健一、サイエンス社、格子場の理論入門:https://shop2.saiensu.co.jp/fs/bookshop/2521
花田政範、松浦壮、講談社、ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門:https://bookclub.kodansha.co.jp/product?item=0000343187
区分求積法の意味・例題・公式の証明 -高校数学の美しい物語:https://manabitimes.jp/math/1545
区分求積法:昔の数学者はこうやって面積を計算した!:https://youtu.be/x-xeLzH5SUw
筑波大学オープンコースウェア、確率論:https://ocw.tsukuba.ac.jp/course/math/probability/p-1/
モンテカルロ法:https://www-np.acs.i.kyoto-u.ac.jp/~harada/education/open_courses/mc/index.html
Hacking C++:https://hackingcpp.com
C++マニアック:http://stlalv.la.coocan.jp