\documentclass[10pt]{article}
\usepackage{psfig}
\usepackage{multicol}

\newlength{\figwidth}
\setlength{\figwidth}{0.8\textwidth}

\pagestyle{empty} 
%\pagestyle{plain} 
%\pagenumbering{arabic}
\renewcommand{\labelenumi}{(\arabic{enumi})}
\setlength{\oddsidemargin}{-0.61cm}
\setlength{\evensidemargin}{-0.61cm}
\setlength{\textwidth}{17.1cm}
\setlength{\topmargin}{-1.8cm}
\setlength{\textheight}{25.5cm}
%\setlength{\footskip}{1.5cm}

\setlength{\columnsep}{0.7cm}

\newcommand{\sign}{\mbox{sign}}
\newcommand{\sn}{\mbox{sn}}
\newcommand{\cn}{\mbox{cn}}
\newcommand{\dn}{\mbox{dn}}
\newcommand{\sech}{\mbox{sech}}
\newcommand{\tr}{\mbox{tr}}
\newcommand{\Tr}{\mbox{Tr}}
\renewcommand{\Re}{\mbox{Re}}
\renewcommand{\Im}{\mbox{Im}}


\title{\vspace*{-1.2cm}
\Large {\bf 格子量子色力学の大規模シミュレーション}
\vspace{-0.3cm} }

\author{
\normalsize  松古 栄夫 \\
\normalsize  高エネルギー加速器研究機構 (KEK)
\vspace{0.3cm}}

\date{
{\Large\bf Large-scale simulations of lattice quantum chromodynamics}
\smallskip\\
\normalsize by \\
\normalsize Hideo Matsufuru 
\vspace{-0.2cm}}

\begin{document}

\maketitle

\setlength{\baselineskip}{4.2mm}

\centerline{\sf\bf ABSTRACT}
%\abstract{
Lattice Quantum Chromodynamics (QCD) is one of the scientific
fields which require most large-scale numerical simulations.
QCD is the fundamental theory of the strong interaction among
quarks and gluons inside hadrons.
Because of difficulty in analytic calculations, numerical
simulations of lattice QCD are important to explore the
properties of hadrons and determine their matrix elements.
In this report, we explain fundamental aspects and algorithms
of lattice QCD, and how its large-scale simulation is performed
on massively parallel computers with an example at KEK.
JLQCD Collaboration is performing dynamical simulations with
overlap fermions, which have theoretically elegant features while
require high numerical cost.
We also introduce international and domestic activities,
ILDG and JLDG respectively, to share the configuration data
produced in lattice QCD simulations.
%}
\bigskip

\begin{multicols}{2}

\setlength{\baselineskip}{5.3mm}

\section*{\normalsize 1. はじめに \vspace{-0.3cm}}

格子量子色力学(Quantum Chromodynamics, QCD)に基づく数値
シミュレーションは、大規模計算を可能にする計算機やアルゴリズムの発展
と歩調を合わせて進展してきた。
筑波大学におけるCP-PACSプロジェクトなどの専用計算機の開発や、
IBM Blue Gene/L の母体となったQCDOC機の開発など、高性能計算の発展
に果たした役割も大きい。
ここでは、格子QCDシミュレーションの原理とアルゴリズムを概説し、
KEKで行っているJLQCD Collaboration のプロジェクトを例として、
大規模シミュレーションがどのように行われているかを紹介する。

シミュレーションアルゴリズムの視点からは、場の方程式を離散化した
ものであり、境界値問題(周期的に取る場合が多い)の解法が中心である。
計算時間のほとんどは、線形方程式の解法に費される。
また、Monte Carlo法を構成する際に、発展方程式が導入される。
ベクトル計算にも、大規模並列計算にも適した性質の問題である。

\section*{\normalsize 2. 格子QCDの原理
\vspace{-0.3cm}}

原子核を構成する陽子や中性子は、クォークと呼ばれる素粒子が、
強い相互作用によって束縛されて出来ている。
量子色力学(QCD)はこの強い相互作用を記述する理論であり、
数学的にはSU(3)群に基づいたゲージ理論として構成される。
これはクォークは「色」の自由度を持つと考え、その変化によって
相互作用を記述することを表し、その際に相互作用を媒介するゲージ粒子
としてグルーオン場が導入される。
QCDを解析的に解いて、陽子や中性子(ハドロンと呼ぶ)の性質を示す
ことは困難である。
これは、量子電磁気学の場合やQCDの高エネルギー現象の場合には、
物理量を結合定数によって展開する摂動論が有効であるが、低エネルギーの
QCDではこの結合定数が大きくなりすぎ、摂動論が破綻するためである。
従って、ハドロンの示す多彩な性質をQCDに基づいて理解し、フレイバー物理
などに必要なハドロンの散乱振幅を精密に計算するためには、
数値的に計算を行う方法が必要である。
現在のところ格子QCDシミュレーションは、低エネルギー領域でのQCDの
性質を調べるための唯一の一般的枠組を与え、
広範に研究が行われている\cite{Aoki2005}。

格子QCDは、4次元Euclid空間上に定義された場の理論であり、
経路積分法を用いて量子化される。
これによって統計力学系と同じ形となるため、Monte Carlo法に
よって数値的な計算が可能となる。
まず時空間を4次元Euclid空間の有限格子で近似し、各格子点(サイト)に
クォーク場を、サイトとサイトを結ぶリンク上にゲージ場
(グルーオン場)を表すリンク変数を配する。
クォークはカラー自由度(R,G,B)を持つと考え、この3次元空間内
での回転によって相互作用を記述する。
リンク変数は$3\times 3$複素行列で、ユニタリー群 $SU(3)$ に属し、
サイトと方向($x,y,z,t$)によってラベルされる。
クォーク場は、カラー自由度に加えてスピノール(スピン$\uparrow$,
$\downarrow$と粒子・反粒子を表す)4成分を持つので、$3\times 4=12$
の自由度を各サイトで持つ。

格子上の場の理論を定義するには、クォークとグルーオンについてそれぞれ
作用を構成する必要がある。
これらは、格子間隔$a\rightarrow 0$の連続極限でQCDと一致するもので
あればよい。
ゲージ場の作用は、格子上の最小のループである四角形に沿って、
リンク変数を掛け合わせ、そのトレースをとったもの(プラケット)
で表される。
プラケットの係数として結合定数(の逆)が導入され、これとクォークの
質量がインプットパラメターとなる。

クォーク場の作用としては、ここではもっとも単純な構造を持つ、
Wilson 作用を考える。
\begin{eqnarray}
  S_F = \sum_{x,y} \bar{\psi}(x) D[U]_{x,y} \psi(y), 
\end{eqnarray}
ここで、$\psi$ はクォーク場を表し、$x$, $y$ はサイトである。
Wilsonフェルミオンのディラック演算子$D[U]_{x,y}$は
次のように与えられる。
\begin{eqnarray}
  D[U]_{x,y} &\!=\!& \delta_{x,y} - \kappa {\textstyle \sum_{\mu}}
 \left[
  (1-\gamma_\mu) U_{\mu}(x) \delta_{x+\hat{\mu},y} \right.
\nonumber\\
 & & \hspace{0.8cm}
+ \left.
 (1+\gamma_\mu) U_{\mu}^\dag(x-\hat{\mu}) \delta_{x-\hat{\mu},y}
 \right]
\label{eq:Wilson_action}
\end{eqnarray}
ここで $U_{\mu}(x)\in\,$SU(3)はリンク $(x,\mu)$ 上のリンク変数、
$\mu=x,y,z,t$は方向を表わし、$\hat{\mu}$は$\mu$-方向の単位ベクトルである。
格子間隔は $a=1$とおいた。
$\gamma_\mu$ は $4\times 4$の行列で$\gamma$-行列と呼ばれ、クォークの
スピノール構造を規定する。
$\kappa$はクォークの質量に関係したパラメターである。

経路積分による量子化によって、物理量$O$の期待値は、次のように表される。
\begin{eqnarray}
 \langle O \rangle  &=& \frac{1}{Z}
 \int {\cal D}U {\cal D}\psi {\cal D}\bar{\psi}
  \,  O(\bar{\psi},\psi,U)
 \nonumber\\
 & &  \times \exp(-S_G[U]-S_F[\bar{\psi},\psi,U]),
\end{eqnarray}
ここで分配関数$Z$は$\langle 1 \rangle =1$となるような規格化因子である。
この表式は、$O$の期待値は、あらゆるゲージ場とクォーク場の値について、
$\exp(-S_G[U]-S_F[\bar{\psi},\psi,U])$の重みをかけて積分すること
によって表されることを表す。
従って、$\exp(-S_G[U]-S_F[\bar{\psi},\psi,U])$の確率で現れる
場の配位を生成することができれば、それらの上で計算した物理量の
平均値として、上の期待値が表されることになる。
これがMonte Carlo法による計算の原理である。

クォーク場は反交換するグラスマン数として扱われるが、このような数を
計算機上で直接扱うのは難しいため、手で積分してしまうと、
\begin{equation}
 \int {\cal D}\psi {\cal D}\bar{\psi}
  \exp(-S_F[\bar{\psi},\psi,U]) = \det D[U]
\end{equation}
となる。
Gauss積分の性質を使って、
\begin{equation}
  \det D[U] =
  \int {\cal D}\phi {\cal D}\bar{\phi}
   \exp(-\phi D^{-1}[U] \phi)
\end{equation}
と書けば、通常のベクトルによってフェルミオンの効果を取り入れること
が出来る。

クォーク場を含む物理量は、クォークの伝播関数を用いて表すことができる。
例えば、クォークの伝播関数を組み合わせてハドロンの相関関数を
構成し、その時間方向への伝播の振る舞いを調べることによって、
ハドロンの質量や崩壊定数を求めることができる。
ハドロン散乱の行列要素などについても、同様の計算が可能である。

クォークの伝播関数は、つぎのような線形方程式を解いて得られる。
\begin{equation}
 {\textstyle \sum_{y}} D[U]_{x,y} x_y = b_x,
\label{eq:linear_eq}
\vspace{-0.2cm}
\end{equation}
Eq.~(\ref{eq:linear_eq})において、$b$はある与えられたベクトル
(サイト$\times$カラー$\times$スピノールの自由度を持つ)、
$x$が解である。
この方程式は巨大連立線型方程式であり、$D[U]_{x,y}$は疎行列であるため、
共役勾配法などによって解を求める。
この線型方程式系の解法が、格子QCDシミュレーションの計算時間のほとんど
を占める場合が多い。

このようにして格子QCDシミュレーションで求めた物理量を、実験値と
比較できるようにするためには、格子間隔 $a$ で次元を持つ量に
直さなければならない。
計算の際に手で指定するパラメターは、結合の強さとクォークの質量で
あり、格子間隔は何らかの物理量、例えば陽子の質量など、
を実験値と比較することによって決める。
現実の物理では、軽い $u$, $d$ クォークとやや重い $s$クォークが
真空偏極に寄与する。
これらを実際に含み、$u$, $d$ クォークに対応する軽い質量での
計算は難しいため、やや重いところからの外挿が行われる。
精密な計算結果を得るためには、これらの軽いクォーク質量への外挿や
連続極限$a\rightarrow 0$への外挿などによる系統誤差をコントロール
する必要がある。


\section*{\normalsize 3. 格子QCDシミュレーション
\vspace{-0.3cm}}

格子QCDシミュレーションは、通常次のような手順に沿って行われる。
\smallskip\\
{\bf (1) リンク変数の配位の生成}\\
リンク変数$U_\mu(x)$の生成は、通常 Hybrid Monte Carlo 法という
アルゴリズムを用いて行う。
この方法では、あるリンク変数の配位(全てのリンク変数がある数値を
持ったような状態)から出発し、分子動力学的にリンク変数を発展させて
ゆくことによって、統計的に独立な配位を作って行く。
この手法では、クォークが真空において対生成・対消滅をくり返す、
真空偏極の効果を取り入れることができる。
真空偏極の効果を計算する際には、クォークの伝播関数を求めるのと同じ
計算を、発展の各ステップで行うことが必要となる。
このため計算コストが非常に大きくなる。
\smallskip\\
%
{\bf (2) クォーク伝播関数の計算}\\
前節で述べたように、クォークの伝播関数を求めることは巨大な線型方程式
を解くことに対応する。
リンク変数の配位を生成する際にも各ステップで同じことが必要であり、
結局この解法が、計算時間の大部分を占める場合が多い。
従って、共役勾配法などの逐次解法のアルゴリズムを改良することは、
シミュレーション時間の短縮に直結する。
\smallskip\\
%
{\bf (3) 物理量の計算}\\
クォークの伝播関数が求まれば、それらを組み合わせて種々の物理量が
計算できる。
例えばハドロンの質量、ハドロン崩壊過程の崩壊定数や形状因子などを計算し、
その統計誤差、系統誤差の見積りを行う。
素粒子物理では、KEKで行われているBファクトリー実験などの
高エネルギー衝突実験に関係するハドロンの遷移確率の計算が重要である。
これらは標準理論を検証し新しい物理を探すために、高精度の値が
必要とされている。
このためには、特に系統誤差をコントロールすることが重要となる。
また最近では、QCDに基づいて核子(陽子、中性子)間に働く核力を理解する
試みも進展している。
\smallskip

格子QCDを並列計算機で扱う場合には、格子を分割し、
各部分格子を各ノードに配する。
現在行われている格子は$16^4$--$24^4$程度のサイズのものが多いため、
$O(1000)$以上の並列ノードで計算する場合には、各ノード上の部分格子
は各方向に数個程度となる。
Eq.~(\ref{eq:Wilson_action})が示すように、となり合うサイト
での相互作用を伝えるために、演算性能に加えて通信性能が重要な要素となる。


\section*{\normalsize 3. KEKのプロジェクト
\vspace{-0.3cm}}

実際の大規模シミュレーションの例として、
高エネルギー加速器研究機構(KEK)のスーパーコンピュータシステムを
用いて我々が行っている、JLQCD Collaboration
のプロジェクト\cite{web_jlqcd}を紹介する。
国内では他に、筑波大学で大規模な研究プロジェクト(PACS-CSプロジェクト)
が行われている。

最近の格子QCDの大きな発展として、格子上でのカイラル対称性の理解が
進んだことがある。
カイラル対称性とは、クォークの質量が非常に軽い場合に成り立つ対称性
であり、現実世界ではこの対称性が自発的に破れることによって、
陽子や中性子の持つ大部分の質量が生じると説明される。
しかしながら格子上でこの対称性を持つ理論を構成することは長い
間困難であった。
最近このカイラル対称性の格子上での理論的な理解が進み、
その結果カイラル対称な格子理論の定式化が可能となった。
そのような理論の一つとして、我々はオーバーラップ・フェルミオン
を採用した。
この定式化は、理論的に優れた性質を持っているが、これに従って
シミュレーションを行うには従来の100倍以上の計算力が必要となる。
これはこの理論が、Eq.~(\ref{eq:Wilson_action})の演算子の
符号関数を含む形で表されるからであり、この符号関数の評価に時間が
かかるためと、Eq.~(\ref{eq:Wilson_action})の固有値がゼロとなる
点で符号関数が不連続に変化するためである。
KEKなどの最新鋭スーパーコンピュータの導入に加え、アルゴリズム
などの改良を経て、我々はオーバーラップ・フェルミオンを用いた
プロジェクトを実行中である\cite{web_jlqcd}。
基本的な演算は Eq.~(\ref{eq:Wilson_action})と同じ演算子から
構成されるため、この演算子をベクトルに作用させる演算を
高速化することによって、シミュレーションの効率を向上すること
ができる。

KEKでは、スーパーコンピュータシステムとして、日立SR11000
(理論演算性能 2.15TFlops, 総メモリ容量512GB)および、
IBM System Blue Gene
(理論演算性能 57.3TFlops, 総メモリ容量5.12TB)
からなる複合システムを2006年3月より運用している\cite{web_keksc}。
このシステムは、素粒子・原子核物理、物性物理、加速器内の
ビーム状態などの
シミュレーションに利用されているが\cite{web_ohgata-s}、
計算時間の大きな部分を格子QCDが占めている。

以下では、Blue Gene で行っている計算について紹介する。
Blue Gene システムは 10ラックからなり、1ラックは1024ノード
(2048プロセッサコア)から構成される。
各ノードは4MBのL3キャッシュを持ち、これを2つのプロセッサが共有して
いる。
ネットワークの単位は1/2ラックで、$8\times 8\times 8$の
トーラス型であり、1ノード内の2プロセッサも二つのノードのように
扱うことが出来るため、$8\times 8\times 8\times 2$
の4次元的トポロジーとなる。
格子サイズは$16^3\times 32$を用いた。
いくつかのパラメターで並列に計算を進めるのが効率的であり、
必要となる演算量、格子サイズ、結果を得るまでの時間などを考慮して
1/2ラックまたは1ラック毎に、一つの格子でのシミュレーションを行っている。

Blue Gene での高速化は、各プロセッサの倍精度複素数演算を
ダブルFPUによって高速化することと、ノード間の通信を最適化する
ことによって行う。
これらを最適化することによって、最も単純なEq.~(\ref{eq:Wilson_action})
のタイプの線型問題に対して、オンキャッシュの場合にはピーク性能の
30\%近くの性能を得ることができた\cite{Doi2006}。
実際のプロジェクトではその他の演算も加わるため、15\%程度の性能で
計算を行っており、現在も改良を続けている。
物理的な結果についても成果が出始めており、
最近では、カイラル対称性が自発的に破れていることを、
QCDから直接的に検証したことが大きな成果として挙げられる
\cite{Fukaya07}。

\section*{\normalsize 4. 配位データの共有にむけて
\vspace{-0.3cm}}

真空偏極の効果を含むシミュレーションを行うには、大きな計算資源が
必要である。
一方、リンク変数の配位を一度作ってしまうと、それを用いて種々の計算が
可能となる。
配位データの公開、共有化を進めることは、同じ計算を行う無駄を防ぎ、
作られた配位を効率的に利用するために有効である。
このような共有化を推進するために、国際的な組織である
ILDG(International Lattice DataGrid)が活動を行っている\cite{web_ildg}。
ILDGはいくつかの地域グリッドから構成され、
半年毎に行われるTV会議によって方針を決定している。
また、メタデータ及びミドルウェアに関する二つのワーキンググループが
活動している。

日本国内では、JLDG(Japan Lattice DataGrid)という組織が活動を行っており、
地域グリッドとして配位データをILDGに提供すると共に、
高速ネットワークSINET3 を通じた遠隔サイト間での高速なデータ転送や、
その効率的運用のためのシステムの開発などを行っている\cite{web_jldg}。
後者としては、産業技術総合研究所で開発されたGfarmを用いたデータ
グリッドの構築を進めている。
現在は、筑波大、KEK、京都大学基礎物理学研究所、大阪大学核物理研究
センター、広島大、金沢大がJLDGに参加している。


\section*{\normalsize 5. 結論
\vspace{-0.3cm}}

格子QCDシミュレーションで行っている計算の概略を紹介した。
最近の理論的進展は著しく、また計算機、アルゴリズムの発展に伴い、
第一原理であるQCDに基づいた定量的な計算が可能になっている。
特にクォーク質量については、従来大きかった外挿による不定性が
著しく改善され、現実的クォーク質量での計算に近付いてきている。
これは素粒子物理学に必要な行列要素の精密な計算に不可欠であり、
新しい物理の探索にも大きく貢献できると期待されている。
また核力の計算などにも応用が始まっている。
シミュレーションで生成された配位データを有効利用するための、
データの共有化の枠組が整備されつつある。



\begin{thebibliography}{99}
\vspace{-0.2cm}

\bibitem{Aoki2005}
  青木慎也「格子上の場の理論」(シュプリンガー・フェアラーク東京, 2005).
\vspace{-0.3cm}

\bibitem{web_jlqcd}
  JLQCD Collaboration, {\tt http://jlqcd.kek.jp/}.
\vspace{-0.3cm}

\bibitem{web_keksc}
  {KEKスーパーコンピュータシステム}, {\tt http://scwww.kek.jp/}.
\vspace{-0.3cm}

\bibitem{web_ohgata-s}
  {KEK大型シミュレーション研究}, {\tt http://ohgata-s.kek.jp/}.
\vspace{-0.3cm}

\bibitem{Doi2006}
土井淳, 寒川光, 松古栄夫, 橋本省二,
情報処理学会 論文誌(トランザクション)コンピューティングシステム
47 No.SIG7 (2006) 114.
\vspace{-0.3cm}

\bibitem{Fukaya07}
  H.~Fukaya {\it et al.}  [JLQCD Collaboration],
  %``Two-flavor lattice QCD simulation in the epsilon-regime with exact chiral
  %symmetry,''
  Phys.\ Rev.\ Lett.\  {\bf 98} (2007) 172001.
%  [arXiv:hep-lat/0702003].
  %%CITATION = PRLTA,98,172001;%%

\bibitem{web_ildg}
 ILDG, {\tt http://www.lqcd.org/ildg/}.

\bibitem{web_jldg}
 JLDG, {\tt http://www.jldg.org/}.

\end{thebibliography}

\end{multicols}

\end{document}

