量子電磁力学に基づく原子分子系における時間発展の計算方法の研究

市川 和秀, 福田 将大, 立花 明知
京大院工マイクロエンジニアリング

熱場の量子論とその応用, 京大基研, 2013年8月28日

近年のミクロスケールの実験の進展
・ナノテクノロジー・アトムテクノロジー
原子・分子レベルで構造を測定・作成・制御
・Cavity QED・スピントロニクス
単一光子・単一電子の生成・測定
スピンの操作・レーザー科学
フェムト秒・アト秒スケールでの物理現象の観測
時間方向・空間方向・粒子の属性というあらゆる意味でより基本的なレベルでの実験・観測が可能となってきている。
より基本的なミクロ法則に基づく理論計算の必要性

従来のシミュレーション方法
・ 時間依存シュレディンガー or ディラック方程式 or 密度汎関数(DFT)法に古典電磁場を入れる。（半古典近似）

e.g.) laser
・ 量子化された電磁場と、モデル化された物質場を用いる。（相互作用などが簡単化されている。）

e.g.) cavity QED
電子は量子化されたディラック場（波動関数ではなく）として扱い、他の電子や原子核との相互作用は光子場（古典電磁場ではなく）によって、つまり量子電磁力学（QuantumElectroDynamics, QED）によって取り扱いたい。
物質と光を記述するうえで最も基本的な理論であるQEDに基づいて系の時間発展を計算機上でシミュレートする方法を確立し、物性や化学の研究に応用したい。

電子ストレステンソルによる化学結合の理解

Tachibana, J. Chem. Phys. 115, 3497 (2001)

i :自然軌道
非相対論極限で、静電ハミルトニアンの定常状態で期待値をとる

このストレステンソルの固有値・固有ベクトルなどをKI, Nozaki, Komazawa,Tachibama, AIP
advances 2, 042195 (2012)

用いた化学結合の研究
10.001

水素分子

00H

H

-0.50x (A)0.51

正の最大固有値→共有結合性

0.001

Li(3)

Liクラスター
2-0.001-0.5-1-1301-0.001

LiLi-0.002

y (A)

y (A)

0.5

⌫i :占有数-0.003-1-0.003-0.004-2-0.004-0.005-3-30(1)-0.002(2)

Li-2-10x (A)123-0.005

負の固有値が縮退→金属結合性

QEDハミルトニアンを使って計算したら？ 時間変化（化学反応時のふるまい）は？

電子ストレステンソルによる化学結合の理解

Tachibana, J. Chem. Phys. 115, 3497 (2001)

i :自然軌道
非相対論極限で、静電ハミルトニアンの定常状態で期待値をとる

このストレステンソルの固有値・固有ベクトルなどをKI, Nozaki, Komazawa,Tachibama, AIP
advances 2, 042195 (2012)

用いた化学結合の研究
10.001

水素分子00H

H-0.50x (A)0.51

正の最大固有値→共有結合性0.001

Li(3)Liクラスター
2-0.001-0.5-1-1

301-0.001

LiLi-0.002

y (A)

y (A)

0.5

⌫i :占有数-0.003-1-0.003-0.004-2-0.004-0.005-3-30(1)-0.002(2)

Li-2-10x (A)123-0.005

負の固有値が縮退→金属結合性

QEDハミルトニアンを使って計算したら？ 時間変化（化学反応時のふるまい）は？

QEDに基づく原子分子系の時間発展シミュレーションを行うためには、計算方法が確立している従来のQEDでは現れなかった問題が生じる。
これら三つの要素をそれぞれ部分的に考慮する場の理論的手法は存在するが、三つ同時に扱う方法はQEDでは知られていない。
(i) 原子・分子系における物質粒子は束縛状態にある。
ベーテ・サルピータ方程式（相対論的な粒子の場合は適用外）(ii) 原子核が存在する。
有効場の理論（散乱理論）
格子場の理論（平衡状態の記述。また、QEDは難しい。）(iii) 系の時間変化を時々刻々追跡する。
閉時間（closed time path, CTP）形式 （相互作用表示による摂動展開にThermo Field Dynamics基づく。束縛状態の取り扱い？）

われわれのアプローチ(i) 原子・分子系における物質粒子は束縛状態にある。
場の演算子を、平面波ではなく外場存在下での量子力学の定常解など波束（束縛基底）で展開して生成消滅演算子を定義する。
（いわゆるファリー表示と類似） Furry, Phys. Rev. 81, 115 (1951)(ii) 原子核が存在する。

Tachibana, J. Chem. Phys. 115, 3497 (2001)

原子核の自由度をシュレディンガー場として取り入れる。
“Rigged”（「追加された」という意味合い）QEDという名前はここから来ている。
核スピンは、核磁気モーメントを生成する「電流」を通常の電流と同様に取り扱うという現象論的な扱いをする。

(iii) 系の時間変化を時々刻々追跡する。

KI, Fukuda and Tachibana, Int. J. Quan. Chem.
113, 190 (2013)

場の演算子をハイゼンベルク表示の演算子として扱い、場の運動方程式に従って時間発展を追う。
（束縛状態を扱うため、摂動のゼロ次にとるべきハミルトニアンが不明。
有効な相互作用表示が不明。）

場の演算子の運動方程式：
光子場 (Maxwell equations) :(Coulomb gauge)
電子場 (Dirac equation) :{ ˆ (x), ˆ† (x )} = (rr)

場の演算子の電荷密度:

交換関係

電流密度:

原子核場 (Schrödinger equation) :[ ˆa (x), ˆ†a (x0 )]± = (~r~r0 )場の演算子の運動方程式：
光子場 (Maxwell equations) :(Coulomb gauge)
電子場 (Dirac equation) :{ ˆ (x), ˆ† (x )} = (rr)

(従来の) QED

電荷密度:
電流密度:

原子核場 (Schrödinger equation) :[ ˆa (x), ˆ†a (x0 )]± = (~r~r0 )原子核部分

場の演算子を有限個の束縛基底で展開する。

KI, Fukuda and Tachibana, Int. J. Quan. Chem.
113, 190 (2013) arXiv:1204.1724

電子の場の演算子→原子核による外場存在下でのディラック方程式の4成分Hartree-Fock方程式を解で展開。
解くことで得られる軌道cf) Furry, Phys. Rev. 81, 115 (1951)

電子消滅演算子

陽電子生成演算子

式を簡単にするために以下の記法を導入。

時間依存性はエネルギー固有値で決めるわけではない（以下を参照）。

{ ˆ (x), ˆ† (x )} = (rr)

光子場は電子場・原子核場と、光子の生成消滅演算子で書ける。

・・

光子の消滅演算子
自由場

光子の生成演算子
遅延ポテンシャル

電荷密度:
電流密度:

偏光ベクトルどんな量を計算するか。
電子電荷密度Ze =演算子:

ここで、1電子励起演算子(たんに座標の関数)
励起演算子の時間発展により、物理量演算子の時間発展が記述される。
物理量演算子の期待値の時間発展が記述するには、励起演算子の期待値(密度行列と呼ぶ)の時間発展がわかればよい。

生成消滅演算子の時間発展の式
電子場
場の演算子の展開を代入して、消滅演算子の時間発展の式を得る。
（生成演算子の式はエルミート共役から得られる。）
原子核場も同様。

以下ではBorn-Oppenheimber近似（原子核を固定）の話に限る。

量子化学でよくでてくる分子積分を係数に持つ微分方程式
運動エネルギー積分
質量エネルギー積分

核引力積分
二電子積分

励起演算子の時間微分を模式的に書くと、d †dê†dê†(ê ê) =ê + êdtdtdt†† †〜 ê ê + ê ê ê ê† ††+ ê ĵT [ê, â]ê + ê â ê + ê â ê
†

両辺の期待値をとってもhê† êi についての閉じた時間発展方程式は得られない。

（さらに高次の演算子の組み合わせが増えてしまうだけ。）

・４次の項を近似的にfactorizeする。
††† †hêêihêShêP êQ êR êS i 〜 P
Q êR i††hêP êR ihêQ êS i 等。

・電流密度演算子 ĵT [ê, ĉ, â] は期待値で置き換える。

†・ ê â ê のような項は、†
方法① hê êihâi とfactorizeする近似。
方法②ê† â ê の微分方程式と連立させて計算。

励起演算子の期待値の時間微分方程式は、のような記号を定義すると、質量エネルギー積分運動エネルギー積分

核引力積分

二電子積分

今回はこの項を無視する。

ê† âêê† â† ê

励起演算子の期待値の時間微分方程式は、のような記号を定義すると、質量エネルギー積分

運動エネルギー積分

核引力積分

二電子積分

今回はこの項を無視する。

ê† âêê† â† ê† hâiとfactorizeする近似。
hêêi
方法①

励起演算子の期待値の時間微分方程式は、のような記号を定義すると、

質量エネルギー積分

運動エネルギー積分

核引力積分

二電子積分

†
方法① hê êihâi とfactorizeする近似。

ここではの single modeとみなせるような光をあらわす、を用いる。

光子の状態がコヒーレント状態の時のみ有限の寄与をする。

水素原子の電子電荷密度の時間発展の計算

ˆe (t, r) =電子の軌道はDIRAC10コードにより Dirac Hartree Fock法で、基底関数STO-3Gを用いて計算した。

初期条件：基底状態| i = ê†1+ |0i
電子の占有軌道
陽電子の全軌道光子の状態は、光子が無い場合と、コヒーレント状態（古典的な外部振動電場が存在することに対応）で計算。

電子

1+1̄+陽電子

11̄水素原子において、核から21 a.u.離れたところでの電荷~1.7x10-4 ~ 2⇡/(2me c )密度の時間変化（初期値からのずれをプロット）。
0という電子質量の二倍に対応する周期の振動が電荷密-14-1x10度に見られる。
（仮想）電子陽電子の対消-14-2x10

滅・対生成に起因する電荷密度の揺らぎと考えられる。

-14-3x1000.001

0.002

0.003

0.004

0.005

1 a.u. = 24.19 as = 2.419 x 10-17 s

光子の初期状態をモードp0（振動電場の振動数に対応）

外部電場による変調のコヒーレント状態にとる。
0p=10無し
0-9-1x10-14-1x10-9-2x10-14-2x10-9-3x10-14-3x1000.001

0.002

0.003

0.004

0.005

0-12

p = 1002x10-81x100.001

0.002

0.003

0.004

0.005

0.002

0.003

0.004

0.005

p0 = 20-121x1000-12-1x10-8-1x10

2 /(cp0 ) = 4.59
000.001

10 30.002

0.003

0.004

0.005

0.005

-12-2x1000.001

励起演算子の期待値の時間微分方程式は、のような記号を定義すると、

質量エネルギー積分

運動エネルギー積分

核引力積分

二電子積分

今回はこの項を無視する。

ê† âêê† â† ê
方法② ê† â ê の微分方程式と連立させて計算。

hê† âêi の時間発展の式

消滅演算子の時間発展はDirac方程式よりわかる（前出）※ほとんど同じ。

0t=0でのみ成り立つ関係式で因子分解をする近似。

np~ : モード(~p, ) の光子の占有数励起演算子の期待値の時間微分方程式は、のような記号を定義すると、真空（全てのモードの占有数がゼロ）でも有限の寄与を与える。

水素原子の電子電荷密度の時間発展の計算（自己エネルギー入り）

ˆe (t, r) =電子の軌道はDIRAC10コードにより Dirac Hartree Fock法で、基底関数STO-3Gを用いて計算した。

初期条件：基底状態| i = ê†1+ |0i
電子の占有軌道
陽電子の全軌道電子

1+1̄+陽電子

11̄

初期状態にphotonは無し。
00<p 10 に１０コ光子のモードは、0  ✓  ⇡ に５コ0  < 2⇡ に４コを各偏光ごとにとっている。

(x,y,z)=(-1,0,0) での電子電荷密度をプロット。 初期値からの変化がプロットされている。
1.2e-132~1.7x10-4 ~ 2⇡/(2me c )no self-energy
with self-energy（10-4 倍している。）

1e-13

charge density (a.u.)

8e-14

自己エネルギー無し

6e-14
4e-14
2e-14
0-2e-14

自己エネルギー有り-4e-14-6e-14-8e-1400.001

0.002

0.003

自己エネルギーの存在によって、振動数が増えている。
→電子質量が大きくなったように見える。

0.004

0.005
time (a.u.)<別の系での例> 水素分子の２つの原子核の中点で、電子電荷密度の初期値からの変化をプロット。

9e-10xH

no self-energy
with self-energy（10-2 倍している。）

8e-10

2~1.7x10-4 ~ 2⇡/(2me c )7e-10

charge density (a.u.)H

自己エネルギー有り

6e-10
5e-10
4e-10
3e-10
2e-10

自己エネルギー無し

1e-10
0-1e-1000.001

0.002

0.003

自己エネルギーの存在によって、振動数が増えている。
→電子質量が大きくなったように見える。

0.004

0.005
time (a.u.)まとめ・ QEDに基づいて原子分子系、すなわち原子核が存在して束縛状態にあるような系の時間発展を時々刻々シミュレートするために、(i) 場の演算子を束縛基底で展開する。
(ii) 原子核場をシュレディンガー場として追加する。
(iii) ハイゼンベルク表示で場の時間発展を追う。という方法を用いて、いくつかの近似の下で数値計算を行った。

・ 生成消滅演算子・励起演算子の時間発展の方程式を導出し、励起演算子の期待値（密度行列）の微分方程式を導出した。光子場において遅延ポテンシャルを無視する近似およびBorn-Oppenheimer近似のもとで、2自己エネルギーの寄与が無い場合は、電子電荷密度に 2⇡/(2me c )~ 1.67x10-4 a.u. の周期の振動が見られた。

・ 自己エネルギーの寄与を入れると、電子陽電子振動の振動数が増加した。この振動は物理的には電子質量で決まるものであると考えられる
ため、電子質量が自己エネルギーの効果によって増加したようにみえる。