1E11 二核遷移金属錯体のジラジカル因子と三次非線形光学物性に対するアキシアル型配位子の効果についての理論的研究(阪大院基礎工)○山田大志、高椋章太、福田幸太郎、廣崎裕多、中野雅由

Theoretical study of axial ligand effect on the diradical characters andthird-order nonlinear optical properties of dinuclear transition-metal compounds(Graduated School of Engineering Science, Osaka University)○Taishi Yamada, Shota
Takamuku, Kotaro Fukuda, Yuta Hirosaki, Masayoshi Nakano【序】超大容量三次元光ディスク、超高速光スイッチングなど将来の光デバイスの基盤材料として、近年、
非線形光学(NLO)物質の分子設計に注目が集まっている。我々は三次NLO物質としてこれまで殆ど検討されていなかった開殻性を示す分子に着目し、量子化学に基づき、開殻性の指標であるジラジカル因子(y)と分子レベルの三次非線形光学特性を表す第二超分極率γとの相関を検討してきた。その結果、「yが中程度の開殻分子系は、サイズや構造が類似の閉殻系や完全開殻系に比べて著しく大きいγを示す」というy-γ相関を見出した[1]。最近、実在の有機開殻分子系において巨大な二光子吸収断面積が実測され、このy-γ相関の妥当性が示された[2]。
	  一方、金属−金属間に多重結合を有する多核遷移金属錯体では、弱い d−d 軌道相互作用に起因する各 dX
軌道[X =σ, π, δ]の結合が開殻性を有しており、開殻 NLO 系の新たなターゲットとして有望であると期待される。以前、我々は金属−金属結合の一重項開殻性が三次 NLO 物性に及ぼす効果を明らかにするため、配位子なしの同核(homo)及び異核(hetero)二核系における y-γ相関を量子化学計算により検討した。その結果、両系ともに有機π電子系とは異なり、中間ジラジカル性を持つ dσ電子がその大きな三次 NLO 物性の起源であるという構造-物性相関を明らかにした[3, 4]。しかし、実在の金属錯体は配位子による安定化の効果を受けて存在しており、その金属−配位子間の電荷移動効果や軌道相互作用が三次 NLO 物性に大きく影響することが予想される。金属二核系における配位子は、結合軸方向に対して垂直に配位するエクアトリアル型配位子と、結合軸方向に沿って外側から配位するアキシアル型配位子に大別され(図 1)、
エクアトリアル型配位子については、既に以前の研究で y とγに及ぼす効果を明らかにした[5]。一方、未検討のアキシアル型配位子は、三次 NLO 物性に主寄与である dσ軌道と大きな軌道相互作用を持つと考えられ、金属−金属結合の開殻性に起因する y やγにより大きな影響を与えると予想される。そこで本研究ではアキシアル型配位子を含む二核遷移金属錯体モデルを構築し、配アキシアル

エクアトリアル

図 1. アキシアル、エクアトリアル	  	
  型配位子を含む二核金属錯体

位子なしの系と比較することによりアキシアル型配位子が y やγに与える効果を検討する。
【モデル系・計算手法】本研究では１つの dσ結合と、２つのdπ結合、２つの dδ結合を持ち形式的に５重結合を形成する
Mo(I)2 系を参照系とし、この系にアキシアル位に CO が２つ配位

図 2. Mo(I)2 と Mo(I)2(CO)2 の構造

した Mo(I)2(CO)2 錯体モデルを検討する(図 2)。Mo(I)2(CO)2 の構造最適化は B3LYP 法を用い、C−O、Mo−C間の距離はそれぞれ 1.113 Å、2.448 Å を得た。また配位子の効果をより強くしたモデルとして Mo−C 間の距
離を 2.0 Å とした系も検討した。ジラジカル因子 y は UHF 計算による自然軌道のスピン射影占有数から求め、
静的γの結合軸方向成分は UCCSD 法で求めた外部電場存在化でのエネルギーから有限場法を用いて算出した。以上の計算では、Mo には内核に Stuttgart/Dresden の ECP を適用した基底関数 SDD を使用し、C, O には
6-311++G**を用いた。実在の二核錯体では配位子の種類により核間距離が変化し、それに対応し y が変化すると考えられるため、本研究では核間距離(R)のみ変化させ、それ以外の構造は固定して y-γ相関を検討する。
【 結 果 】 図3、図4にそれぞれMo(I)2とMo(I)2(CO)2におけるy(dX)のR依存性を示している（Mo−C間は2.0Å）。
図3よりMo(I)2はRが2.5 Å以上でジラジカル因子の大きさの順が、y(dδ) > y(dπ) > y(dσ)となることが明らかになった。Rが2.4 Å付近で見られるy(dσ)増大はdσ結合を構成するdz2軌道の逆位相のリングの影響であると考察される。一方、図4より、配位子を導入したMo(I)2(CO)2は全てのRでジラジカル因子の大きさの順がy(dδ) >y(dσ) > y(dπ)となっており、y(dσ)とy(dπ)の大きさの順が配位子なしのものと逆転することが分かった。また、同じRにおいて両系で各y(dX)の結果を比較すると、y(dπ)とy(dδ)は1.0

ほぼ同じ値をとるのに対し、y(dσ)は配位子を導入することで値が0.8

Mo(I)2(CO)2のy(dX)のR依存性の結果については、Mo−C間が2.0Åの場合よりもy(dσ)の増大は小さく、y(dσ)とy(dπ)の大きさの順の逆転

y(dX) [-]

増大することが明らかになった。またMo−C間が2.448 ÅのはRが2.5Å以上でしか見られなかった。これはMo−C間の距離が小

0.6
0.4
0.2

さいほど配位子によるy(dσ)の増大効果が大きいことを示唆していy (dσ)y (dπ)y (dδ)

る。このリガンドによるy(dσ)の増大効果の原因を考察するため、

0.0
1.5 2.0 2.5 3.0 3.5 4.0 4.5
R [Å]Mo−C間が2.0Åの場合について、R=2.0 ÅにおいてUHF計算で得ら

図 3. Mo(I)2 における y の R 依存性

1.0

れたdσ型の結合性自然結合の分布を比較した（図5）。その結果、

0.8

方向の自然軌道の成分（黄色）の寄与が減少する一方、dz2軌道のリングに由来する逆位相の成分（青色）の寄与が増大していることが明らかになった。以上の結果より、アキシアル配位子導入に

y(dX) [-]

配位子を導入した系では、金属−金属結合の強さと相関する結合軸

0.6
0.4

y (dσ)y (dπ)y (dδ)

よるy(dσ)の増大は、dσ軌道の結合軸方向の成分の寄与の減少に起0.2因すると考察される。当日はy(dσ)軌道の成分の変化の詳細な考察

0.0
1.5 2.0 2.5 3.0 3.5 4.0 4.5
R [Å]に加え、γの結果についても合わせて報告する。
【参考文献】[1] M. Nakano et al., J. Phys. Chem. A 109, 885 (2005);
Phys. Rev. Lett. 99, 033001 (2007). [2] K. Kamada et al. Angew. Chem.
Int Ed. 46, 3544 (2007). [3] H. Fukui et al., J. Phys. Chem. Lett. 2, 2063

図 4. Mo(I)2(CO) 2 における y の R 依存性

Mo(I)2
Mo(I)2(CO)2(2011); J. Phys. Chem. A 116, 5501 (2012). [4] T. Yamada et al., Chem.
Phys. Lett. 579, 73 (2013). [5] Y. Inoue et al., Chem. Phys. Lett. 570, 75(2013).

図 5. R=2.0 Å における Mo(I)2 と Mo(I)2 (CO)2（Mo−C 間は 2.0Å）の結合性 dσ NO 図（黄色：結合軸成分、青色：逆位相のリング）

1E-12開殻因子の制御に基づく新規なシングレットフィッション分子の提案(阪大院基礎工) ◯伊藤聡一、中野雅由

Proposal of new candidate molecules for singlet fission
based on the tuning of diradical character(Graduate School of Engineering Science, Osaka University) ◯Soichi Ito, Masayoshi Nakano【序】シングレットフィッション(singlet fission: SF)は、多重励起子生成の一種であり、太陽電池の光
電変換効率を劇的に向上させる可能性があるとして近年注目を集めている[1, 2]。これは主に分子集合体において起きる現象であり、その場合、初期過程として、光を受けて一重項励起状態となった分子(一重項励起子)が、近くの分子との相互作用を経て二つの分子上に三重項励起子をそれぞれ一つずつ、合計で二つ生成する。古くは 1960 年代にアントラセン分子結晶での遅延蛍光の測定で発見され、2009
年に Hanna ら[2]が太陽電池への応用の可能性を示唆したことにより研究が盛んに行われるようになった。多くの研究グループがそのメカニズムについて実験・理論両面から研究を進めているが、現在のところ、その対象分子は極めて少なく、SF に適した分子の理論設計指針の構築と具体的な分子の設計が喫緊の課題である。SF が効率的に生じる条件は完全には明らかになっていないものの、最低限の条件として、以下のエネルギー準位に関する適合条件が提案されている[1]。
2E(T1) – E(S1) ~ 0 or < 0(1),

2E(T1) – E(T2) < 0(2)

ここで、E(S1)は一重項最低励起エネルギー、E(T1)及び E(T2)は、各々、三重項の第一、第二励起エネ
ルギーである。式(1)は一つの一重項励起子から二つの三重項励起子を生成するための条件、式(2)は生成した三重項励起子が再結合してより高いエネルギー準位の三重項励起子になることを防ぐための条件である。一般に、多くの有機分子においては E(T1)が高く、この式(1)が満たされないことが研究対象の分子種を狭めている原因の一つである。我々は、最近の研究により、条件式(1)、(2)が開殻因子（化学結合の弱さの指標）yi[3]と強い相関があることを明らかにした[4-6]。すなわち、条件式(1)、(2)を満たす系は、比較的小さいか中程度の y0 (> 0)と、非常に小さい y1 (<< y0)を持つことを見出した。これは、
y0 が現れることにより、E(T1)の急激な低下が起こること、一方、y1 が極めて小さいことにより E(T2)の低下が抑制されることから説明される。
	  一方、太陽電池への応用の観点からは、(i) 式(1)ができるだけ~0 に近いこと、及び、(ii) SF で生じた三重項励起子が界面で電荷分離を効率的に起こすこと、の二点が共に重要である。条件(i)は比較的小さな y0 の系を探索することで満たされるが、条件(ii)を満たすには、電荷分離の駆動力となる E(T1)の大きさを、式(1)、(2)を満たした上で、できるだけ大きくする必要がある。最近の我々の研究で、条件(ii)を満たすためには、上述の開殻因子の制御に加えて、最高占有軌道(HOMO)-最低非占有軌道(LUMO)
間の空間分布の重なりが大きいことが重要であることを明らかにした[4]。
	  本研究では、これらの結果を基に、開殻因子 yi および HOMO−LUMO 空間分布の二つの因子に注目
して、条件(1)、(2)及び(i)、 (ii)を満たす新規な SF 分子を設計する。図 1 に本研究で対象とする分子を示す。これらの分子は双性イオン−ジラジカルの二つの共鳴構造を考えることができるため、開殻性を持つことが予想される。また、開殻因子はπ共役長にも左右されるので、分子サイズ n = 3、4 の二種類を検討した。また、開殻因子に影響を与える共鳴構造の寄与の違いを調べるため、芳香環への窒素原子導入の効果も検討した。これら二種類の化学修飾により開殻因子を調節し、SF 分子の設計を行った。
【計算方法】分子構造の最適化を spin-flip 時間依存密度汎関数法(spin-flip TDDFT)で行い、各励起エネルギーを extended general multi-configuration quasi-degenerate perturbation theory (X-GMCQDPT)法により算出した。電荷分布及び芳香族性の算出にはスピン制限/非制限の長距離補正
BLYP(LC-RBLYP/LC-UBLYP)汎関数を用いた。芳香族性は nucleus independent chemical shift (NICS) (1)により評価した。開殻因子の算出には近似スピン射影スピン非制限 Hartree-Fock(PUHF)法を用いた。基
底関数は、構造最適化には 6-311G*を、それ以外にはすべて 6-31G*を用いた。
【結果】各モデル分子の開殻因子 yi と各エネルギーの結果を表 1 に示す。まず開殻因子を分子の種類で比較すると c、b、a の順に、分子サイズで比較すると n = 3、4 の順に、y0 と y1 が共に増大する。前者は芳香環への N 原子の導入数が増えるほど HOMO–LUMO のエネルギーギャップεH-L が大きくなる
ことにより、また後者は π 共役の増大で同ギャップが小さくなることにより理解できる。一方、開殻因子の上昇と同時に E(S1)、E(T2)、E(S1)のいずれもが小さくなった。特に y0 の増大に対して E(T1)が急激に減少した。SF のエネルギー適合条件(1)を調べると、a-3、a-4、b-4、c-4 において満たされることがわかった。ただし、a-4、b-4 については、y0 が比較的大きく、E(S1)に比べ E(T1)が低下し過ぎるため、
一重項励起子から三重項励起子へ分裂する際のエネルギーロスが大きい[4-6]。一方、a-3 及び c-4 は、
y0 が比較的小さく、このエネルギーロスが小さくなるため、SF として望ましい性質を持っている。さらに、三重項励起子の再結合を抑制するための条件(2)についてもこの二つの分子は満たしている。三重項励起エネルギーと HOMO−LUMO 空間分布との

分子	
 

y0

y1

εH–L

E(S1)

E(T1)

E(T2)

関係等、その他の結果については当日報告する。

a-3

0.32

0.00

7.2

2.32

1.16

2.37

b-3

0.19

0.00

7.8

2.51

1.40

2.43

c-3

0.15

0.00

8.2

2.66

1.65

2.69

N

H
N

•

N

X2
X1

H
N

N

•

N
n-2

N

H
N

X2

H
N

N

X1

N

a-4

0.63

0.03

5.7

1.55

0.30

1.53

n-2

b-4

0.53

0.02

6.4

2.17

0.77

1.90

c-4

0.27

0.01

7.1

2.25

1.09

2.54

N

図 1. モデル分子の双性イオン−ジラジカルの共鳴構造(a-n: X1=X2=CH, b-n: X1=N, X2=CH, c-n: X1=X2=N)【参考文献】[1] M. B. Smith and J. Michl, Chem. Rev. 110, 6891 (2010).
[2] M. C. Hanna and A. J. Nojik, J. App. Phys. 100, 074510(2006).

[3] K. Yamaguchi, Self–Consistent Field. Theory

表 1. 図 1 モデル分子の開殻因子 y0 [-]、
HOMO–LUMO エネルギーギャップ[eV]、各一
重項、三重項状態の励起エネルギー [eV]and Applications; Carbo, R., Klobukowski, M., Eds.;
Elsevier: Amsterdam, (1990).

[4] T. Minami and M. Nakano, J. Phys. Chem. Lett. 3, 145 (2012).

Ito and M. Nakano, J. Phys. Chem. Lett. 4, 2133 (2013).
19729 (2012).

[5] T. Minami, S.

[6] S. Ito, T. Minami and M. Nakano, J. Phys. Chem. C. 116,

1E-13高周期典型元素を含むパラ及びメタキノジメタン誘導体の開殻性と	
 第二超分極率に関する理論的研究	
 (1 阪大院基礎工、2 筑波大数理物質系) ○福田	
 幸太郎 1、廣崎	 裕多 1、松井	 啓史 1、	 岸	 亮平 1、野澤	
 竹志 2、一戸	 雅聡 2、関口	
 章 2、中野	 雅由 1	 	
 Theoretical investigation of the open-shell nature and the second
hyperpolarizabilities of para/meta-quinodimethane derivatives involving maingroup elements(1Graduate School of Engineering Science, Osaka University, 2Graduate School of Pure andApplied Sciences, University of Tsukuba) ○Kotaro Fukuda1, Yuta Hirosaki1, Hiroshi Matsui1,
Ryohei Kishi1, Nozawa Takeshi2, Masaaki Ichinohe2, Akira Sekiguchi2, Masayoshi Nakano1【序】非線形光学（NLO）現象は外部電場に対する応答として生じる分極の非線形項に起因する現象で、三次元メモリなど将来の光エレクトロニクスへの応用が期待されている。我々は、量子化学に基づき、世界に先駆けて、開殻性と NLO 物性との間に強い相関関係が存在す
ること、そして特に中間的な開殻性を持つ系で NLO 物性が増大すること（y − γ 相関、ここ
で第二超分極率γは分子レベルの三次 NLO 物性を示す）を明らかにし[1]、一重項開殻性を持
つ分子系が優れた NLO 物質となることを提唱した。開殻性は、量子化学計算によって算出されるジラジカル因子 y（0 ≤ y ≤ 1; y = 0:完全閉殻、y = 1:完全開殻）により定量的に表すことができ、分子が有する幾何構造やサイズ、構成元素の違いに応じて大きく変動する。この y – γ
相関を検討する際に用いた系の一つに一重項基底状態を持つパラキノジメタン（PQM）がある。PQM はキノイド構造（閉殻）とベンゼノイド構造（開殻）の二つの共鳴構造を描くことが可能であり（Fig. 1）、原理的には基底状態において中間的な開殻性を示す可能性がある。しか

H

H
C

H

H

C

H
C

H

C

H

H

Figure 1. Resonance structures of PQMしながら、量子化学計算の結果、平衡構造では小さなジラジカル因子しか示さず、ほぼ閉殻の分子であることがわかった。また、環外 C-C 結合長を変化させた仮想 PQM モデルでは、
1) C-C 結合伸長に伴い環外炭素の位置に奇電子（不対電子）が出現し、一重項開殻性が生じること、2)中間開殻性を示す長い環外 C-C 結合の場合にγが極大値をとること、が量子化学計
算により判明した[1]。以上の結果から、NLO 物性が増大する中間的な開殻性を示す長い環外C-C 結合を持つ PQM は非常に不安定であり、この分子骨格をもつ純粋炭化水素では中間開殻性を実現することは困難であることがわかった。
一方、高周期典型元素は第二周期の炭素原子に比べて混成が生じにくいため、通常の炭素−炭素結合よりも弱い結合を形成する傾向がある[2]。すなわち、PQM 骨格の一部を高周期典型元素で置き換えることにより、平衡構造で大きな開殻性を示す PQM
類似系の実現が期待される。実際に、最近、PQM 及びメタキノジメzH

タン（MQM）骨格の環外炭素をシランで置き換えた PQM、MQM 類
似系が合成され、X 線結晶構造解析から、こうした類似系は通常のPQM や MQM が持つ環外 C-C 結合に比べて遥かに長い環外 Si-C 結合を持つことが明らかになってきた[3]。そこで本研究では、PQM、

H
X

X

H

H

1a. para-quinodimethane(PQM, X = C)
2a. PQM-like compound (X = Si)
3a. PQM-like compound (X = Ge)

MQM の環外炭素原子を同族の高周期典型元素である Si、Ge で置き換えた PQM、MQM 類似体に着目し、その構造、開殻性 y 及び静的
第二超分極率 γ の相関関係の解明を行い、高周期典型元素導入の効

H

H

H

X

X

H

果について検討する。
時間依存密度汎関数法を用いて行った。一重項開殻性は APUHF スキ

1b. meta-quinodimethane(MQM, X = C)
2b. MQM-like compound (X = Si)
3b. MQM-like compound (X = Ge)ームによるジラジカル因子 y により評価し、長軸（z 軸）方向成分のFigure 2. Model structures【理論計算】系の構造最適化は BHandHLYP 汎関数を用いた spin-flip

静的第二超分極率( γzzzz)は UCCSD(T)を用いて有限場（Finite Field）法により算出した。基底関数としては、環外 X 原子に LANL2DZ(d,p)、
その他の原子に 6-311+G**を用いた。以上全ての計算には Gaussian 09 を用いた。
【結果】各 PQM 類似体の一重項基底状態における構造を比較すると、環外 X が炭素からより高周期の元素となるにつれて環外 C-X 結合が伸長（1a : 1.349 Å, 2a: 1.823 Å, 3a : 1.917 Å）
し、中央ベンゼン環の結合交代（BLA）が小さく（1a : 0.109 Å, 2a : 0.034 Å, 3a: 0.023 Å）なっており、すなわちキノイド構造の寄与が小さく、ジラジカル構造の寄与が大きくなることがわかった。こうした平衡構造の変化は一重項開殻性の発現を示しており、実際に 2a、3a では中間的な開殻性(2a : y = 0.684, 3a: y = 0.747)を持つ。次にγzzzz を比較したところ、ほぼ閉殻の1a（γzzzz = 24800 a.u.）に比べ、中間的な開殻性を持つ 2a（γzzzz = 185000 a.u.）、3a（γzzzz = 175000
a.u.）では約 7 倍のγzzzz 値を示し、y − γ 相関に従う NLO 物性の増大が見られた。これらの結果は実在分子である 2a、3a が優れた NLO 材料となることを示唆するだけでなく、分子骨格への高周期典型元素導入による開殻性及び NLO 物性の制御を可能にする新しい設計原理の開拓に繋がると期待される。MQM 類似体に関する結果及び電荷移動を含めたより詳細な解析は当日報告する。
Table: External C-X bond length (R), y and γzzzz for singlet ground state of each PQM derivative

X

R [Å]

BLA [Å]y [-]γzzzz [a.u.]

C (1a)

1.349

0.109

0.143

24800

Si (2a)

1.823

0.034

0.684

185000Ge (3a)

1.917

0.023

0.747

175000【参考文献】	
 
[1] M. Nakano et al., J. Phys. Chem. A 109, 885 (2005); Phys. Rev. Lett. 99, 033001 (2007); J. Chem.
Phys. 133, 154302 (2010); J. Chem. Phys. 138, 244306 (2013).
[2] F. Breher, Coord. Chem. Rev. 251, 1007 (2007).
[3] T. Nozawa, A. Sekiguchi et al., J. Am. Chem. Soc. 133, 5773 (2011).

1E14
量子散逸系ダイナミクスの確率的リウビル-フォン・ノイマン方程式
を用いた解法：スペクトル密度形状への依存性(東北大院理) ○今井甫、大槻幸義、河野裕彦
A method to calculate the dynamics of quantum dissipative systems
with a stochastic Liouville-von Neumann equation:
Dependence of the dynamics on the shape of its spectral density(Tohoku Univ.) ○Hajime Imai, Yukiyoshi Ohtsuki, Hirohiko Kono
[序]大きな系において、その全系を興味ある系とその環境に分けて考えることで、環境から影響を受けた興味ある系の動的な過程について考えることができる。このときの興味ある系が量子散逸系と呼ばれ、物理化学や、固体物理、生物物理、量子測定などの分野の発展において避けることのできない困難で重要な問題を含んでいる。このような量子散逸系の研究では系と環境の間の複
雑な相互作用をより簡単な形で記述することに興味が持たれており、さまざまな近似に基づく定式化や量子統計力学に基づいた定式化が試みられている。
例えば確率的リウビル-フォン・ノイマン方程式に基づいて、量子散逸系に作用する環境の影響を確率過程で表された有色ノイズとして定式化することができる[1]。この手法は揺動を表す確率過程を含むハミルトニアンに基づく確率的リウビル方程式[2]に散逸項を表すもう一つの確率過程を加えた、揺動散逸定理を満たす形式となっており、有限の温度での時間発展を表すことができる。
[理論]全系のハミルトニアンを Caldeira-Leggett モデルで表す：

ここで、は系の座標、運動量、質量、ポテンシャルを、は調和振動子の集まりで表された環境のモード の座標、運動量、質量、周波数を、また、 は系と環境のモードのカップリング強度をそれぞれ表す。
このとき、系の縮約密度行列
きることが示される[1]：が以下の確率過程を含む時間発展方程式によって計算で

ここで、は系が環境から受ける作用を表す確率過程の 番目のサンプルであり、は 2K 個の独立な白色ガウスノイズであり、関数と置ける。ただし

は系と環境の全モードのカップリング強度の情報を含むスペクトル密度関数：と温度によって定まる、影響汎関数の積分核の周波数成分である。環境が十分多数のモード からなるとき、スペクトル密度は連続関数で表される。

[応用]この手法の適用例として、光合成において光化学反応中心へエネルギーを伝達する役割をもつといわれている、FMO タンパク質複合体の中の色素間の励起エネルギー移動系を２状態でモデル化した系を考える。FMO タンパク質複合体のスペクトル密度は実験的に求められているが[3]、従来の励起エネルギー移動ダイナミクスの研究においては、解析的に扱いやすい Drude-Lorentz型を仮定していた（Fig. 1）。我々の手法はスペクトル密度関数として任意の関数型を用いることができ、とくに Drude-Lorentz 型スペクトル密度の場合、励起エネルギー分布の時間発展の計算
結果は階層型量子マスター方程式の結果[4]を再現した。また、さまざまなスペクトル密度関数に対する分布の時間発展を計算し（Fig. 2）、その形状を特徴づけるパラメータへの系統的な依存性を見出した。

Fig. 1.スペクトル密度関数の例.

Fig. 2. 励起子分布の時間変化のスペクトル密度関数依存性.

[1] J. T. Stockburger, Chem. Phys. 296, 159 (2004).
[2] R. Kubo, J. Math. Phys. 4, 174 (1963).
[3] M. Wendling, et al., J. Phys. Chem. B 104, 5825 (2000).
[4] A. Ishizaki & G. R. Fleming, J. Chem. Phys. 130, 234111 (2009).

1E15

Pederson–Khanna DFT 法による零磁場分裂テンソルのスピン軌道項計算の解析

12( 阪市大院理・ 阪大院基礎工) ○杉崎研司 1・豊田和男 1・佐藤和信 1・塩見大輔 1・北川勝浩 2・工位武治 1

Analysis of the Pederson–Khanna DFT method for the calculations of the spin–orbit
term of zero-field splitting tensors(1 Osaka City University, 2 Osaka University) ○Kenji Sugisaki,1 Kazuo Toyota,1 Kazunobu Sato,1 Daisuke
Shiomi,1 Masahiro Kitagawa2 and Takeji Takui1【序】 三重項以上のスピン多重度を持つ分子系の電子スピン構造を特徴づける物理量である零磁場分
裂テンソル(D テンソル)の第一原理計算が近年注目されている。非相対論的 Schrödinger 方程式から出発する摂動論では、スピン–スピン双極子項(DSS テンソル)とスピン軌道項(DSO テンソル)がそれぞれ摂動
一次および二次の項として現れる。我々は、DSO テンソルの高精度計算手法として、ハイブリッド
CASSCF/MRMP2 法を提案し、様々な分子系に適用してきた[1]。本手法により、有機高スピン分子のDSO テンソルの定量計算が可能となったが、単分子磁石など、大きな分子系の DSO テンソル計算には、
CASSCF 法などの多配置理論よりも DFT 法の方が計算量の観点から適していると考えられる。
DFT 法による DSO テンソル計算法は、Pederson–Khanna (PK)法[2]、quasi-restricted orbital (QRO)法[3]、coupled perturb (CP)法[4]の 3 手法が知られているが、このうち最も頻繁に用いられ、実験値とよい一致を示す傾向があるのは PK 法である。しかし、PK 法が高精度 ab initio 法と大きく異なる DSO テンソルを与える分子系も幾つか報告されている[5]。その原因について、DFT 法ではスピン反転を伴う励起の記述が不正確であると指摘している論文はあるが[6]、具体的にどの軌道からどの軌道への励起が DSO テンソルに大きく寄与するかなど、詳細な解析は行われてこなかった。今回我々は、軌道対の電子占有数に基づく
軌道領域分割による解析法を提案し、DSO テンソル計算における PK 法の振舞いについて議論する。
【Pederson–Khanna 法】 PK 法の解式は、以下の形である。
















1 occ vir  i h  a  a h  i
1 occ vir  i h  a  a h  i
D 



 a   i a   i4S 2 i a
4S 2 i a
SO

1 occ vir  i h  a  a h  i
1 occ vir  i h  a  a h  i
 2 
 2 
 a   i a   i
4S i a
4S i a(1)

ここで、h はスピン軌道ハミルトニアンの空間部分、φ は分子軌道、ε は軌道エネルギーである。(1)式により計算した DSO 値(D = DZZ − (DXX + DYY)/2)は、S2, Se2 などでは実験値の半分程度と過小評価するが、vanWüllen は、電子スピンの量子性を考慮し、Ŝx と Ŝz が交換しないことから、(1)式の各項の前係数に 1/S2 の代わりに 1/S(S − 1/2)を用いることで改善がされると報告している[7]。これまで、PK 法による DSO テンソルの解析は、(1)式の各項で部分和をとるものが主流であった。しかし、例えば SOMO から空軌道への励起を考えると、 → と → というスピン配置が可能であり、DSO テンソルには逆符号で寄与する。すなわち、
もしも空軌道が、軌道と軌道で空間分布と軌道エネルギーが一致していれば、SOMO から空軌道への励起の寄与は完全に打ち消しあう。従来のスピン配置に基づく解析は、このような本来打ち消しあいが起こる配置を分けて部分和を取っているため、各項から見掛け上大きな DSO テンソルの寄与が出てくる。例
として、Table 1 に S2 分子の 1 3g−基底状態の DSO テンソル解析結果(UBOP/cc-pVTZ level, van Wüllen

の提案した前係数を使用)を示す。なお、実験値は DSO = 23.6 cm−1 である。
【軌道対の占有数に基づく軌道領域分割解析法】 まず、Kohn–Sham 分子軌道において、軌道と
軌道で最も重なり積分の大きなペアを軌道対として定義し、軌道対の電子占有数が 2, 1, 0 の領域を、それぞれ二電子占有領域(doubly occupied region: DOR)、一電子占有領域(singly occupied region: SOR)、
非占有領域(unoccupied region: UOR)とし(Figure 1)、例えば SOR から UOR への励起（あるいは軌道対か
ら軌道対への励起）について部分和をとることで解析を行う。各領域間の励起で許容されるスピン配置はTable 2 にまとめた。Table 3 には、Table 1 と同じ S2 分子の DSO テンソルの、軌道領域分割による解析結果を示す。Table 1 で見られたような、見掛け上大きな DSO テンソルへの寄与は見られず、SOR から SOR への 励 起 が DSO テ ン ソ ル に 主 に 寄 与 し て い る こ と が 一 目 瞭 然 で あ り 、 こ れ は full valence
SOS-CASSCF/def2-TZVPP による計算結果[1]とも一致している。当日は、Figure 2 に示す、PK 法が ab
initio 法と異なる DSO テンソルを与えた系の DSO テンソル解析や、PK 法における交換相関汎関数依存性
および基底関数依存性についても議論する。
Table 1. S2 分子の従来法による DSO テンソル主値の解析結果
DSOXX/cm−1

DSOYY/cm−1

DSOZZ/cm−1

→

−30.42796−29.33822−20.04893

→−29.16052−28.10046−20.39088

→

33.93117

32.86437

47.74227

→

25.14973

24.06664

17.90904

Total−0.50758−0.5076725.21150

Table 2. 軌道領域間の励起で可能なスピン配置
Allowed excited spin configurations
DOR → UOR

 → 

 → 

SOR → UOR

 → 

→

DOR → SOR

 → 

SOR → SOR

→

 → 

→
→
Figure 1. 軌道の領域分割の概念図

Table 3. S2 分子の軌道領域分割による DSO テンソル主値解析結果
DSOXX/cm−1

DSOYY/cm−1

DSOZZ/cm−1

DOR → UOR

0.01057

0.010480.00030

SOR → UOR

−0.30853−0.30853−0.01249

DOR → SOR

−0.20962−0.20962−0.00566

SOR → SOR

0.00000

0.00000

25.22935

Figure 2. 採り上げる分子系[1] K. Sugisaki, et al. in EPR of Free Radicals in Solids I, 2nd Edition, pp. 363–392., and references therein.
[2] M. R. Pederson, S. N. Khanna, Phys. Rev. B 1999, 60, 9566–9572.
[3] F. Neese, J. Am. Chem. Soc. 2006, 128, 10213–10222.
[4] F. Neese, J. Chem. Phys. 2007, 127, 164112.
[5] 杉崎・豊田・佐藤・塩見・北川・工位 第 16 回理論化学討論会, 3L15 (2013).
[6] S. Ye, et al. Inorg. Chem. 2010, 49, 977–988.
[7] C. van Wüllen, J. Chem. Phys. 2009, 130, 194109.

1E17

反応経路網の全面探索：GRRM 法の超並列化（量子化学探索研究所, 東北大院理） ○大野公一 1,2、岸本直樹 2、岩本武明 2
Global reaction route mapping and development of its massively parallel control(Institute for Quantum Chemical Exploration1: Graduate School of Science, Tohoku University2)
Koichi OHNO1,2, Naoki Kishimoto2, Takeaki Iwamoto2【序】ポテンシャル表面の非調和下方歪(ADD)に着目することによって個々の化学組成の安定構
造・遷移状態・固有反応経路を全面自動探索する GRRM 法[1]について最初の報告を行ってから、
１０年が経過した。その間に計算機の機能は大幅に向上し、多数の演算装置を装備したマルチコ
ア・マルチノード環境が一般的となった。今回は、マルチコア環境における反応経路網全面探索
結果の特徴と GRRM 法を超並列環境で利用可能とするマルチノード化の試みについて報告する。

【方法】量子化学計算には Gaussian プログラム（g03, g09）を使用し、ADD に着目する超球
面探索を並列化して行う GRRM プログラム(GRRM11) [2]を用いて、反応経路網の全面探索を行った。GRRM11 そのものの並列処理は１つの計算ノード内に限定されているが、さらに多数のノードまたは計算機を併用して GRRM の並列処理を行うための制御プログラムとして
NeoGRRM を開発し[3]、多数のコアを利用して反応経路網の全面探索を行った。対象として、
水素、窒素、水、一酸化炭素、二酸化炭素、シアン化水素、アセチレン、アンモニア、メタンなどの基本的な分子を数個程度含む閉殻系を中心に約 100 種類の化学式を選定した。
【結果・考察】反応経路網全面探索は、未知の新反応・新構造を発見する可能性を秘めている。水性ガスシフト反応（H2O+CO→
H2+CO2）について、３０年以上にわたる長期間の研究により、ギ酸の２つの異性体を経る３段階の反応機構が知られていたが、１つの TS を経る直接機構はまったく知られていなかった。このほど、GRRM 法の適用により、１つの TSを経る直接機構が見出された[4]。
同じ化学式にどれほどの異性体数（EQ 数）があるか、興味が持たれるが、調査結果によると、図２に示すように、 図１ 水性ガスシフト反応の新機
原子数に対する相関は極めて悪い。例えば、原子数 6 の 構（左枠）と多段反応機構（矢印）H4CO や H4C2 の異性体は 2 種類のみだが、HC2NO2 の異性体数は 217 にもなる。また、原子数 7 の H5CNの異性体は 2 種類のみであるが、H3C2NO では 201
種類、H2C3O2 では 207 種類にもなり、同じ原子数でも、異性体数（EQ 数）は 100 倍も違う。原子数 10 以上の調査例はまだ少ないが、原子数 12 のベンゼン（調査中）で 1500 を超える異性体がすでに見出されているので、同原子数での異性体数のばらつきは原子数の増加とともにさらに拡大する
ことがうかがわれる。
自動探索される TS 数と EQ 数の関係を見ると、

図２ EQ（異性体）数と原子数の相関図３に示したように、おおむね EQ 数の増加につれてTS 数も増加し、TS 数≧EQ 数の関係にあるが、その比は、1 倍程度から 10 倍程度まで、大きな違いがあり、安定構造(EQ)の周りにどの程度の数の反応経路が存在するか、一概に推計できないことがわかる。
図４に、EQ 数に対する計算時間を(hour/core)単位で示した。EQ 数が 20 未満のデータは、計算時間がこのスケールでは無視できるので省略した。１年間=365×

図３ EQ 数と TS 数の相関

24=8760 時間であることを考慮してこの図を見ると、EQ 数が 100 を超えると、1 core での計算時間は、半年～１年以上かかることがわかる。
GRRM 法では、個々の EQ の周囲の反応経路探索を
並行して行うことができるため、多数のコアを利用できれば、計算時間を大幅に短縮できる。ただし、GRRM11 を通常のやり方で用いる場合、１つのノードのコア数の範囲内でしか並列度を高められない。現在開発中の超並列制御プログラム NeoGRRM を用いると、ノードの壁（計算機の箱の壁）を超えて多数のコアを

図４ EQ 数と計算時間の相関有効に利用することができる[3]。NeoGRRM では、
Main-JOB としてまず GRRM11 を起動するが、同時に利用できる他のノード(計算機)に、GRRM11を FirstOnly モードで Sub-JOB として多数投入する。FirstOnly モードでは、１つの EQ 周りの探索しか行わない。NeoGRRM は、Sub-JOB の結果をMain-JOB の結果に統合するととともに、新たな図５ NeoGRRM による GRRM の超並列化

Sub-JOB を、MAIN-JOB の TOTAL EQ-List 中でその周囲が未探索の EQ を Donel List を参照して見つけ出し、空いている計算機に FirstOnly モードで新たな Sub-JOB として投入する。
【結論】GRRM の反応経路網全面探索は、予期せぬ新反応・新構造を与える可能性があり、他の手法にない貴重な特色を有する。見つかる構造(EQ)や TS の数は、原子数との相関が極めて低く、系によって大幅に変化するため予想困難である。計算時間は、構造数(EQ 数)とともに増加する傾向があるが、同じ構造数でも十倍程度の違いがみられる。GRRM の並列処理に
より、計算時間は大幅に短縮される。GRRM の並列処理は、超並列制御プログラム NeoGRRMを用いることでノード（計算機）の壁を超えて可能である。
[1] 大野公一, 前田理, 松本真, 分子構造総合討論会, 京都, 2Tp14 (2003); K. Ohno, S. Maeda, Chem.Phys.Lett. 384,
277 (2004); S. Maeda, K. Ohno, J.Phys.Chem.A 109, 5742 (2005); K. Ohno, S. Maeda, J.Phys.Chem.A 110, 8933 (2006).
[2] 大野公一, 長田有人, 前田理, 諸熊奎治, 理論化学討論会, 岡山, 2D1b (2011).
[3] 大野公一, 理論化学討論会, 仙台, 1P17 (2012); 理論化学討論会, 福岡, 1P05 (2013).
[4] 大野公一, 日本化学会第 83 春季年会, 草津, 2G3-10 (2013); 大野公一, 原渕祐, 前田理, 武次徹也, 化学反応討論会, 仙台, 2P5 (2013).

1E18

グラフェンやカルビンを含む低次元炭素の構造予測（和歌山大システム工 1，和歌山大院システム工 2，量子化学探索研究所 3，
東北大院理 4）○山門 英雄 1， 澤田 裕 2，大野 公一 3, 4

Structure prediction of low dimensional carbon materials including
graphene and carbyne（Faculty of Systems Engineering, Wakayama Univ. 1，Graduate School of Systems
Engineering, Wakayama Univ.2，Institute for Quantum Chemical Exploration 3，Graduate
School of Science, Tohoku Univ. 4）○Hideo Yamakado1，Yu Sawada2，Koichi Ohno3, 4【序】量子化学計算によって、原子や分子を構成要素とする結晶の構造を、予め何の構造情報も与えずに多形も含めて予測できれば、新物質の設計やその物性予測に、非常に有益な手段が提供される。しかし構造を仮定するやり方では、試行回数が膨大となるため、原子の組み換えに対応できる一般的な結晶構造予測法は確立されていない。我々は結晶構造予測に超球面探索法(Scaled
Hypersphere Search method: SHS 法)[1]を適用し、炭素結晶やイオン結晶について、多形を含んだ探索を行ってきた[2]。近年、一般的な 3 次元結晶構造に対して、無限に広がった 2 次元、および 1 次元物質はその低次元性がもつ特異な物性から注目されている。またそれらがどのような構造を取り得るかといった問題は、新規材料の開拓の観点からも興味深い。低次元物質の例として、
炭素の無限平面構造から成るグラフェン、同じく炭素の直線状ポリマーのカルビン等が挙げられる。今回、一般化した超球面探索法[3]を構造予測に用いることにより、炭素の 1 次元および 2 次元構造の探索を行い、また炭素の 3 次元結晶構造予測の結果とも比較した。
【方法】探索の速度を向上させるために半経験的なパラメーターを用いる SCC-DFTB 法でエネルギー値の計算を行った。パラメーターには固体用の pbc-3-0 を使用した。1 次元構造の探索は原子 3,4,5 個を 1 単位とした構造探索をそれぞれ行い、また 2 次元構造の探索は 2,3,4 個を 1 単位とした構造探索をそれぞれ行った。以下では、得られた平衡構造を EQ と呼び、得られたエネルギ
ー値の安定な方から（最安定構造を０番として）通し番号を割り当てた。
【結果と考察】1 次元 C3/chain の探索では 3 種類の構造が得られ、EQ0 は均一な一次元鎖で、他の 2 つは 3 個単位の炭素が並んだ構造であった（図 1）。一方 C4/chain では均一な一次元鎖ではなく、原子間距離が交互に変化する構造が得られている。これは C3 では 2 原子周期の交代が表現できないが、C4 でそれが可能になるためであると考えられる。
続いて、2 次元 C2/sheet ではグラフェンが最安定となり、現実の結果と一致した。また、その他のＥＱとしては、1 次元鎖同士が一定間隔で並ぶ構造が見られた。

図 1 C3/chain で得られた 1 次元構造

次に、平面単位格子内に 3個の炭素原子を置く計算(C3/sheet)を行った（図 2）。
炭素原子 8 個で構成される、変形した六角形がグラフェンのように無限につながる構造が最安定であった。類似の構造が 2 種類独立に見つかり、ジグザグ部分の結合長が均一になっているも

図 2 C3/sheet の探索で探索で得た、８員環を員環を持つ 2 つの構造
つの構造（
構造（数字は
数字は結合距離(Å）。ジグザグ：均一；
不均一
合距離 Å)）。ジグザグ部
ジグザグ部の結合長 EQ0：
均一；EQ1:不均一のがエネルギー的に安定（EQ0）で、交互に長さの異なる結合長を持つものがエネルギー的に不安定(EQ1)である。これはπ電子が構造全体に非局在化した事によって、全エネルギーが安定となっているためと考えられる。また単位格子に 3 原子を入れた場合はグラフェンを表現できないため、この構造が安定構造として出てくるものと考えられる。よりエネルギー的に不安定な構造として C2/sheet と同様に、直線状ポリマーが面
内で平行に並ぶ構造も得られた。
C4/sheet では図 3 のような特徴的な 4 員環-8 員環及び 3
員環-9 員環構造が得られた。既に行った炭素の 3 次元結晶
構造予測において、4 員環-8 員環の層状構造が見つかっており、その構造は 8 員環の真上に 4 員環が、4 員環の真上に 8 員環が来るように平行に重なっていた。この構造は最安定構造のグラファイトよりも相対的に高エネルギーであり、他の多形と比べてより不安定であった。一方で 3 員
環-9 員環構造は 3 次元結晶構造の探索では類似構造を含め見つかっておらず、2 次元平面構造探索での特徴的な構造であるといえる。
【結論】SHS 法を用いて炭素の低次元構造を探索した。単位構造として与える原子の数によって、表現できない構造

図 3 C4/sheet の探索で探索で得た特徴的な
特徴的な 2 つなく、4 員環-8 員環及び 3 員環-9 員環といった特徴的な構 の構造 (a)4 員環-8
員環 員環 (b)3 員環-9
員環 員環
造も取り得る事がわかった。
も存在することが見いだされた。また、6 員環構造だけで[1] K. Ohno, S. Maeda, Chem. Phys. Lett., 348, 277 (2004); S. Maeda, K.Ohno, J. Phys. Chem. A109, 5724
(2005); K. Ohno, S. Maeda, J. Phys. Chem. A110, 8933 (2006).
[2] 山門英雄、時子山宏明、前田理、大野公一、分子科学討論会 2009、2P133; H. Tokoyama, H. Yamakado,
S. Maeda, and K. Ohno, WATOC2011 (17-22 July 2011, Santiago de Compostela, Spain) PIII-065; Yu Sawada,
Hiroaki Tokoyama, Hideo Yamakado, Satoshi Maeda, and Koichi Ohno, 14th ICQC (25-30 June, 2012, Boulder,
Colorado, USA), IV.63 他[3] 大野公一、長田有人、前田理、分子科学討論会 2010、１E15

1E19
スピン渦理論（京大院工）立花 明知

Quantum electron spin vorticity principle(Kyoto Univ.) Akitomo TACHIBANA【序】電子ストレステンソルの起源を Einstein の一般相対性理論に求め [1]、Einstein が発見した「時空の曲率が物体に働く重力を生み出す」、という測地線原理に加えて、「時空のねじれが電子スピンに働くトルクを生み出す」、という量子電子スピン渦原理を発見した[2]：(1)ε Aµν + τ Aµν ( g ) =0
Quantum electron spin vorticity principledsAntisymmetric electronic stress tensor drives the electron spin through vorticity

d
rots


tA
Fig. 1. Quantum electron spin vorticity principle.

本報告では、超重力から導かれるスピン渦とその計算例を相対論的量子力学の範囲で示す。
【理論】超重力ダイナミクスにおける量子電子スピン渦原理は単純超対称性のもとで

ε Aµν ( SUGRA ) + τ Aµν ( SUGRA ) =0となる。弱重力の極限を取れば、超カレントは

=
Θµ

(

(2)

))

(

i
−4 2 ( Φ n †∂ µ Φ n − Φ n ∂ µ Φ n † ) − i ( DΦ n † ) γ µ ( DΦ n ) × c
∑
12 n

(3)

カイラル超場は1
 1+ γ 5  1+ γ 5 F=φn − 2  θ
ψ n  + θθ  Fn − i (θγ 5∂/φnθ )
n
2
2
2

 

2
1
 1− γ 5 1i (θγ 5θ )  θ+∂/ψ n  +  2 (θγ 5θ ) φn
2
2

 8

(4)

拡張ゲージ原理は




(5)
Φ n ( x, θ ) → ∑ exp  i ∑ t AΩ A ( x− )   Φ m ( x, θ )m 
 A  nm(6)Γ ( x, θ ) → exp  −i ∑ t AΩ A ( x− )  Γ ( x, θ ) exp  +i ∑ t AΩ A† ( x− ) 
 A

 A

将来的に超対称粒子の存在が確認されれば、超対称量子電子スピン渦原理が実証される[3]。

【考察】Eq. (2)において Minkowski 時空の極限をとれば、電子スピン s の時間発展は電子ス



トレステンソル τ Π の反対称成分 τ A のみによって導き出されることが示せる：

(

)

d d d
d
d t 1

d
s ( t , r ) − s ( t0 , r ) =2r × ∫ ∫ divt A ( t ', λ r ) λ d λ dt '
t0

0

(7)







この時、電子スピン渦度 rots が重要な役割を演ずる。 rots は電子の運動量 Π を補いもする：d∂cP 0 +c 2 divP =
0
∂t
1

M + h.c.) + H γ 
(



1

µ
µ
2
P=
Pe µ + PEM=
, P + rots + G 

c
2




1
( M + h.c.) 

 1 
H
0
2
, Pe = P + rots , PEM 0 = γ , PEM = G
Pe =
2
c
c
S
Π

∇ν T ν 0 =→0さらに、電子ストレステンソル τ の対称成分 τ

(8)

(9)

から、運動方程式が導かれる：∂ d ∇ν T ν k =→0
P + div (σ − t S ) =
0
(10)
∂t
∂   S(11)
Pe= L + t∂t
同様にして、角運動量 J 保存則が導かれる：
∂ d
d  
(12)
∂ λ M λk = 0 →J + div r × (σ − t S ) = 0
∂t


(13)
J= r × P
dSS
ここに、 σ は Maxwell のストレステンソル、 τ = divτ は張力、 L は Lorentz 力である。

()

【計算例】簡単のために、平面電磁場

 
Ammm
= A (fff) , = k ⋅ x = k 0ct − k • r , lim A ( ) = 0
→ inf
ff

のもとで Dirac 電子がスピン第 3 固有値の漸近値 ζ = ±(14)1
 を持つとする。
2

さらに簡単のために

Aµ = ( 0, Ax , 0, 0 )(15)

k = ( k , 0, 0, k )
µ

0

(16)

0



と仮定すると、第 3 軸方向に運動する電子の電荷密度 j 0 、スピン密度 s 、スピン渦度 rots は
2

1 0
1
2
q
N=
j = 1+ 0 0  ( Ax )cq
2 p ( p − pz )  c 1  1 qs=±   0 Ax , 0,1 − ( N − 1) 
2 p c0
k q dAx 
1 drots =,0±   0, − 0
p c dφ 
2 

(17)
(18)
(19)

で与えられる。
参考文献[1] A. Tachibana, “General relativistic symmetry of electron spin torque,” Journal of Mathematical
Chemistry 50, 669–688 (2012).
[2] A. Tachibana, “Electronic Stress with Spin Vorticity,” In Frontiers in Theoretical Chemistry:
Concepts and Methods, Taylor & Francis / CRC Press, Chap. 12, pp. 235-251 (2013).
[3] A. Tachibana, to be published.

1E20
Rigged QED に基づくシミュレーションにおける効率的な光子相互作用の記述（京大院・工）瀬波 大土, 立花 明知

Eﬃcient computation of interaction with photon in simulation based on Rigged QED(Kyoto Univ.) Masato Senami, Akitomo Tachibana
QED にゲージ不変性を保った形で原子核を追加した Rigged QED[1] に基づいて系の時
間発展を調べる計算コード QEDynamics[2] の開発を行っている。これまで、QED は摂動計算を用いた自由粒子の散乱過程や磁気双極子の予言等で非常に大きな成功をもたらしている。ただし、それらの成功は自由粒子が漸近場として存在するとの仮定が非常に良い場合に限られている。しかし、原子や分子中の電子や原子核は束縛状態としてしか記述できない。束縛状態の記述については、これまでベーテ・サルピーターの方法や NRQED[3] による記述等があり一定の成果を上げているが、まだまだ満足のいく段階に達していない。
そのような状況を踏まえて、我々は場の理論の枠組みで束縛状態の電子・原子核に加えて光子も共に生成消滅演算子で記述して時間発展を取り扱う計算コードの開発に取り組んでいる。本研究では電子と原子核は 2 成分シュレディンガー場として記述することとす
る、この枠組みは Primary Rigged QED[4] と呼ばれている。ローレンツ対称性は用いる
Lagrangian の展開項を無限に取り入れることにより回復できる。
電子・原子核場の時間発展は Heisenberg 描像に基づいて生成・消滅演算子の時間発展として記述される。電子場・原子核場はそれぞれ、Ψ̂, χ̂ として以下のようにあらわす。
∑
Ψ̂(t, ⃗x) =êm (t)ψm (⃗x),(1)
mχ̂a (t, ⃗x) =∑fˆam (t)χam (⃗x),(2)

m

このとき、êm , fˆam が電子・原子核の消滅演算子である。すなわち、電子・原子核の時間発展は生成消滅演算子の時間発展として記述される。展開関数としてば、平面波ではなくFurry 描像で用いられる方法に習って束縛状態の固有関数を用いる。
光子場は次のように分解して表現する、Âµ (t, ⃗x) = Âµrad (t, ⃗x) + ÂµA (t, ⃗x) + ÂµM (t, ⃗x),(3)

ここで Âµrad は輻射光子を、ÂµA,M は系と外部環境の相互作用光子を表す。本研究ではクー
ロンゲージを用いて光子場を表現することにする。輻射光子は光子の生成消滅演算子 â, â†
を用いて、
√ ∑∫
1
c
ν
ν
µ
(4)
Ârad (t, ⃗x) = √d3⃗k √ (â(⃗k, σ)ϵµ e−ikν x /h̄ + â† (⃗k, σ)ϵ∗µ eikν x /h̄ ),2k0
π 2h̄ σ kとあらわされる。一方で ÂµA,M は電子・原子核の生成消滅演算子により展開される。

∫∫
Ψ̂† (t, ⃗s)Ψ̂(t, ⃗s) ∑χ̂† (t, ⃗s)χ̂a (t, ⃗s)Â0A,M (t, ⃗x) =Ze e
d ⃗s
+
Za e
d3⃗s a,|⃗x − ⃗s|
|⃗
x
−
⃗
s
|
A,M
A,M
a
∫
∫
ˆ
⃗
⃗ˆ
1
s) 1
s)
ˆ
3 jeT (u, ⃗3 jN T (u, ⃗⃗AA,M (t, ⃗x) =
+
,
d ⃗s
d ⃗s
c A,M
|⃗x − ⃗s|c A,M|⃗x − ⃗s|3

ここで u = t − |⃗r − ⃗s|/c は遅延効果を表している。

(5)

(6)

電子・原子核場の時間発展は、以下の運動方程式に従って発展する。
[
]
(
)2
∑∫
d
Ze e i13†−ih̄∂i −ih̄ êm (t) =d ⃗xψm (⃗x)
Â+ Ze eÂ0 ψn (⃗x)ên (t),
dt
2me
c
n
]
[
(
)2
∫
∑
d
Za e i1
ih̄ fˆam (t) =d3⃗xχ†am (⃗x)−ih̄∂i −Â+ Za eÂ0 χan (⃗x)fˆan (t).
dt
2m
c
a
n

(7)
(8)

この式は演算子を含む式であり、数値計算に取り入れるために t = 0 における演算子の多項式を用いて時刻 t の演算子を展開し、その多項式の係数を各時刻で計算するという手法を用いる。時間発展のステップ数に応じて、多項式中の高次の項の組合せの寄与が入ってくるため、展開に用いなくてはならない多項式の種類が急速に増える。このため、1 ス
テップを計算する際に要する演算時間も増大するし、それらの情報を記録するデータ自体も巨大に膨れ上がる。そのため本研究では同種粒子の演算子を 3 次まで、これより高次の項や異種粒子の演算子は期待値を取って定数によって置き換えて計算を行う。
本研究では相互作用する光子特に ÂiA をどのように計算するべきかを報告する。ÂiA の計算には過去の演算子の情報も必要なため正確な計算のためには膨大な計算量が必要となる。そのため近似を用いて効率的に計算を行う。具体的には ÂiA の計算は期待値として行う、またその期待値も空間を離散化した点上で行う。時間的にも各計算ステップごとにÂiA の計算を行わずに、数ステップに一度計算を行い、その間は同じ値を使用して計算量を減らすこととする。この方法の利点は計算量を上げて、計算間隔を無限に小さくすることによって正確な計算が行える点である。
QED としての計算を行うには、QED に対して定義されたハミルトニアンを作るためのthermalization 過程が必要である。これまでの研究の結果として、現行の thermalization
過程には膨大な計算時間が必要であることがわかった。この原因の一つとして、Â0 は無限の光子交換により平均化され thermalize されたポテンシャルであるが、Âi は ⟨Âi ⟩ = 0
から摂動的な相互作用を繰り返して thermalize を行うものである。この摂動的相互作用を無限回繰り返すことが thermalization であり、現行の計算方法では非常に長時間の計算が必要であり、現実的には thermalization を完了できない。これは、非摂動系である束縛状態に対して、ÂiA の計算に摂動的手法を用いているのが原因と言える。このため ĵ i を Âiから、ĵLi を Â0 から計算すると、現状の計算では ĵ i の中に、ĵLi が含まれていないという矛盾した結果となってしまっている。この問題を解決するために摂動的な thermalization
ではなく新たな非摂動手法を用いる。その方法では ĵ i = ĵLi + ĵTi の関係式が正しく再現されることを目的とする。本講演ではこの新たな thermalization 方法について報告する。

参考文献[1] A. Tachibana, in Fundamental World of Quantum Chemistry, A Tribute to theMemory of Per-Olov Lödin, ed. by E. J. Brändas and E. S. Kryachko, (Kluwer
Academic, Dordrecht, 2003) Vol. 2, 211; J. Mol. Modelling 11, 301 (2005); J. Mol.
Struct. (THEOCHEM), 943, 138 (2010); J. Math. Chem. 50, 669 (2012).
[2] K. Ichikawa, M. Fukuda, A. Tachibana, Int. J. Quant. Chem. 113 (2013) 190; M.
Senami, T. Miyazato, S. Takada, Y. Ikeda, and A. Tachibana, to be pulished in J.
Phys. Conf. Ser.; M. Senami, Y. Ogiso, T. Miyazato, F. Yoshino, Y. Ikeda, and A.
Tachibana, to be pulished in Trans. Mat. Res. Soc. Jpn.; QEDynamics, M. Senami,
K. Ichikawa, and A. Tachibana, (http://www.tachibana.kues.kyoto-u.ac.jp/qed/)[3] W. E. Caswell and G. P. Lepage, Phys. Lett. 167B, 437 (1986).
[4] A. Tachibana, “Electronic Stress with Spin Vorticity”, In Concepts and Methods in Modern Theoretical Chemistry: Electronic Structure and Reactivity (Atoms,
Molecules, and Clusters), CRC Press (2013) 235.

1E21
Rigged QED による時間発展シミュレーションにおけるベクトルポテンシャルの効果について（京大院工）○市川 和秀，福田 将大，立花 明知

Eﬀect of vector potential in time evolution simulation of Rigged QED(Kyoto University) ○ Kazuhide Ichikawa, Masahiro Fukuda, Akitomo Tachibana

われわれは原子分子系と光の相互作用の時間発展を量子電磁力学（QED：Quantum ElectroDynamics）に基づいて場の量子論的に計算するための定式化と数値計算法を研究している。Rigged
QED とは電子場・光子場に加えて原子核場を量子場として取り入れたものである [1]。これまでに、この理論のもとで電子や原子核の生成消滅演算子の時間発展の式を密度行列の時間発展で近似するという定式化および数値計算を行っており、電荷密度に、仮想電子陽電子対生成と対消滅に起因していると考えられる「電子陽電子振動」という高速の振動が見られることを見いだした [2, 3]。
通常の QED と異なる点として、電子が束縛状態にある設定で時々刻々物理量の時間発展を計算したいということがある。われわれはこれらの問題を、電子場を束縛基底で展開して生成消滅演
算子を定義し、ハイゼンベルク表示で時間発展を計算することで対処している。具体的には、電∑ D ∑子陽電子を表すディラック場を ψ̂(ct, ~r) = Nr) のように、外場中での Diracn=1
a=± êna (t)ψna (~方程式の規格化された解 ψna (~r)（a = +, − がそれぞれ n 番目の電子、陽電子軌道を表す）で展開
して電子の消滅演算子 ên+ (t) と陽電子の生成演算子 ên− (t) を定義する。これらを用いて励起演算†子 Êpc qd ≡ êpc êqd を定義すると、電荷密度演算子など物理量の演算子は励起演算子で表されるが、この励起演算子について、時間発展をハイゼンベルク表示で計算する。実際には、励起演算子を初期状態 |Φi で期待値をとった密度行列 Epc qd (t) ≡ hΦ|Êpc qd (t)|Φi の時間発展の式を解いて、物理量演算子の期待値の時間発展を計算している。
光子場については、われわれはクーロンゲージでの正準量子化を採用している。ハイゼンベル
ク表示では、光子場のベクトルポテンシャル部分は自由な輻射場が量子化された部分と相互作用している物質場からの寄与の部分の遅延ポテンシャルの和で表される。静電場極限のハミルトニ
アンはこのベクトルポテンシャルを落としたもので、実在系を表現するためにはこの寄与を取り入れることが必要である。具体的には以下のようになる。
~ˆrad (ct, ~r) + A~ˆA (ct, ~r)
= A
√
∫
]
0
0
4π~2 c ∑d3 p~ [√Âkrad (~r) = √â(~p, σ)ek (~p, σ)e−icp t/~ ei~p·~r/~ + â† (~p, σ)e∗k (~p, σ)eicp t/~ e−i~p·~r/~(2π~)3 σ=±1
2p0∫
~ˆ~ˆA (ct, ~r) = 1 d3~s jT (cu, ~s) , u = t − |~r − ~s|Ac|~r − ~s|
c
~ˆ
A(ct,
~r)

ここで â(~p, σ) は光子の消滅演算子で、p~ は光子運動量、σ は左右の円偏光を表し、ek は偏光ベクˆˆトルである。~jT は電流演算子の横波成分（div~jT (~r) = 0）である。遅延時間 u については、物質場が t < 0 でゼロであるとして、次のように変形して扱う [4]。
~ˆA (ct, ~r) =
A

1
c2 π

∫ t

∫ ∞
du

0

−∞

(
)
dα exp iα(u − t)2∫()(~r − ~s)2ˆ~d ~s jT (cu, ~s) exp −iαc23~ˆA や ~ˆjT の演算子を、|Φi で期待値をとったものとして近似的に扱う。
また、以下では Aこれにより、電子の密度行列の時間発展方程式は、ボルンオッペンハイマー近似のもとで、自由輻射場の寄与を省略して書くと次のようになる。
dEna mbdt= Ona mb + On† a mb

i~Ona mb=

ND ∑∑r=1 e=±+

ND∑

hmb re Ena re +∑()(mb re |pc q d ) Ena re Epc qd − Ena qd Epc rer,p,q=1 e,c,d=±

ND ∑∑IjT mb re (t)Ena re

r=1 e=±

ここで、
hna mb= Tna mb + Mna mb +~ a ),(Za e)Vna mb (Ra=1∫
Tna mbNn∑d ~r ψn† a (~r)γ 0 γ i ∂i ψmb (~r),= −i~c3∫~
Vna mb (R)= (Ze e)d3~s∫(na mb |pc q d ) = (Ze e)2∫
Mna mb = me c2

d3~r ψn† a (~r)γ 0 ψmb (~r)ψn† a (~s)ψmb (~s)~|~s − R|

d3~r d3~s ψn† a (~r)ψmb (~r)

1ψ †c (~s)ψqd (~s)|~r − ~s| p∫ ∞ND ∑ ∫ t()
1 ∑IjT mb re (t) = − 3dudα exp iα(u − t)2c π p,q=1−∞c,d=± 0{}dEpc qd
× Ijj,mb re pc qd Epc qd (u) + IjE,mb re pc qd
(u)
dt
(
)
∫
3
∑
(~r − ~s)2kIjj,mb re pc qd =d3~r d3~s jmr)jpkc qd (~s) exp −iαb r e (~
c2
k=1
(
)
3 ∫
∑
(~r − ~s)2
3
3kk
IjE,mb re pc qd =d ~r d ~s jmb re (~r)Epc qd (~s) exp −iα
c2
k=1
∫
]
[
Ze e
(~t − ~s)kjpka qb (~r) = Ze e c ψp†a (~r)γ 0 γ k ψqb (~r) , Epkc qd (~s) = −d3~t ψp†c (~t)ψqd (~t)
4π|~t − ~s|3

本発表では光子場のうちでベクトルポテンシャル、特に遅延ポテンシャルの寄与を計算するための方法を議論する。近似方法や数値計算方法を説明し、初期値でゼロと置いたベクトルポテン
シャルがどのように時間発展し、電子場から構成される電流密度と自己無撞着になるかを議論する。また、電子電荷密度などの物理量の時間発展にどのように影響するかを議論する。

参考文献[1] A. Tachibana, J. Chem. Phys. 115, 3497 (2001); in Fundamental World of Quantum Chemistry,
A Tribute to the Memory of Per-Olov Löwdin, ed. by E. J. Brändas and E. S. Kryachko, (Kluwer
Academic, Dordrecht, 2003), Vol. 2, p. 211; J. Mol. Struct. (THEOCHEM), 943, 138 (2010).
[2] K. Ichikawa, M. Fukuda and A. Tachibana, Int. J. Quant. Chem. 113, 190 (2013)[3] QEDynamics, M. Senami, K. Ichikawa and A. Tachibana
http://www.tachibana.kues.kyoto-u.ac.jp/qed/index.html[4] “Electronic Stress with Spin Vorticity,” A. Tachibana, In Concepts and Methods in Modern Theoretical Chemistry: Electronic Structure and Reactivity (Atoms, Molecules, and Clusters): CRC
Press, Chap. 12, pp. 235-251 (2013).