"シミュレーションで探る基礎科学：素粒子の基本法則から元素の生成まで"2021年度第２回HPCIC計算科学フォーラム 2022/03/28, Zoom online

CI(Configuration Interaction)計算とその発展的手法による大規模原子核構造計算
東京大学理学系研究科 附属原子核科学研究センター
清水 則孝
阿部喬（理研）、宇都野穣（原研）、大塚孝治(東大理・理研）、角田佑介（東大原子核セ）、吉田亨（RIST）
T. Otsuka, T. Abe, T. Yoshida, Y. Tsunoda, N. Shimizu, N. Itagaki, Y. Utsuno, J. Vary, P. Maris and H.
Ueno, accepted.
T. Abe, P. Maris, T. Otsuka, N. Shimizu, Y. Utsuno, and J.P. Vary, Phys. Rev. C 104, 054315 (2021).
N. Shimizu, Y. Tsunoda, Y. Utsuno, and T. Otsuka, Phys. Rev. C 103, 014312 (2021).

富岳成果創出加速プログラム"シミュレーションで探る基礎科学：素粒子の基本法則から元素の生成まで"

2

Contents• 研究目的• 原子核CI計算とは• モンテカルロ殻模型・準粒子真空殻模型– エネルギー分散期待値による外挿法– 角運動量射影– 密行列積を利用したアルゴリズム実装• 研究成果– 軽い核の第一原理計算と炭素１２のホイル状態• クラスター状態の解析– ネオジム同位体の”形の相転移”• ニュートリノレス二重ベータ崩壊核行列要素
3

研究目的

豊かな核構造の記述、予言クラスター構造、形の共存、形の相転移、魔法数の消滅と新魔法数、不安定核の存在限界、…
12C ホイル状態

加速器施設による実験研究

新物理探索実験に必要な核データの創出
ニュートリノレス二重ベータ崩壊、シッフモーメント、 …

理研・RIビームファクトリー

核子間相互作用大規模計算による核構造研究

元素合成過程の理解に必要な核データの創出
中性子過剰核、核準位密度、ガンマ線強度関数、 …

4

原子核の存在限界（中性子ドリップライン）とMg同位体の「反転の島」と不安定核の存在限界@ Oakforest-PACS, 京計算機

N. Tsunoda, T. Otsuka, K. Takayanagi, N. Shimizu, T. Suzuki, Y. Utsuno,
S. Yoshida, and H. Ueno, Nature 587, 66 (2020)5軽い核の第一原理CI計算 (no-core MCSM)• 調和振動子ポテンシャルの固有状態で空間を展開する。
• 主殻の数（Nshell）を設定し、その中で許されるあらゆる配位を考える。

| Ψۧ =

𝑣1 |𝑚1 ۧ+𝑣2 |𝑚2 ۧ

Nshell=7
富岳＋
MCSMでNshell=7の計算が可能に+

𝑣3 |𝑚3 ۧ + ・・・

6

原子核における第一原理CI計算• 原子核の波動関数は、多数の「配位」の線形結合で表される。

|Ψۧ =𝑣1 |𝑚1 ۧ+𝑣2 |𝑚2 ۧ+𝑣3 |𝑚3 ۧ + ・・・

シュレディンガー方程式は、ハミルトニアン行列の固有値問題に帰着する

H  E m H m' v  Ev
m'mm'

許される基底(|𝑚𝑖 ۧ)の数がハミルトニアン行列の次元に相当する。
大次元(≲ 1011 )疎行列の固有値・固有ベクトルをランチョス法により求める。
量子化学における FCI (Full Configuration Interaction)に相当。
多様な量子多体問題（強相関電子系、分数量子Hall系など）に類似な手法。
7

核図表：原子核理論模型の適用領域
第一原理CI計算

原子核殻模型計算は、適用領域を中重核から大きく広げつつある。

中重核のCI計算（殻模型計算）
密度汎関数

模型空間を広げて、原子核の存在限界や、エキゾチックな殻構造の理解、高励起状態の統計的性質、さらには応用に向けた研究を進めている。
第一原理CI計算Taken from UNEDF project, 2009

中重核のCI計算

8ハミルトニアン行列： 大次元疎行列
殻模型ハミルトニアン行列24Mg (16O を芯とする) 28,503 dim.

KSHELL+京1012 次元の疎行列：
行列要素はオンザフライ生成
１つの状態は1012 次元のベクトルで表現⇒ 8TB/状態
N. Shimizu, T. Mizusaki, Y. Utsuno, and Y. Tsunoda, Comput. Phys. Comm. 244, 372 (2019). 9軽い核のCI計算とハミルトニアン行列次元• 対称性を仮定しなければ、対角化すべき行列の次元は𝑁𝑠𝑝 (𝜋)𝐶𝑍 × 𝑁𝑠𝑝 (𝜈)𝐶𝑁– 模型空間と粒子数に依存して、組み合わせ爆発で増えていく。
2820 4070 112 168 240 330

𝑁𝑠𝑝

直接対角化計算の限界

T. Abe, P. Maris, T. Otsuka, N. Shimizu, Y. Utsuno, and J. P. Vary, Phys. Rev. C 104, 054315 (2021)10

モンテカルロ殻模型・準粒子真空殻模型• CI計算：多数の単純な多体基底の線形結合で波動関数を表現。
• より効率的な少数の基底の線形結合で表現

モンテカルロ殻模型(MCSM)

角運動量射影（回転対称性の回復）
パリティ射影Nb

Ii 1

K  II   f i  g K PMK
P  ( D (i ) )複数基底の重ね合わせ

準粒子真空殻模型(QVSM)
Nb

Ii 1

K  I

核子のSlater行列式
Hartree-Fock計算に使われる形

陽子数・中性子数射影I   f i  g K PMK
P P Z P N  i
準粒子真空基底11Hartree-Fock-Bogoliubov 計算に使われる形

モンテカルロ殻模型(MCSM)の概念図• 効率的な基底を選択することによって、少数の基底によって固有状態を精度良く表現する。

巨大
しばしば1000億以上

単純な基底による表現

１００程度「最適化された」基底による表現
Review : N. Shimizu, T. Abe, M. Honma et al., Phys. Scr. 92, 063001 (2017)12

「近似の良さ」の指標としてのエネルギー分散期待値• シュレディンガー方程式を固有値問題として解く。
𝐻𝑣 = 𝐸𝑣 (𝑣 𝑇 𝑣 = 1)

• エネルギー期待値
𝐸 = 〈𝐻〉 = 𝑣 𝑇 𝐻𝑣
• エネルギー分散期待値〈Δ𝐻2 〉 = 𝐻2 − 𝐻 2 = 𝑣 𝑇 𝐻2 𝑣 − 𝐸 2• もし 𝑣が厳密に𝐻の固有ベクトルなら、𝐻2 𝑣 = 𝐸 2 𝑣より、 Δ𝐻2 = 0• 𝑣が固有ベクトルからわずかにずれていると、 Δ𝐻2 は小さな値をとる。エネルギー分散期待値は「近似の良さ」の指標となる。
13MCSMにおけるエネルギー分散期待値による外挿法
変分計算において、基底数を増やしていくと、エネルギー期待値がゆるやかに下がっていく。エネルギー分散期待値 Δ𝐻2 = 𝐻2 − 𝐻 2 を求める
ことにより、収束の判定や外挿による精密な物理量の推定を可能とした。
12C, N

2次多項式によるフィット＋外挿

shell=7, Daejeon16 int.

エネルギー(MeV)

1019 次元の行列の固有値問題に相当
0+3𝐽𝜋 = 0+30+2𝐽𝜋 = 0+2𝐽𝜋 = 01+

MCSM基底数
1019 基底目への外挿は困難

01+

エネルギー分散 (MeV2)
Δ𝐻2 = 𝐻2 − 𝐻 2
厳密解では Δ𝐻2 = 0

Δ𝐻2 → 0 への外挿が容易
同様な手法は強相関電子系における経路積分繰り込み群法、変分モンテカルロ法などに使われている14

量子数射影法
NbIi 1K  II   f i  g K PMKP  ( D (i ) )• 角運動量射影：– 「変形」して回転対称性を破った状態(スレーター行列式)を用意する。
– あらゆるオイラー角(実際には30,000点程度)に回転させ、重ね合わせることにより、角運動量が良い量子数となる空間へ射影する。

• 「自発的対称性の破れ」と、破れた対称性の回復。
• 核構造計算において、古くから使われてきた。
– 強相関電子系や量子化学の計算にも、スピン射影が用いられる。
e.g. T. Mizusaki and M. Imada, Phys. Rev. B 69 125111 (2004).
15

計算のボトルネック：
二体相互作用行列要素の計算
𝛼 𝛼𝛼
𝜓𝑎 𝑉𝑃 𝐽 𝜓𝑏 = ෍ 𝑊𝛼 ෍ 𝑣𝑖𝑗𝑘𝑙 (𝜌𝑘𝑖𝜌𝑙𝑗 − 𝜌𝑙𝑖𝛼 𝜌𝑘𝑗)𝛼

初期の疑似コード

𝑖<𝑗,𝑘<𝑙𝑣𝑖𝑗𝑘𝑙 : ハミルトニアンの対称性を反映した疎行列

do a = 1, np (≈O(104))
𝛼: 角運動量射影の分点rho(:,:) = make_rho_p(a)do n = 1, size(idx,2) (≈O(104))i = idx(1,n)
j = idx(2,n)
k = idx(3,n)
l = idx(4,n)
v = v_ijkl(n)x = x + w(a) * v * ( rho(k,i) * rho(l,j) – rho(l,i)*rho(k,j) )end doend do16モンテカルロ殻模型コードのアルゴリズム改良
ボトルネック：二体相互作用の期待値
𝛼 𝛼𝛼
𝜓𝑎 𝑉𝑃 𝐽 𝜓𝑏 = ෍ 𝑊𝛼 ෍ 𝑣𝑖𝑗𝑘𝑙 (𝜌𝑘𝑖𝜌𝑙𝑗 − 𝜌𝑙𝑖𝛼 𝜌𝑘𝑗)

𝛼𝑖<𝑗,𝑘<𝑙1𝛼𝛼𝜓𝑎 𝑉𝑃 𝜓𝑏 = ෍ 𝑊𝛼 ෍ 𝜌(𝑘,𝑖)𝑣 𝑘,𝑖 ,(𝑙,𝑗) 𝜌(𝑙,𝑗)2𝐽𝛼(𝑘,𝑖)ボトルネック部分の実行性能の測定値＠富岳 single core

期待値計算を、密行列積(DGEMM)で表現することにより、実効性能を大きく向上。
富岳にて、理論性能比で実行性能 15%程度、1000ノードまでの良好なストロングスケーリング並列性能を達成。
17
Y. Utsuno, N. Shimizu et al., Comput. Phys. Comm. 184, 102 (2013)

原子核構造計算とモンテカルロ殻模型、その発展ストロングスケーリング＠富岳
モンテカルロ殻模型, 4He, Nshell=7対角化すべきハミルトニアン行列の次元ストロングスケーリング＠富岳
準粒子真空殻模型, 68Ni, pfg9d5 space

提案手法による研究成果

富岳において、モンテカルロ殻模型は1,000ノードまで準粒子真空殻模型では15,000ノードまでの良好な並列性能を達成。

18

京・富岳成果創出加速プログラムによる成果 Ⅰ
軽い核の第一原理計算とクラスター構造
12C ホイル状態：

元素合成過程に重要なトリプルアルファ反応に重要な役割を果たす。
第一原理CI計算

実験室系の密度分布(8Be 基底状態)

物体固定座標系の密度分布
𝐽

角運動量射影 𝑃𝑀𝐾

モンテカルロ殻模型の波動関数
角運動量射影（回転対称性の回復）
パリティ射影NbJi 1K  JJ   f i  g K PMKP  ( D (i ) )複数基底の重ね合わせ

Slater行列式

重ね合わせにより、「形の揺らぎ」を記述する。
Wave function = f1+ f2+ f3

強相関電子系における Path Integral Renormalization Group method,
Resonating Hartree-Fock などに類似。

+ ......
20

京・富岳成果創出加速プログラムによる成果 Ⅱ
芯を仮定したCI計算

Nd同位体の核構造と「形の相転移」
Nd同位体の励起エネルギー、E2遷移確率、四重極能率: 実験値-: 計算値●模型空間(112Zr core)

中性子数150Nd21

ニュートリノレス二重ベータ崩壊の核行列要素
通常の二重ベータ崩壊（2νββ） ニュートリノレス二重ベータ崩壊探索実験が限られた特別な原子核は、２つのニュートリノと２つの電子を放出して崩壊する。

ニュートリノレス
二重ベータ崩壊（0νββ)
ニュートリノを放出しない特別なモード存在は未確認。

世界中で行われているが、いまだ発見されていない。半減期が計測されたとすると、その値からニュートリノ有効質量を求めるには、原子核構造計算によって核行列要素と呼ばれる物理量を精密に求める必要がある。
半減期:

位相因子

核行列要素
(NME)

ニュートリノ
有効質量ニュートリノ粒子がマヨラナ粒子であるかを明らかにし、物質優勢宇宙の解明のカギとなる探索実験の対象となっている原子核:
48Ca, 76Ge, 82Se, 96Zr, 100Mo, 116Cd, 124Sn, 130Te,136Xe, 150Ndなど
22

二重ベータ崩壊の核行列要素：様々な理論計算
EDF: 密度汎関数法
QRPA: 準粒子乱雑位相近似

Half life(exp.)

核行列要素

effective
neutrino mass

IBM: 相互作用ボゾン模型

SM:殻模型計算

本研究の計算結果

A. Gando et al., Phys. Rev. Lett. 117, 082503 (2016)

48Ca, 76Ge, 82Se, 96Zr, 100Mo, 116Cd, 124Sn, 130Te, 136Xe, 150Nd

Engel & Menéndez, Rep. Prog. Phys. 80, 046301 (2017)より改変
23Summary• 原子核のCI計算(殻模型計算)は、陽子・中性子を構成要素とする原子核を量子多体問題として解く最も有力な数値計算手法の一つ。
• 1千億次元を超えるような巨大な疎行列の対角化が必要となる。この困難を回避するためにモンテカルロ殻模型（MCSM)、準粒子真空殻模型（QVSM)を開発、密行列積を用いるように計算アルゴリ
ズムを改良し、京・富岳を用いて研究を進めている。
– エネルギー分散期待値による外挿法– Lebedev quadrature の導入による角運動量射影の計算量削減– 密行列積利用によるアルゴリズムの改善• 従来困難であった、軽い核の第一原理CI計算によるクラスター構造の記述や、芯を仮定したCI計算による質量数１５０程度の重い核の系統的な研究が可能となりつつある。

24