2022
19

5

16
3

Journal of Atomic Collision Research

Vol.19 Issue32022

原⼦衝突学会｜www.atomiccollision.jp

原子衝突学会賛助会員（五十音順）

アイオーピー・パブリッシング・リミテッド（IOP 英国物理学会出版局）
http://iopscience.iop.org/
アドキャップバキュームテクノロジー株式会社
http://www.adcap-vacuum.com
有限会社イーオーアール

Electronics Optics Research Ltd.

http://www.eor.jp/

イノベーションサイエンス株式会社
http://www.innovation-science.co.jp/
株式会社オプティマ
http://www.optimacorp.co.jp/
クリムゾン インタラクティブ・ジャパン
http://www.enago.jp/
http://ulatus.jp/
http://www.voxtab.jp/
コスモ・テック株式会社
https://www.cosmotec-co.jp/
株式会社サイエンス ラボラトリーズ
http://www.scilab.co.jp/
真空光学株式会社

http://www.shinku-kogaku.co.jp/

スペクトラ・フィジックス株式会社
http://www.spectra-physics.jp/
ソーラボジャパン株式会社
http://www.thorlabs.jp/
株式会社東京インスツルメンツ

http://www.tokyoinst.co.jp/
株式会社ナバテック
http://www.navatec.co.jp/

仁木工芸株式会社

http://www.nikiglass.co.jp/
伯東株式会社
http://www.g5-hakuto.jp/

フォトテクニカ株式会社
https://www.phototechnica.co.jp/

第 19 巻 第 3 号

目

次

巻頭言

田沼肇

… 34

解説
高強度レーザー場中の多電子ダイナミクスのための時間依存多配置理論
および時間依存結合クラスター理論

佐藤健

… 36

行事委員会 … 49

第 29 回原子衝突セミナー報告
第 29 回原子衝突セミナー参加報告

手塚智哉 … 50

第 29 回原子衝突セミナー参加報告野村絢也 … 50
庶務委員会 … 51

2022 年度第 1 回運営委員会（文書持回）開催報告

広報渉外委員会 … 52

国際会議 IVC-22（日本表面真空学会主催）のお知らせ【訂正】原子衝突学会第 46 回年会 ホットトピック講演の学生会員発表者

顕彰委員会 … 52
「しょうとつ」原稿募集
ユーザー名とパスワード

編集委員会事務局 … 53
… 53

解

説

高強度レーザー場中の多電子ダイナミクスのための時間依存多配置理論
および時間依存結合クラスター理論
佐藤健 1, 2, 3
1
2

東京大学大学院工学系研究科原子力国際専攻 〒113-8656 東京都文京区本郷 7-3-1

東京大学大学院工学系研究科附属光量子科学研究センター 〒113-8656 東京都文京区本郷 7-3-1
3

東京大学光量子科学連携研究機構 〒113-0033 東京都文京区本郷 7-3-1
sato@atto.t.u-tokyo.ac.jp
令和 4 年 2 月 5 日原稿受付

高強度レーザー場中の原子・分子におけるトンネルイオン化や高次高調波発生などの高強度場現象を第一原理的に記述するための，時間依存多配置理論および時間依存結合クラスター理論を概観する．
これらの理論では，全電子波動関数を配置間相互作用展開や結合クラスター展開で表現し，展開の基底となる一電子軌道関数も時間発展させることにより，系統的かつ柔軟な近似が可能で，束縛状態から連続状態への遷移を記述でき，ゲージ不変性とエーレンフェストの定理が満たされる．これらの長所を活かし，高強度場現象を精度良く計算するために欠かせないプログラム実装についても，数値計算例を交えながら詳しく解説する．

1．はじめにべきである．また実験と直結するイオン化確率

高強度レーザー場中の原子・分子の超閾イオ
ン化，トンネルイオン化，高次高調波発生（HHG），

や光電子スペクトルを抽出できることが望ま

非逐次多重電離（NSDI）などの高強度場現象をレーザーのパラメータ）や計算資源（許される研究する分野が高強度場物理（強光子場科学）
である[1–3]．高強度場現象を第一原理的に記述

計算時間やメモリサイズ）に応じて最適な計算するには，レーザー場中の多電子系に対する時間依存シュレーディンガー方程式（TDSE）を解く必要がある．
𝜕) (𝑡)Ψ(𝑡)(1)𝑖 Ψ(𝑡) = 𝐻𝜕𝑡) (𝑡)はハミルトニ
ここでΨ(𝑡)は系の波動関数，𝐻しい．(i)は計算対象（原子・分子種や照射する手法を柔軟に選んだり，計算精度を系統的に改善することで計算結果の信頼性を吟味したりする際に重要である．また異なる近似レベルでの計算結果を比較することで実験では得難い
物理的洞察が得られることもある．(ii)は高強度レーザー照射下のイオン化過程を記述する

アンである．しかし電磁場下の多電子系でため，(iii)は外部電磁場中の観測量を適切に定TDSE を厳密に解くのは困難で，近似的な時間

義・計算するために重要である．

依存多電子理論が必要である．
近似的な時間依存多電子理論に求められる要件として以下が挙げられる．すなわち(i)

系統的かつ柔軟な近似が可能で，(ii) 束縛状態から連続状態への遷移を記述
でき，(iii) ゲージ不変性とエーレンフェストの定

これらの条件を満足する理論的枠組のひとつが，時間依存多配置自己無撞着場法（TDMCSCF）である[4–9]．TD-MCSCF では全電子波動関数をスレーター行列式（電子配置）Φ! (𝑡)の重ね合わせで表現する．
Ψ(𝑡) = + 𝐶! (𝑡)Φ! (𝑡)

(2)

!

理を満たす

原子衝突学会誌しょうとつ 第19巻第3号 (2022)36困難がある．式(4)の展開を打ち切ることにより（たとえば𝐶1% (𝑡)まで）計算コストを削減でき

sha1_base64="za/pIuqwAhBCibXdAmT1Cff26bE=">AAACbHichVG7SgNBFD1ZXzG+4qMQRBBDJDZhVoKKVdDGMlHjAxPC7jrq4L7YnQRi8AdsLSzUQkFE/Awbf8DCTxDBRsHGwpvNgqiod5iZM2fuuXNmRndN4UvGHiJKS2tbe0e0M9bV3dPbF+8fWPWdimfwguGYjreuaz43hc0LUkiTr7se1yzd5Gv63kJjf63KPV849oqsubxkaTu22BaGJonaKLq+KGdScrIcT7A0C2LsJ1BDkEAYOSd+hSK24MBABRY4bEjCJjT41DahgsElroQ6cR4hEexzHCBG2gplccrQiN2jcYdWmyFr07pR0w/UBp1iUvdIOYYku2fX7IXdsRv2yN5/rVUPajS81GjWm1rulvsOh5ff/lVZNEvsfqr+9CyxjdnAqyDvbsA0bmE09dX945fluaVkfYJdsCfyf84e2C3dwK6+Gpd5vnSCGH2A+v25f4LVqbQ6nVbzmUR2PvyKKEYwjhS99wyyWEQOBTrXwhFOcRZ5VoaUEWW0mapEQs0gvoQy8QFBcIzW</latexit><latexit⋮
4 (t)1,2,···X,k(t) = C1(t)るが，打ち切った CI 展開ではサイズ無矛盾性

sha1_base64="t6sDHQLYFpr9i9940UrWh/6JRSs=">AAACbHichVHLSgMxFD0dX7W+6mMhiFAsim5KpoiKq6Ibl621VlEpM2PU4LyYSQu1+ANuXbhQFwoi4me48Qdc9BNEcKPgxoW30wFRUW9IcnJyz81Jorum8CVj9YjS0trW3hHtjHV19/T2xfsHVn2n7Bm8YDim463pms9NYfOCFNLka67HNUs3eVHfX2zsFyvc84Vjr8iqy7csbdcWO8LQJFHrm64vSulJOVWKJ1mKBZH4CdQQJBFG1olfYxPbcGCgDAscNiRhExp8ahtQweASt4UacR4hEexzHCJG2jJlccrQiN2ncZdWGyFr07pR0w/UBp1iUvdImcA4e2A37IXds1v2yN5/rVULajS8VGnWm1rulvqOhvNv/6osmiX2PlV/epbYwVzgVZB3N2AatzCa+srByUt+fnm8NsEu2RP5v2B1dkc3sCuvxlWOL58iRh+gfn/un2A1nVJnUmpuOplZCL8iihGMYZLeexYZLCGLAp1r4RhnOI88K0PKiDLaTFUioWYQX0KZ+AA9aozU</latexit>

が満たされない[10]．つまり電子数が増えると<latexit

i1 ,i2 ,··· ,iN

Ci1 i2 ···i(t) det [ i1+C(t)3(t)i2 (t) · · ·+ ⋯+CiN (t)]
N 2(t)

sha1_base64="Y5cQ5oELkVM+/O8oUjcSITnjopI=">AAACbHichVG7SgNBFD1ZXzE+Eh+FIIIYFG3CrIiKVdDGMj5iIomE3XXUIftidxKIwR+wtbBQCwUR8TNs/AELP0EEGwUbC282C6Ki3mFmzpy5586ZGd01hS8Ze4goLa1t7R3RzlhXd09vPNHXv+E7Fc/gWcMxHS+vaz43hc2zUkiT512Pa5Zu8pxeXmrs56rc84Vjr8uay7csbdcWO8LQJFGbRdcXJXVSTpUSSZZiQYz+BGoIkggj4ySuUMQ2HBiowAKHDUnYhAafWgEqGFzitlAnziMkgn2OA8RIW6EsThkasWUad2lVCFmb1o2afqA26BSTukfKUYyze3bNXtgdu2GP7P3XWvWgRsNLjWa9qeVuKX44tPb2r8qiWWLvU/WnZ4kdzAdeBXl3A6ZxC6Opr+4fv6wtrI7XJ9gFeyL/5+yB3dIN7OqrcbnCV08Qow9Qvz/3T7AxnVJnU+rKTDK9GH5FFMMYwyS99xzSWEYGWTrXwhFOcRZ5VgaVYWWkmapEQs0AvoQy8QE7Z4zT</latexit><latexit

2 (t)

!

0,1
X

n1 ,n2 ,··· ,nk
1 (t)Cn1 n2 ···nk (t)|n1 i ⌦ |n2 i · · · ⌦ |nk i

初期軌道

計算精度が劣化してしまう．そこで次式のクラスター展開を考える．
Ψ(𝑡) = 𝑒 &'()) Φ(𝑡)(5)原子・分子物理や量子化学分野の結合クラスタ

レーザー照射ー（CC）法で用いられるのと同型の波動関数で

図 1: TD-MCSCF 波動関数の模式図．

ある．CC 法は展開を打ち切ってもサイズ無矛

原子・分子物理や量子化学の配置間相互作用（CI）法で用いられるのと同型の波動関数である．TD-MCSCF の特徴は，CI 係数𝐶! (𝑡)とスレーター行列式Φ! (𝑡)を構成する軌道関数-𝜓" (𝑡)/をともに時間発展させる点である（図 1 参照）．
式(2)（図 1）で，与えられた数の軌道関数に電子を配置するすべての組み合わせを取り入
れ る 手 法 が 多 配 置 時 間 依 存 Hartree-Fock（MCTDHF）法であり，軌道関数の数を増やす
ことにより計算精度が系統的に改善され
TDSE の解に収束する（要件 i）．CI 係数だけでなく軌道関数も時間発展させることにより，固定した軌道関数を用いる場合よりも少ない軌
道関数でイオン化過程を記述できる（要件 ii）．
また時間依存変分原理（TDVP）に基づく運動方程式を採用することにより，ゲージ不変性とエーレンフェストの定理を満たすことができる（要件 iii）．本稿では TDVP を第 2 節で，TDMCSCF を第 3 節で紹介する．
式(2) （図 1）の CI 展開は，参照配置Φ(𝑡)(図1 の最初の配置)とΦ(𝑡)からの一電子励起配置，

盾であるという著しい利点を持つ．例えば式(4)で𝐶1% (𝑡)のみとっても，指数演算子を Taylor
展開するとわかるように，
1
𝑒 &'! = 1 + 𝐶1% + 𝐶1%% + ⋯2(6)積の形で多電子励起配置が取り込まれる．これにより電子数に依存しない計算精度と物理量の示量性が担保されるのが CC 法の強みである．私達は式(5)に基礎を置く時間依存最適化
結合クラスター法（TD-OCC）を開発している[11]．その特徴は，TD-MCSCF と同様に励起演
算子𝐶1+ (𝑡)と軌道関数をともに時間発展させる点である．TD-OCC は第 4 節で紹介する．
TD-MCSCF や TD-OCC の長所を活かし，高強度場現象を，実際に，精度良く計算するには高品質な計算プログラム開発が欠かせない．特に，軌道関数の柔軟で効率的な時空間離散化が鍵を握る．第 5 節では，実際の数値計算例を交えながらプログラム実装について詳しく解説する．
本稿で紹介する TD-MCSCF，TD-OCC 以外

二電子励起配置，…と分類して次のように書きにも，時間依存の多電子理論が多数提案されている．高強度場現象の理論と計算手法に関する直すことができる．
Ψ(𝑡) = -𝐶# (𝑡) + 𝐶1 (𝑡)/Φ(𝑡)(3)

日本語の入門的な解説が[12–14]に，より専門

𝐶1 (𝑡) = 𝐶1$ (𝑡) + 𝐶1% (𝑡) + ⋯(4)

ここで𝐶# (𝑡)は参照配置の複素振幅，𝐶1$ (𝑡)は一電子励起演算子，𝐶1% (𝑡)は二電子励起演算子…

的なレビューが[15]にある． [16]では時間依存
電子状態理論が総括的にレビューされている．
また[17]ではフェルミ粒子・ボーズ粒子の混合である．N 電子系で全 N 電子励起まで取り入

多体系の時間依存理論の概要と網羅的な引用文献が与えられている．以下断らない限りれる MCTDHF 法は，高精度だが計算コストがHartree 原子単位系を用いる．

電子数に対して指数関数的に増加するという

原子衝突学会誌しょうとつ 第19巻第3号 (2022)37

2．時間依存変分原理の時間微分𝑑〈𝑂8〉/𝑑𝑡に対してエーレンフェスト

2.1 ゲージ不変性とエーレンフェストの定理
採用し，古典的な光と相互作用する非相対論的

の定理
8
𝑑〈𝑂8〉
) P<ΨN + QΨR 𝜕𝑂 RΨS= −𝑖MΨ<O𝑂8, 𝐻
𝜕𝑡𝑑𝑡

電子からなる原子・分子を考える．系の時間発が成立する．電気双極子近似を超える扱いにつ展は式(1)の TDSE で厳密に与えられ，ハミル

いては，たとえば，私達の最近の研究[18]とそ

トニアンはこでの引用を参照されたい．

本解説ではボルン・オッペンハイマー近似を) (𝑡) = ℎ8(𝑡) + 𝐻
)%
𝐻

(15)

(7)
2.2 時間依存変分原理-ℎ8(𝑡) = + ℎ(𝒓, , 𝒑, , 𝑡),(8),

厳密な TDSE と異なり，近似理論では一般にゲージ不変性が満たされず，エーレンフェスト-)% = +
𝐻
,/.

1

(9)

<𝒓, − 𝒓. <である．ここで𝑁は電子数，𝑡は時刻，𝒓, , 𝒑, はそれぞれ𝑖 番目の電子の位置と正準運動量，ℎ(𝒓, 𝒑, 𝑡)は光との相互作用を含む一電子ハミルトニアンである．外部電磁場との相互作用で電気双極子近似を採用すると，長さゲージ(LG)の定理も成立しない．これらを満たす理論的枠組みを提供するのが時間依存変分原理(TDVP) である[19]．
TDVP では次の action 汎関数𝑆を導入し，)") (𝑡) − 𝑖 𝜕 RΨ(𝑡)S 𝑑𝑡
𝑆 = U QΨ(𝑡)R𝐻
(16)
𝜕𝑡
)
#

境 界 条 件 δΨ(𝑡# ) = δΨ(𝑡$ ) = 0 を 満 た す 変 分

および速さゲージ（VG）において，ハミルト)01 = ℎ801 + 𝐻
)% ，𝐻
)21 =
ニアンはそれぞれ𝐻δΨ(𝑡)に対して𝑆を停留化して次式を得る．
) − 𝑖 𝜕 [ RΨS = 0(17)
Re QδΨR Z𝐻
𝜕𝑡)% ，ℎ821 + 𝐻

式(17)で任意の変分δΨを許せば TDSE が再現ℎ801 (𝒓, 𝒑, 𝑡) =|𝒑|%

+ 𝑬(𝑡) ∙ 𝒓 + 𝑉345 (𝒓)2|𝒑 + 𝑨(𝑡)|%ℎ821 (𝒓, 𝒑, 𝑡) =+ 𝑉345 (𝒓)2(10)

されるが，近似理論を構築するにはパラメータ(11)Ψ(𝜶(𝑡))を考え，与えられたパラメータ空間内で与えられる．ここで𝑨(𝑡)はベクトルポテンシ
ャル，𝑬(𝑡) = −𝑑𝑨/𝑑𝑡は電場，𝑉345 (𝒓)は原子核

で𝑆を停留化する．すなわち，δΨ )𝜕
Re + δ𝛼.∗ Q
R Z𝐻 − 𝑖 [ RΨS = 0
δ𝛼.
𝜕𝑡

(18)

.

からのポテンシャルである．
長さゲージのハミルトニアンと速さゲージのハミルトニアンはゲージ変換
𝜕𝜕
)21 − 𝑖 = 𝑒 6,𝑨∙𝒓: F𝐻)01 − 𝑖 G 𝑒 ,;∙𝒓:𝐻𝜕𝑡𝜕𝑡𝜶(𝑡) = (𝛼$ (𝑡), 𝛼% (𝑡), ⋯ ) に 依 存 す る 波 動 関 数

より𝜶(𝑡)の運動方程式を導出できる．
上で述べた TDVP の他に McLachlan の変分(12)（𝒓H = ∑で関係づけられる．厳密な TDSE,<$ 𝒓, ）の解では，長さゲージの波動関数Ψ01 (𝑡)と速さゲージの波動関数Ψ21 (𝑡)はユニタリ変換
Ψ21 (𝑡) = 𝑒 6,;())∙𝒓: Ψ01 (𝑡)(13)

原理 [20]，Dirac-Frenkel の変分原理[21]と呼ばれる定式化もある．TDVP の特徴は，時間に依存しないハミルトニアンに対してエネルギーが保存し，ゲージ不変性とエーレンフェストの定理を満足する近似理論を構築できる点である．一方， McLachlan 形式にも特有のメリッで結ばれ，観測量の期待値はどちらのゲージでトがあり，量子コンピュータを用いた変分的時計算しても同じになる（ゲージ不変性）．またTDSE では，観測量𝑂8の期待値〈𝑂8〉 = MΨ(𝑡)<𝑂8<Ψ(𝑡)N(14)

間発展の観点からも注目されている[22]．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)

38

3. 時間依存自己無撞着場法

次式のようにアクティブ軌道{𝜙$ 𝜙% ⋯ 𝜙J }を 2

3.1 概要グループに分類し，TD-MCSCF においてただひとつの行列式を用いるのが TDHF 法である．閉殻系で空間軌(22)

各グループの占有数𝑁$ と𝑁% に次のような制限

道を-𝜙" /と書くと，TDHF 波動関数は%
Ψ>?@A (𝑡)：𝜙$% 𝜙%% ⋯ 𝜙-/%ΨKLCGH ：{𝜙$ 𝜙% 𝜙M }-" {𝜙N 𝜙O 𝜙J }-!

(19)と表記できるだろう．𝜙,% はスピンを対にした二電子に占有された空間軌道関数𝜙, をシンボを設けることにより4 ≤ 𝑁$ ≤ 6
0 ≤ 𝑁% ≤ 2

(23)
(24)

参照配置𝜙$% 𝜙%% 𝜙M% からの一電子励起および二電リックに示している．一方，与えられた𝑛 ≥子励起を考慮した CI 波動関数（量子化学で𝑁/2個の空間軌道で全𝑁電子を完全相関させる

CISD と呼ばれる）に基づく TD-MCSCF を構築

方法が MCTDHF であり，できる．

ΨCD>?@A (𝑡)：{𝜙$ 𝜙% ⋯ 𝜙E }-(20)と書ける．MCTDHF は当時東北大の加藤・河3.2 運動方程式
TD-MCSCF の action 汎関数は

野[4]と，当時ウィーン工科大の Scrinzi グルー)"

𝑆 = U f𝐿(𝑡) − 𝑖 + 𝐶!∗ 𝐶!̇ i 𝑑𝑡

プ[5]によって同時期に独立に開発された．与)#えられた数（𝑛）の軌道を用いた変分的に最善の解を保証するが，計算コストが𝑁および𝑛にと書くことができ，ラグランジアン𝐿は) − 𝑖𝑋8l<ΦP N𝐶P
𝐿 = + 𝐶!∗ MΦ! <j𝐻

対して指数関数的に増大する．

"

"Q

Q

1
"S QR
+ ℎQR 𝜌"S2

ℎQ" = U 𝜓"∗ (𝑥$ )ℎ𝜓Q (𝑥$ )𝑑𝑥$(21)

コア軌道はさらに時間発展する動的コア軌道と時間発展しない凍結コア軌道に柔軟に分類できる[6, 23]．式(21)のコア軌道とアクティブ
軌道の分割を柔軟に調整することにより，計算コストを削減したり，実験では得難い物理的洞察を得たりすることができる．
（具体例は 5.3 節
参照.）
より大きな系を扱うために，アクティブ軌道
空間内の打ち切った CI 展開に基づく手法も開発されている[7–9]．TD-ORMAS[9]では，アクティブ軌道を任意の数の軌道グループに分類
し，各グループの占有数に制限を課して柔軟な配置選択を可能にする．例えば 6 電子閉殻系で原子衝突学会誌しょうとつ 第19巻第3号 (2022)

(27)

"QSR

と二通りの等価な表現で与えられる．ここで内で完全相関させる．
-%
ΨDGHHDA (𝑡)：𝜙$% ⋯ 𝜙E%$ -𝜙E$I$ ⋯ 𝜙E /"𝐿 = +jℎQ − 𝑖𝑋Q l𝜌" +𝑛 − 𝑛F 個のアクティブ軌道に分類し，𝑁; = 𝑁 −2𝑛F 個のアクティブ電子をアクティブ軌道空間

(26)

!P

そこで TD-CASSCF[6]では，𝑛個の空間軌道
を常に二重占有される𝑛F 個のコア軌道と𝑛; =

(25)

!

(28)

"Q

ℎSR =
𝜓"∗ (𝑥$ )𝜓Q∗ (𝑥% )𝜓S (𝑥$ )𝜓R (𝑥% )(29)
UU
𝑑𝑥$ 𝑑𝑥%|𝒓$ − 𝒓% |（𝑥, = {𝒓, , 𝑠, }は空間𝒓, とスピン 𝑠, についての
電子座標）は一体および二体のハミルトニアン
Q

QR

行列要素，𝜌" ，𝜌"S は-𝜓" /を基底とする一体おQ
よび二体の縮約密度行列，𝑋8および𝑋" は軌道関数に対する時間微分演算子とその表現行列
𝜕𝜓"
𝜕𝜓Q
"
(30)
= 𝑋8𝜓" , 𝑋Q = Q𝜓" R 𝜕𝑡 S
𝜕𝑡である．𝐿の表現に式(26)を採用し，CI 係数について TDVP を適用すると𝑖

𝜕𝐶!
) − 𝑖𝑋8l<ΦP N 𝐶P
= +MΦ! <j𝐻
𝜕𝑡

(31)

P

39

を得る．これが CI 係数の運動方程式である．
𝐿の表現に式(27)を採用し，軌道関数について
TDVP を適用して次式を得る．
𝜕𝜓") l𝜓" + 𝑖 + 𝜓Q 𝑋"Q= 𝑄8jℎ8 + 𝑊
𝜕𝑡) = ℎ8 + 𝐻)% と𝑖𝑋8 = ℎ8 +みに他ならない．さらに𝐻
𝑖𝑋8% を式(39)に代入すると)TUU = ℎ8 + Π)𝐻)% + ij1 − Π) l𝑋8%
𝐻
(41)

(32)

) と比較すると，ゲージ依存のを得る．これを𝐻ℎ8は厳密に扱われていることがわかる．従って

これが軌道関数の運動方程式である．ここでTD-MCSCF はゲージ不変である．
次に一電子演算子𝑂8（𝜕𝑂8/𝜕𝑡 = 0とする）の期

𝑖

Q

QR 6$ )) 𝜓" ≡ + 𝑊RS 𝜓Q 𝜌)S
(𝜌 )"
𝑊

(33)

QSR)

𝑊RS (𝑥$ ) = U

𝜓S∗ (𝑥% )𝜓R (𝑥% )𝑑𝑥%|𝒓$ − 𝒓% |待値とその時間微分を考えよう．密度演算子"

𝜌H$ = +<𝜓" N𝜌Q ⟨𝜓Q |

(34)

(42)

"Q

は電子間相互作用の軌道関数ダイナミクスへの寄与，𝑄8 = 18 − ∑Q<𝜓Q NM𝜓Q <は一電子ヒルベルを導入すると，𝑂8の期待値とその時間微分はそ

ト空間において軌道関数が張る空間の補空間(43)への射影演算子である．軌道時間微分行列は
𝑖𝑋8 = ℎ8 + 𝑖𝑋8%""

𝑖 +-(𝑋% )S 𝜌QS − 𝜌S (𝑋% )SQ / + 𝚤̇
S

(35)
𝜕𝜌Q"
𝜕𝑡
(36)

"R SR
"R SR
= +jℎS) 𝜌Q)− 𝜌S) ℎQ)l
SR)

れぞれ𝜌H$ および𝜕𝜌H$ /𝜕𝑡によるトレース〈𝑂8〉 = tr𝑂8𝜌H$
𝑑〈𝑂8〉
𝜕𝜌H$
= tr𝑂8𝑑𝑡𝜕𝑡で与えられ，
𝜕𝜌H$
=
𝜕𝑡
"
𝜕𝜌Q
<𝜓" N z(𝑋𝜌)"Q − (𝜌𝑋)"Q +{ ⟨𝜓Q |
𝜕𝑡

(44)

(45)

である．ここで TDVP に基づく式(32)，(35)，で与えられる．

(36)が成立していれば
3.3 変分的時間発展の重要性
TD-MCSCF 波動関数[式(2)]を直接時間微分

(37)

!

に CI 係数および軌道関数の運動方程式[式(31)，(32)]を代入すると次式を得る．
𝜕Ψ)TUU Ψ𝑖=𝐻
𝜕𝑡)TUU = Π)𝐻) + 𝑖j1 − Π
) l𝑋8
𝐻

(46)

となり，これを式(44)に代入して
𝑑〈𝑂8〉) P<ΨN(47)= −𝑖MΨ<O𝑂8, 𝐻
𝑑𝑡を示すことができる．よってエーレンフェストの定理が成立する．

(38)

このように，CI 係数と軌道の変分的時間発(39)

展によってゲージ不変性とエーレンフェストの定理が保証される．ゲージ不変性によって初

ここで) (𝑡) = +|Φ! (𝑡)⟩⟨Φ! (𝑡)|ΠSR)−𝑖-(ℎ𝜌)"Q − (𝜌ℎ)"Q / − -(𝑋𝜌)"Q − (𝜌𝑋)"Q /して得られる次式𝜕Ψ𝜕𝐶!
=+
Φ (𝑡) + 𝑋8Ψ
𝜕𝑡
𝜕𝑡 !

𝜕𝜌Q"
"R SR
"R SR
= −𝑖 +jℎS)𝜌Q) − 𝜌S)
ℎQ) l
𝜕𝑡

(40)

!

めて外部電磁場存在下での物理量が意味を持
ち，エーレンフェストの定理によって物理量の) はその補空間は与えられた CI 空間への，1 − Π)TUU は TD-MCSCF 有への射影演算子である．𝐻

計算が著しく簡略化される．

効ハミルトニアンと解釈でき，その第一項は与4. 時間依存結合クラスター法

えられた CI 空間内のダイナミクスを，第二項

4.1 背景は軌道の時間発展による CI 空間自体のダイナ

一般に時間依存結合クラスター法では，CIミクスを表す．後者の存在が TD-MCSCF の強展開の代わりにクラスター展開(5)を用いる[11,原子衝突学会誌しょうとつ 第19巻第3号 (2022)40

24–26]．慣例に従って励起演算子を記号𝑇8で書き直しておこう．

4.2 時間依存軌道最適化結合クラスター法
一方我々は，よりシンプルな時間依存最適化

VW ())

Ψ(𝑡) = 𝑒 Φ(𝑡)𝑇8(𝑡) = 𝑇8$ (𝑡) + 𝑇8% (𝑡) + ⋯(48)

結合クラスター法（TD-OCC）を定式化した[11]．

(49)

𝑇8$ = 𝑇,X 𝑎HXY 𝑎H, ≡ 𝑇,X 𝐸8,X

(50)

TD-OCC は次の実値 action 汎関数)")𝑒 6VW Z𝐻) − 𝑖 𝜕 [ 𝑒 VW RΦS 𝑑𝑡
𝑆 = Re U QΦRΛ𝜕𝑡)

𝑇8% = 𝑇,.XZ 𝑎HXY 𝑎HZY 𝑎H. 𝑎H, ≡ 𝑇,.XZ 𝐸8,.XZ

(51)

ここで

は生成 (消滅) 演算子𝑎H"Y (𝑎H" )で書かれた一電子励起演算子および二電子励起演算子である．
式(48)を用いる結合クラスター法は定常状態の量子化学で非常によく用いられており，数十原子までからなる分子の高精度計算は結合クラスター法の独壇場となっている．ところが
時間依存の軌道関数を用いた時間依存結合クラスター法は最近まで開発されていなかった．
その最大の理由は，式(49)を打ち切ったとして&)𝑒 VW の展開が閉じないため，式(16)のも𝑒 VW 𝐻action 汎関数の計算に MCTDHF と同じ指数関数コストを要することである．一方，相似変換) 𝑒 VW の展開は有限項で閉
ハミルトニアン𝑒 6VW 𝐻
じ，1) 𝑒 VW = 𝐻) + O𝐻) , 𝑇8P + €O𝐻) , 𝑇8P, 𝑇8•
𝑒 6VW 𝐻
2!
1
) , 𝑇8P, 𝑇8• , 𝑇8„+ ƒ€O𝐻
3!
1) , 𝑇8P, 𝑇8• , 𝑇8„ , 𝑇8†+ …ƒ€O𝐻
4!

𝑇8 = 𝑇8% + 𝑇8M + ⋯ = + 𝑇!; 𝐸8!;

(55)

!,;

)=1+Λ
)% + Λ
) M + ⋯ = 1 + + Λ!; 𝐸8;!
Λ

(56)

!,;

で与えられる励起および脱励起演算子である．
XZF⋯=
た だ し 多 体 テ ン ソ ル の 添 字 を 𝑇,.+⋯,.+⋯

𝑇!; , ΛXZF⋯ = Λ!; のように略記した．式(55)，(56)で一電子励起・脱励起を考慮せず，代わりに軌
道関数を最適化する（変分的に時間発展させる）のが本手法のポイントであり名の由来である．
Action 汎関数(54)は次のように変形でき，
)"
𝑖
̇ ;∗
𝑆 = U f𝐿 − +jΛ!; 𝑇̇!; − Λ!∗; 𝑇! li 𝑑𝑡
2
)#

(57)

!;

TD-MCSCF の場合と同様，ラグランジアン𝐿は

"

"

Q

𝐿 = +jℎQ − 𝑖𝑋Q l𝜌" +
"Q

1
"S QR
+ ℎQR 𝜌"S
2

(58)
(59)

"QSR

𝐿の表現に式(58)を用い，式(57)をΛ!; と𝑇!; につ

量子化学ではこの定式化が用いられている．
そこで Kvaal[26]は次の action 汎関数に基づく擬変分的手法を提案した．
)")𝑒 6VW Z𝐻) − 𝑖 𝜕 [ 𝑒 VW RΦS 𝑑𝑡
𝑆 = U QΦRΛ𝜕𝑡)

)はに基礎を置く．ここで𝑇8およびΛ

等価な二通りの表現で与えられる．
1)𝑒 6VW j𝐻) − 𝑖𝑋8l𝑒 VW <ΦN + 𝑐. 𝑐
𝐿 = MΦ<Λ
2
(52)

(54)

#

(53)

#

いて停留化して次式を得る．
𝜕) − 𝑖𝑋8l𝑒 VW <ΦN
𝑖 𝑇!; = MΦ!; <𝑒 6VW j𝐻𝜕𝑡𝜕 !
)𝑒 6VW Oj𝐻) − 𝑖𝑋8l, 𝐸8!; P𝑒 VW |Φ⟩Λ = 𝑖⟨Φ|Λ
𝜕𝑡 ;

(60)
(61)

)は脱励起演算子である．Kvaal の手法では式Λ

これらが𝑇!; とΛ!; の運動方程式である．一方，縮(53)を解析的複素関数として扱う．数理的に興約密度行列を用いた表式(59)は，TD-MCSCF の味深い理論構造を持つが，非直交な軌道関数，式(27)と同じである．従って軌道の運動方程式

非エルミートな密度行列，虚値の物理量などがは式(32)で与えられる．実値 action 汎関数(54)現れて扱い難い．

を基礎とすることにより，エルミートな密度行列が自然に導かれ，TD-MCSCF との高いアナロジーが得られる．さらに 3.3 節と同様の議論
によりゲージ不変性と実値物理量に関するエーレンフェストの定理が成立する．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)41

式(55)，(56)で全ての二電子項を取り込む手

5. 実装の概略と数値計算例

法が TD-OCCD，さらに全ての三電子項を取りTD-MCSCF や TD-OCC のプログラム実装で込む手法が TD-OCCDT である．アクティブ電は，軌道関数が与えられたときその時間微分を子数𝑁; に対して TD-OCCD の計算コストは

返す軌道 module と，CI 係数や CC 係数が与え

𝑂(𝑁;J )，TD-OCCDT の計算コストは𝑂(𝑁;] )の依られたときその時間微分を返す配置 module を存性を持つ．

統合する．両者間でハミルトニアン行列要素と4.3 時間依存軌道最適化二次摂動法

縮約密度行列要素を通信する必要がある（図前節で導入した TD-OCCD を出発点とする

2）．CI や CC の演算は量子化学プログラムに

よ り 低 コ ス ト の 計 算 手 法 を 考 え よ う ． TD-

倣って容易に実装でき，焦点は軌道関数の運動

OCCD について式(60)は) % lj𝐻) − 𝑖𝑋8l<ΦN
2𝐿DD? = MΦ<j1 + Λ) % lO𝐻) − 𝑖𝑋8, 𝑇8% P<ΦN
+ MΦ<j1 + Λ
1)% €O𝐻) − 𝑖𝑋8, 𝑇8% P, 𝑇8% • ŠΦS + 𝑐. 𝑐
+ QΦŠΛ
2方程式の実装となる．我々はイオン化過程を記述したいので，Gauss 型関数など局在化した基(62)

ッドで直接離散化している[23, 28]．原子核近と書き下せるが，Møller-Plesset 摂動論に倣って) − 𝑖𝑋8 = j𝑓1 − 𝑖𝑋8l + 𝑉8
𝐻
(63)
"
".
".
"
𝑓1 = ŒℎQ + +jℎQ. − ℎ.Q l• 𝐸8Q

底による展開は用いず，軌道関数を実空間グリ(64)

傍で必要な高解像度と高強度レーザー中の大振幅電子運動をカバーできる広い計算領域の両立が中心的な課題である．
原子と直線偏光レーザーとの相互作用に特化したコード[23]では，CI 係数[式(2)]または.

) % について二次までと分割し，𝐿DD? で𝑉8 , 𝑇8% , ΛCC 係数[式(55)，(56)]の保持に必要なメモリサの項を残すと次式を得る．
) − 𝑖𝑋8l<ΦN
𝐿C^% = MΦ<(1 + 𝑥)j𝐻
) − 𝑖𝑋8, 𝑇8% P<ΦN
+MΦ<O𝐻1)+MΦ<Λ% O𝑓 − 𝑖𝑋8, 𝑇8% P<ΦN + 𝑐. 𝑐イズが計算の可否を決める．大雑把には，研究室レベルのワークステーションで TD-OMP2(65)

この時間依存 Hylleraas 汎関数から出発して，基底状態の軌道最適二次摂動法（OMP2）に対応する時間依存軌道最適化二次摂動法（TD-OMP2）を導出できる[27]．TD-OMP2 はは𝑁; 〜1000，TD-OCCD は𝑁; 〜100，TD-OCCDT
は𝑁; 〜20 までの計算を実行できる．計算時間は軌道関数を離散化するグリッド点数𝑁_ に比
例し，照射するレーザー電場の強度𝐼# や波長𝜆#に強く依存するj𝑁_ ∝ “𝐼# 𝜆%# l．具体的に，Sn 原子の Ar 殻を凍結コア近似し，𝑁; ＝13 個の価電摂動論に基づくが，軌道が時間発展し外場は無摂動ハミルトニアン𝑓1に取り込まれているため

子を𝑛; =18 個のアクティブ軌道で相関させる

高強度場現象にも適用できる．実値変分原理よW/cm2, 中心波長𝜆# = 800 nm の foot-to-foot 8,.∗
り𝑇,.XZ = ΛXZ が導かれ，𝑇,.XZ の運動方程式は

𝑖

𝜕 XZ
XZ
XZ
𝑇 = 𝑢,.
− 𝑝,. 𝑓.+ 𝑇,++ 𝑝XZ 𝑓FX 𝑇,.FZ
𝜕𝑡 ,.

TD-CASSCF 計 算 で は ， 強 度 𝐼# = 2 × 10$N

時間発展アルゴリズム(66)（𝑝"Q は置換演算子）と著しく簡略化される．
TDHF 計算が可能な系では TD-OMP2 も可能であることが多く，大規模系の相関ダイナミクスに適用可能である．

ü 陽的Runge-Kutta法
ü Operator splitting法 [23]ü Exponential integrator [40]%

軌道module
ü 極座標 [23]
ü 直交座標 [28]
ü 曲線座標 [40]

!

"%/"$
"

"$

ℎ! , ℎ!#""!

(! , ($#
"
)!

"!/"$

配置moduleü CASSCF [6]ü ORMAS [9]ü OCC [11]図 2: TD-MCSCF/TD-OCC のプログラム実装の概
要．要素技術の詳細は[括弧内]に与えた文献を参照されたい．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)42

サイクルパルスを照射した場合に収束したと書くことにする．ここで右辺の各列は軌道関

HHG スペクトルを得るには動径グリッド点数

数(32)および CI 係数(31)または CC 係数(60)，450，球面調和関数展開の最大角運動量 47，等

(61)の運動方程式右辺から成る．𝑡, と 𝑡,I$ =間隔 80,000 時間格子点が必要で，クロック周

𝑡, + 𝛿𝑡をそれぞれ𝑖番目と𝑖 + 1番目の時間格子

波数 2.4 GHz の Intel Xeon 40 プロセッサを用い

点，𝛿𝑡を時間刻み幅とし，式(67)を区間[𝑡, , 𝑡,I$ ]た場合この計算におよそ 15 時間を要した．

で積分すると次式を得る．

原子・分子やレーザー電磁場に対称性を仮定
しない多重解像度グリッド[28]や曲線座標を)'("

𝒖(𝑡,I$ ) = 𝒖(𝑡, ) + U

𝒇j𝒖(𝑡)l 𝑑𝑡

)'

(68)

用いたコードではより大きなメモリ容量や長上式右辺で𝒇j𝒖(𝑡)l = 𝒇j𝒖(𝑡, )lと近似すると陽い計算時間を要するが，アラニン（重原子 1s 軌

的 Euler 法の時間発展公式が得られ，

道を擬ポテンシャル近似して 36 電子系）や同𝒖(𝑡,I$ ) = 𝒖(𝑡, ) + 𝛿𝑡𝒇j𝒖(𝑡, )l程度の分子量の分子の多電子ダイナミクス計

これは高次の陽的 Runge-Kutta 法に容易に拡張算が可能になりつつある．

できる．陽的 Runge-Kutta 法は実装が簡単で時(69)

間ステップ毎の計算コストが小さいのが長所
5.1 時間発展アルゴリズム

だが，高い空間解像度（細かいグリッド）や特

実空間で離散化された軌道関数と CI 係数ま異点のある座標系（例えば極座標）での計算でた は CC 係 数 を ま と め て 列 ベ ク ト ル 𝒖 =は極端に小さな時間刻み幅が必要になってしV(𝝍, 𝑪) とし，その運動方程式を形式的に
𝜕𝒖(67)= 𝒇(𝒖)
𝜕𝑡

まう．
そこで式(67)を次のように書き換える．
𝜕𝒖(70)= 𝐿𝒖 + 𝑾(𝒖)𝜕𝑡dipole acceleration (a.u)

dipole acceleration (a.u)

ここで右辺第一項は小さな時間刻み幅を要す
0.015
0.010

N=5000N=10000
N=20000
N=40000

るが計算が容易な線形項であり，右辺第二項は計算コストは高いが滑らかな非線形項である．

0.005
0.000
-0.005-0.010(a)0.0N=2000
0.5
0.015 N=3000
N=40000.010 N=5000

1.0
1.5
2.0
time (optical cycle)

2.5

3.0

0.000
-0.0050.0(71)

上式を区間[𝑡, , 𝑡,I$ ]で積分すると次式を得る．
𝒖(𝑡,I$ ) = 𝑒 a)𝑳 𝒖(𝑡, )
)'("
(72)
+ U 𝑒 ()'("6))𝑳 𝑾j𝒖(𝑡)l 𝑑𝑡

0.005-0.010式(70)は次式と等価であり，
𝜕 6)𝑳(𝑒 𝒖) = 𝑒 6)𝑳 𝑾(𝒖)
𝜕𝑡

)'

(b)

上式右辺で𝑾j𝒖(𝑡)l = 𝑾j𝒖(𝑡, )lと近似すると
0.5

1.0

1.5

2.0

2.5

3.0

time (optical cycle)

142

図 3: 波長 800 nm，強度 8 ×10 W/cm ，foot-tofoot 3 サイクルレーザーパルスに照射された Ne 原子の双極子加速度．1s 軌道を凍結コアとし，価電
子 8 電子を 13 軌道で相関させた TD-CASSCF による計算結果．(a) 二次の Operator splitting 法，(b) 四次の Exponential Runge-Kutta 法を用いた
場合について，レーザー電場 1 周期毎の時間ス
テップ数 N を変えた計算結果の比較．
Exponential Runge-Kutta 法 に よ っ て Operator
splitting 法よりも早く時間ステップについて収束した結果が得られていることがわかる．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)

Exponential Euler 法の時間発展公式が得られ，
𝒖(𝑡,I$ ) = 𝒖(𝑡, ) + 𝛿𝑡𝜙$ (𝑳)𝒇j𝒖(𝑡, )l(73)𝑒 a)𝑳 − 1(74)
𝛿𝑡𝑳
これは高次の Exponential Runge-Kutta 法に容易
𝜙$ (𝑳) =に拡張できる．Exponential Runge-Kutta 法は陽的 Runge-Kutta 法と比べて時間ステップ毎の計
算コストは大きいが，高い空間解像度や特異点のある座標系での計算でも reasonable な時間刻43

である．図 3 に高強度レーザー照射下の Ne 原子の加速度応答の計算結果を示す．よく用いられる二次の Operator splitting 法と比較して，Exponential Runge-Kutta 法はより大きな時間刻みで安定かつ正確な計算が可能であることがRadial distribution function

み幅で安定に時間発展できる点が大きな利点!!"#

6x10-55x10

4x10-5mask
irECS
large box

3x10-5
-5

2x10

-5

1x10

0
0

20

見てとれる．式(70)における線形項と非線形項の分割は任意だが，私達の経験では，軌道関数の時間発展に係る一電子ハミルトニアンを線形項とし (𝐿 = −𝑖ℎ) ，残り全てを非線形項とするのが最も効率的である．
5.2 吸収境界条件
イオン化を伴うシミュレーションでは，計算領域（box）の端に到達した軌道の反射を防ぐ

!$

-5

40
60
80 100 120
Radial coordinate (a.u)

140

図 4: 波長 800 nm，強度 8 ×1014 W/cm2 の footto-foot 5 サイクルレーザーパルス照射後の He 原子の動径波動関数（5 軌道 MCTDHF）．半径𝑅! = 64 か ら マ ス ク 境 界 条 件 を 用 い た 場 合（mask），同じスケーリング半径𝑅! = 64 の irECS
境界条件を用いた場合（irECS），および反射が無視できる大きな box を用いた場合（large box）の結果の比較．非スケーリング領域では irECS と
large box の結果がよく一致している．𝑅"#$ (≤ 𝑅! )は本文で定義されているイオン化半径の設定例．

ための吸収境界条件が不可欠である．このために，ある半径𝑅# から外側で 1 から 0 に漸減する関数を時間ステップ毎に軌道に乗じるマスQ

U 𝑑𝑥𝜑Q∗ (𝑥)𝜑" (𝑥) = 𝛿"(77)

ク関数法や，複素吸収ポテンシャル法がよく用

また非スケーリング領域（|𝑟| < 𝑅# ）で正確ないられるが，我々はより効率的な吸収境界条件（全空間で解いた場合の結果と必要な精度でである外部複素スケーリング（ECS）法の多電一致する⇔𝑅# について収束した）波動関数を，子系への適用に初めて成功した[29]．簡単のた

驚くほど小さな𝑅# で得られる．具体例として高め極座標を考えると，多電子 ECS では
𝑟 → 𝑅(𝑟)𝑟(𝑟 < 𝑅# )
(75)
=¢
,b
𝑅# + (𝑟 − 𝑅# )𝑒(𝑟 ≥ 𝑅# )強度レーザー照射直後の He 原子における電子のようにスケーリング半径𝑅# の外側で動径座
標を角度𝜂だけ複素平面に解析接続し，軌道の運動方程式(32)を式(75)によって変換する．さらに，スケーリング領域を無限遠まで拡張し，[𝑅# , ∞)で指数関数的に減衰する離散値表現基
底を用いる無限範囲 ECS（irECS）により，原理的には全空間をカバーするシミュレーショ

密度の動径分布関数を図 4 に示す．irECS とマ
スク関数法の結果を比較すると，irECS が非ス
ケーリング領域で large box の結果をより良く再現できていることがわかる．
5.3 高次高調波スペクトルの計算
エーレンフェストの定理と正準交換関係を
用 い る と ， 電 子 座 標 𝑧 成 分 の 期 待 値 〈𝑧̂ 〉 =∑"Q 𝜌Q" M𝜓Q <𝑧<𝜓" Nの二階微分である加速度の期

我々の irECS の実装では，イオン化とともに

待値 〈𝑎H〉 = 𝑑% 〈𝑧̂ 〉/𝑑𝑡 % は次式で与えられる．
𝜕𝑉〈𝑎H〉 = + 𝜌Q" Q𝜓Q R Z− 3 [ R𝜓" S(78)
𝜕𝑧

次式の計算 box 内規格直交性は破られ得るが，高強度フェムト秒レーザー照射下の原子・分子ンが可能となる．

(𝑆cde )Q" = U|S|gh#

𝑑𝑥 𝜑Q∗ (𝑥)𝜑" (𝑥) ≠ 𝛿"Q(76)

式(75)の逆変換を適用した軌道の全空間重な
り積分についての規格直交性が，時間発展の精度の範囲内で保存される．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)"Q

について〈𝑎H〉の時間 Fourier 変換から HHG スペ
クトルを計算できる．原理的には〈𝑧̂ 〉 からもHHG スペクトルを得られるが，遠方で発散す
る位置演算子の直接評価は数値的困難を伴う．
エーレンフェストの定理によってこれが局在化した核引力−𝜕𝑉3 /𝜕𝑧に置き換えられ，5.2 節

44図 5: (a) 波長 800 nm，強度 1015 W/cm2，foot-to-foot 2 サイクルレーザーパルスに照射された Ne 原子からの高次高調波スペクトル．時間依存密度汎関数法（TDDFT），TDHF，および TD-CASSCF の計算結果の比較．TD-CASSCF では 1s 軌道を凍結コアとし，価電子 8 電子の記述に 8 軌道[(1,0,8)とラベル]または 13 軌道[(1,0,13)とラベル]を用いた．(b) TD-ORMAS 法による波長 770 nm，強度 3×1014 W/cm2，footto-foot 4 サイクルレーザーパルスに照射された Mn 原子からの高次高調波スペクトル． Ne コアと 3s 軌道の 12 電子は凍結コアとして扱い，初期キャラクター3p, 3d, 4s の 13 電子をアクティブに計算した場合（active 3p とラベル）と，3p 軌道まで凍結コアとし初期キャラクター3d, 4s の 7 電子をアクティブに計算した場合（frozen 3p とラベル）の結果の比較．(c) (b)で用いたアクティブ軌道の模式図．図示された参照配置
からの一電子励起および二電子励起までを考慮に入れた．

の吸収境界条件によって正確な HHG スペクトで与えられる．ただし{𝑥, , 𝑥,I$ , ⋯ , 𝑥- }をまとめルの計算が可能になる．

て𝑦, と書き，d𝑦, = 𝑑𝑥, 𝑑𝑥,I$ ⋯ 𝑑𝑥- とした．一般

図 5 に Ne 原子と Mn 原子の HHG スペクトルの計算例を示す．Ne の例では価電子軌道の8 電子をアクティブとし，アクティブ軌道数𝑛;を𝑛; = 8，𝑛; = 13と変えた計算がよく一致すに𝑝電子密度分布関数は次式で与えられる．
𝜌" j𝑥$ , ⋯ , 𝑥" , 𝑡l
𝑁!
(80)
=
U 𝑑𝑦"I$ |Ψ(𝑦$ , 𝑡)|%(𝑁 − 𝑝)!

ることから計算結果の収束を確認することがレーザー場中の原子・分子のイオン化を定量できる．Mn の場合，初期キャラクター3p の軌化するには，イオン化半径𝑅id3（|𝒓| ≥ 𝑅id3 で近道関数をアクティブとした場合と凍結コアとした場

似的に自由電子と見做せる半径）を導入し，次合の計算結果の比較から，3p-3d 巨大共鳴が光式で𝑘電子イオン化確率を定義すると便利で子エネルギー50 eV 付近の高調波強度増大をもある．

たらすことを見出した[30]．これらの例が示すように，系統的かつ柔軟な近似理論の適切な実装

𝑃+ = U 𝑑𝑥$ ⋯ 𝑑𝑥+ U 𝑑𝑦+I$ |Ψ(𝑦$ )|%
/

(81)

g

により，計算結果の信頼性を確実にしたり，実

ここで> (<) を付した積分では空間積分をイ

験では得難い物理的洞察を得たりすることがオン化半径の外側|𝒓| ≥ 𝑅id3 （内側|𝒓| < 𝑅id3 ）できる．

に制限する．式(81)を直接計算できるのは少数

5.4 イオン化過程の記述

電子系に限られるが，形式的にTD-MCSCF や TD-OCC を解いて得られる時間依存の多電子波動関数を直接解析するのは多くの場合容易ではなく，注目する電子以外の情報を縮約すると見通しが良くなる．例えば一電子密度分布関数は-U 𝑑𝑦+I$ = ¯ °U 𝑑𝑥, − U 𝑑𝑥, ±
g

,<+I$

として二項定理を用いると，-

𝑃+ (𝑡) = +j<+𝜌$ (𝑥$ , 𝑡) = 𝑁 U 𝑑𝑥% ⋯ |Ψ(𝑥$ , 𝑥% , ⋯ , 𝑡)|%(79)
≡ 𝑁 U 𝑑𝑦% |Ψ(𝑦$ , 𝑡)|%

原子衝突学会誌しょうとつ 第19巻第3号 (2022)

(82)

/

(−1)j6+𝑇 (𝑡)𝑘! (𝑚 − 𝑘)! j

𝑇j (𝑡) = U 𝑑𝑥$ ⋯ 𝑑𝑥j 𝜌j (𝑥$ , ⋯ , 𝑥j , 𝑡)

(83)

(84)

/

45

0.5

(83)の展開を打ち切ることで効率的な計算が0.4可能である．例えば，式(83), (84)で多電子密度を軌道表示に変換し，𝑃% 以降，𝑃M 以降，𝑃N 以降

probabilityを得る．この多電子密度による展開では，式TDHF
TD-OCCD
TD-OCCDT
TD-CASSCF

0.3
0.2

が無視できる場合の𝑃$ をそれぞれ𝑃$k , 𝑃$kk , 𝑃$kkk

0.1

とすると

0.0
0.041.0

(85)

1 ""
Q
Q
𝑃$kk = 𝑃$k − 𝜌Q""Q!! (𝑆/ )""" (𝑆/ )"!!
2

(86)

" " "
Q
Q
Q
𝑃$kkk = 𝑃$kk + 𝜌Q""Q!!Q)) (𝑆/ )""" (𝑆/ )"!! (𝑆/ )"))

(87)

1.5
2.0
2.5
TDHF
TD-OCCDtime (optical cycle)
TD-OCCDT
0.03 TD-CASSCF
probability"
Q
𝑃$k = 𝜌Q"" (𝑆/ )"""交条件(77)より次式が成立し，(𝑆/ )Q" = 𝛿"Q − U|𝒓|gh*+,

𝑑𝑥 𝜑Q∗ (𝑥)𝜑" (𝑥)(88)|𝒓| < 𝑅id3 の軌道のみから(𝑆/ )Q" を得ることができる．図 4 のように𝑅id3 ≤ 𝑅# とし，5.2 節の

3.0

0.02
0.01

(b) Ar2+が得られる．ここで(𝑆/ )Q" はイオン化半径の外側|𝒓| ≥ 𝑅id3 でとった重なり行列だが，規格直(a) Ar+

0.00
1.0

1.5

2.0
2.5
time (optical cycle)
143.0

図 6: 波長 800 nm，強度 6 ×10 W/cm2 の footto-foot 3 サイクルレーザーパルス照射時間の Ar
原子からの一重電離確率(a)および二重電離確
率 (b) ． TD-OCCD ， TD-OCCDT ， お よ び TDCASSCF では Ne コアを凍結コアとし，価電子軌
道の 8 電子を 13 個のアクティブ軌道内で相関させた．TDHF によるアセスメントの後，式(86)（および𝑃% の対応する式）を用いて計算した．

吸収境界条件を適用することで有限の box サ
イズで式(88)を精度よく評価できる．
どれだけの電子がイオン化するかは実験や= 𝑁 U 𝑑𝑦% Ψ(𝑥$ , 𝑦% , 𝑡)Ψ∗ (𝑥$k , 𝑦% , 𝑡)

シミュレーションを行うまでわからないので，が実空間表示の一電子縮約密度行列である．こ一連の公式(85)–(87)を評価して収束を観る必

れを用いて光電子密度行列を次式のように導要がある．幸いこの収束性は計算手法に強く依存しないので，精度とコストの低い手法（たと入しよう．
𝜌^l (𝑥; 𝑥 k ) = 𝜃id3 (𝒓)𝜃id3 (𝒓k )𝜌$ (𝑥; 𝑥 k )えば粗いグリッドで TDHF）でアセスメントしここでた後に精度とコストの高い手法（たとえば細か

𝜃id3 (𝒓) = ¢いグリッドで TD-CASSCF）で計算することが(90)

0 (|𝒓| < 𝑅id3 )1 (|𝒓| ≥ 𝑅id3 )(91)

できる．具体例として，高強度レーザー照射中

はヘビサイド関数である．運動量空間の情報をの Ar 原子における一重電離および二重電離確

得るには光電子密度行列を Fourier 変換する．
𝜌^l (𝒌; 𝒌k , 𝑡)
(92)
= U 𝑑𝑥 U 𝑑𝑥 k 𝜒𝒌∗ (𝒓)𝜌^l (𝑥; 𝑥 k , 𝑡)𝜒𝒌- (𝒓k )率の時間変化の様子を図 6 に示す．式(85)–(87)を吟味して時空間離散化についての収束を確認できており，これによって初めて，電子相関効果（計算手法による結果の相違）の物理的なここで𝜒𝒌 (𝒓)は外場中の平面波（Volkov 波）で，議論が可能となる[11]．

し，非対角項も考慮した𝜌$ (𝑥$ ; 𝑥$k , 𝑡)

原子衝突学会誌しょうとつ 第19巻第3号 (2022)

Long
time

Pulse

0

(89)

!!"#

5

"

10

15

20

position

25

30

350

photoelectron

Re[wave function]

前節で議論した一電子密度分布関数を拡張!!"#
Re[wave function]

5.5 光電子スペクトルRe[wave function]

!!"#

405

10

"

15

20

25

position

30

35

40
0

5

10

15

"

20

25

30

35

40

position

図 7: tSURFF 法の模式図．

46

6

𝜒𝒌 (𝒓) = (2𝜋)M
% 𝑒 ,𝒌∙𝒓In())

(93)

𝑖 )
𝜙(𝑡) = − U 𝑑𝜏|𝒌 + 𝑨(𝑡)|%
2 6o

(94)

これを用いて，光電子エネルギースペクトル（PES）を求めることができる．

|𝒑 + 𝑨|%
𝜕
𝜒+ =
𝜒+
𝜕𝑡
2
|𝒑 + 𝐴|%𝜕Q𝑖 𝜑" ≈
𝜑" − 𝜑Q 𝑌"
𝜕𝑡
2

(95)

ここでΩは𝒌の方向角，𝑡は光との相互作用が完了して十分経過した後の時刻である．
TD-MCSCF や TD-OCC の計算から光電子ス
ペクトル（PES）を抽出するために，式(92)の縮約密度行列を軌道表示に変換すると(99)を用いて光電子振幅の運動方程式を導出し，𝑖𝜌^l (𝐸) = U 𝑑Ω |𝒌|% 𝜌^l (𝒌; 𝒌, 𝑡)

(98)

𝑖

𝜕
Q
𝑎 (𝒌, 𝑡) = − + 𝑎Q (𝒌, 𝑡)𝑌"𝜕𝑡 "Q

+ Q𝜒+ R ƒ𝜃id3 ,

(100)

%

|𝒑 + 𝑨|„ R𝜑" S2

これを数値的に時間積分する．ここで𝑌"Q =) l<𝜑" N − 𝑖𝑋"Q である．式(100)第一項
M𝜑Q <jℎ8 + 𝑊
は|𝒓| < 𝑅id3 の情報のみから計算でき，ヘビサ
イド関数の微分（𝒑との交換子）はデルタ関数

𝜌^l (𝒌; 𝒌k , 𝑡) = 𝑎" (𝒌, 𝑡)𝜌Q" 𝑎Q∗ (𝒌k , 𝑡)

(96)𝑎" (𝒌, 𝑡) = U 𝑑𝑥𝜒𝒌∗ (𝒓)-𝜃id3 (𝒓)𝜓" (𝑥, 𝑡)/(97)だから，式(100)第二項は境界面|𝒓| = 𝑅id3 での面積分に帰着する．図 4 のように𝑅id3 ≤ 𝑅# とし，5.2 節の吸収境界条件を適用することで有を得る．式(97)の光電子振幅を直接求めるには，
イオン化領域に入った電離波束𝜃id3 (𝒓)𝜓" (𝑥, 𝑡)を全て保持できる巨大な計算領域が必要となる（図 7）．我々はこれを避けるために，多電子系の時間依存 surface flux 法（tSURFF）を開発した[31]．tSURFF では，|𝒓| ≥ 𝑅id3 で成立する
次式

限の box サイズで PES を精度よく計算できる．
具体例として Be 原子からの光電子断面積の実験結果と計算結果の比較を図 8 に示す．

6. まとめと展望
本稿では，高強度レーザー場中の多電子ダイ
ナミクス計算について，計算手法の理論的背景を概観した後，計算プログラム実装について数値計
算例を交えながら解説した．後者は思い切って掘り下げた内容まで含めたが，数値計算現場の臨場感を読者にお伝えできていれば幸いである．

(a) Experiment

私達はここで紹介した計算手法以外にも，時間に依存しない軌道を用いるがゲージ不変性を満

cross section (arb. units)

0.008

37fs

0.007

75fs

112fs

たす gauge-invariant TDCIS [33, 34]，TD-OCCD

0.006に 対 す る 近 似 で あ る TD-OCEPA0[35], TD-

0.005

OCCDT に対する近似である TD-OCCDT(4)[36]，

0.004
0.003
0.002
0.001

多体波動関数を介さず二電子縮約密度行列を直(b) Calculated
9.5

10

10.5

11

11.5

12

12.5

13

photon energy (eV)

図 8: 光子エネルギー22 eV の広帯域パルスに照射された Be 原子の光電子断面積．(a)実験結果[32]と(b)TD-CASSCF 法による計算結果の比較．
13 eV 近傍の振動構造は 2pnl 自動電離状態を介したイオン化経路の干渉に由来する．パルス終了
後の時間経過（37 fs，75 fs，および 112 fs）とともに成長する干渉構造は自動電離の進行を反映
し，スペクトル形状は実験結果と良く一致している．

原子衝突学会誌しょうとつ 第19巻第3号 (2022)

接実時間発展させる TD-2RDM[37, 38]などを開発してきた．
本稿では電磁場と電子の相互作用に電気双極子近似を仮定したが，例えば，光学異性体識別のための二色円偏光高次高調波分光では磁気相互作用が主役を演じ，電気双極子近似を超えた定式化が必須となる[39]．重い元素や内殻電子を取り扱う場合にはスピン軌道相互作用などの相47

対論効果の取り込みが望ましいケースもあり[40]，多電子時間依存波動関数理論で相対論効果を取り込める枠組みの開発が必要である．一方，任意のフェルミ粒子とボゾン粒子の複合系についての基礎理論が成熟しつつあり[17, 41, 42]，高強度レーザーと分子との相互作用において電子だけ
でなく原子核も量子力学的に扱う，非断熱分子ダ
イナミクス計算の発展が期待されている．タンパク質の構成ユニットであるアミノ酸や同規模の分子について非断熱分子ダイナミクス計算が実現すれば，電子運動のアト秒制御とそれに基づく新物質科学，という夢の実現に大きく貢献できると思われる．

参考文献[1] M. Protopapas, C. H. Keitel, and P. L. Knight,
Rep. Prog. Phys. 60, 389 (1997).
[2] F. Krausz and M. Ivanov, Rev. Mod. Phys. 81,
163 (2009).
[3] 大森賢治 他，アト秒科学（化学同人，2015）．
[4] T. Kato and H. Kono, Chem. Phys. Lett. 392,
533 (2004).
[5] J. Caillat et al, Phys. Rev. A 71, 012712 (2005).
[6] T. Sato and K. L. Ishikawa, Phys. Rev. A 88,
023402 (2013).
[7] H. Miyagi and L. B. Madsen, Phys. Rev. A 87,
062511 (2013).
[8] D. J. Haxton and C. W. McCurdy, Phys. Rev. A
91, 012509 (2015).
[9] T. Sato and K. L. Ishikawa, Phys. Rev. A 91,
023417 (2015).
[10] T. Helgaker, P. Jørgensen, and J. Olsen,
Molecular Electronic-Structure Theory (Wiley,
2000).
[11] T. Sato et al, J. Chem. Phys. 148, 051101(2018).
[12] 石川顕一,佐藤健, レーザー研究, 44, 784(2016).
[13] 佐藤健,石川顕一,アンサンブル, 18, 236(2016).
[14] 石川顕一,佐藤健, 日本物理学会誌, 74, 40(2018).
[15] K. L. Ishikawa and T. Sato, IEEE J. Sel. Topics
Quantum Electron. 21, 8700916 (2015).
[16] X. Li et al, Chem. Rev. 120, 9951 (2020).
[17] A. U. J. Lode et al, Rev. Mod. Phys. 92,
011001 (2020).
[18] R. Anzaki et al, Phys. Rev. A 98, 063410(2018).
[19] P. Kramer and Marcos Saraceno, Geometry of
the time-dependent variational principle inquantum mechanics (Springer, 1981).

原子衝突学会誌しょうとつ 第19巻第3号 (2022)[20] A.D. McLachlan, Mol. Phys. 8, 39 (1964).
[21] J. Frenkel, Wave mechanics, advanced general
theory (Clarendon Press Oxford, 1934).
[22] X. Yuan et al, Quantum 3, 191 (2019).
[23] T. Sato et al, Phys. Rev. A 94, 023405 (2016).
[24] J. Arponen, Ann. Phys. 151, 311 (1983).
[25] C. Huber and T. Klamroth, J. Chem. Phys. 134,
054113 (2011).
[26] S. Kvaal, J. Chem. Phys. 136, 194109 (2012).
[27] H. Pathak, T. Sato, and K. L. Ishikawa, J.
Chem. Phys. 153, 034110 (2020).
[28] R. Sawada, T. Sato, and K. L. Ishikawa, Phys.
Rev. A 93, 023434 (2016).
[29] Y. Orimo et al, Phys. Rev. A 97, 023423 (2018).
[30] I. S. Wahyutama, T. Sato, and K. L. Ishikawa,Phys. Rev. A 99, 063420 (2019).
[31] Y. Orimo, T. Sato, and K. L. Ishikawa, Phys.
Rev. A 100, 013419 (2019).
[32] R. Wehlitz et al, Phys. Rev. A 68, 052708(2003).
[33] T. Sato, T. Teramura and K. L. Ishikawa, Appl.
Sci. 8, 433 (2018).
[34] T. Teramura, T. Sato, and K. L. Ishikawa, Phys.
Rev. A 100, 043402 (2019).
[35] H. Pathak, T. Sato, and K. L. Ishikawa, J.
Chem. Phys. 152, 124115 (2020).
[36] H. Pathak, T. Sato, and K. L. Ishikawa, J.
Chem. Phys. 154, 234104 (2021).
[37] F. Lackner et al, Phys. Rev. A 91, 023412(2015).
[38] F. Lackner et al, Phys. Rev. A 95, 033414(2017).
[39] R. Anzaki et al, Phys. Rev. A 98, 063410(2018).
[40] R. Santra, R. W. Dunford, and L. Young, Phys.
Rev. A 74, 043403 (2006).
[41] R. Anzaki, T. Sato, and K. L. Ishikawa, Phys.
Chem. Chem. Phys. 19, 22008 (2017).
[42] Y. Li, T. Sato, and K. L. Ishikawa, Phys. Rev.
A 104, 043104 (2021).

48

2022 年度 役員・委員会
会長
田沼肇（東京都立大学）
幹事
木野康志（東北大学）〔副会長〕

金安達夫（九州シンクロトロン光研）

酒井康弘（東邦大学）

間嶋拓也（京都大学）

的場史朗（高エネルギー加速器研究機構）
運営委員
東俊行（理化学研究所）

岡田邦宏（上智大学）

小田切丈（上智大学）

加藤太治（核融合科学研究所）

金安達夫（九州シンクロトロン光研） 木野康志（東北大学）
高口博志（広島大学）

酒井康弘（東邦大学）

高峰愛子（理化学研究所）

土田秀次（京都大学）

中村信行（電気通信大学）

古川武（東邦大学）

間嶋拓也（京都大学）

的場史朗（高エネルギー加速器研究機構）

山崎優一（東京工業大学）

山下琢磨（東北大学）

常置委員会
編集委員会

委員長：間嶋拓也（京都大学）

行事委員会

委員長：金安達夫（九州シンクロトロン光研）

広報渉外委員会
顕彰委員会

委員長：酒井康弘（東邦大学）
委員長：木野康志（東北大学）

庶務委員会

委員長：的場史朗（高エネルギー加速器研究機構）

編集委員

加藤太治，歸家令果，久間晋，熊谷嘉晃，永田祐吾，
中野祐司，間嶋拓也，山下琢磨，渡邉昇

原子衝突学会誌しょうとつ
第 19 巻 第 3 号（通巻 106 号）
Journal of Atomic Collision Researchⓒ原子衝突学会 2022
http://www.atomiccollision.jp/
発行: 2022 年 5 月 16 日
配信: 原子衝突学会事務局 <acr-post@bunken.co.jp>