文部科学省科学研究費新学術領域「コンピューティクスによる物質デザイン：複合相関と非平衡ダイナミクス」
平成25年度第１回研究会
東京大学本郷キャンパス武田先端知ビル
平成２５年７月８－９日

極限的パルス光と結晶の相互作用を記述する
第一原理計算筑波大学数理物質系物理学域
計算科学研究センター
矢花一浩
共同研究者：
佐藤駿丞（筑波大数理）
谷口億宇（筑波大計科セ）Contents
公募研究課題名：
光と電子のダイナミクスを記述する
第一原理マルチスケールシミュレーション法の開発

１．高強度・超短パルスと結晶の相互作用と、非線形電子ダイナミクス（フェムト秒・アト秒科学）
２．Maxwell+TDDFT 第一原理マルチスケール計算
３．シリカ薄膜を透過するパルス光：計算と実験
４．汎関数の改善：メタGGA、ハイブリッド汎関数（現状報告）

物質（結晶）と高強度・超短パルス光の相互作用
強い光

短い光

光電場の強度と、物質中で電子が感じる電場が同程度となる。

1eV（赤外から可視領域）の光の周期と同程度のパルス光

非線形電子ダイナミクス

フェムト・アト秒ダイナミクス

フェムトテクノロジー:レーザーパルスを用いた非熱的加工

アト秒X線パルスを用いたレーザーパルス波形の実時間記録

R.R. Gattass, E. Mazur, Nature Photonics 2, 220 (2008).

Joint LMU-MPQ Laboratory of Attosecond

アト秒科学のフロンティア： 原子分子から固体へ

F. Kransz’s group (MPQ),
Nature 493, 70(2013, Jan.3)

パルス光により誘電体（SiO2）中に自由キャリアが生じ、電流が流れる。

1013 W/cm 2

物質中のパルス光伝播をミクロに記述する
巨視的電磁場の運動を記述する
マクスウェルの方程式
  1 ∂B=0
∇× E +c ∂t∇B = 0∇D = 4πρ  1   ∂D ∇ × H =  4πj +c∂t 

物質の性質は構成方程式により取り入れられる。
物質を特徴づけるのは、（非）線形分極率   D(r , t ) = E (r , t ) + 4πP(r , t )ttt(1)( 2)(t − t ' , t − t ' ')Eβ (r, t ')Eγ (r, t ') + (t − t ')Eβ (r, t ') + ∫ dt '∫ dt ' 'χαβγ
Pα (r , t ) = ∫ dt ' χαβ

極限的な高強度パルス光では、もはや摂動展開は使えない。
与えられた電場により生じる分極（電場の局所汎関数）をどうやって求めるか？
  
P(r , t ) = P E (r , t )[]結晶中の電子ダイナミクスを記述する時間依存Kohn-Sham方程式（長波長極限）時間に依存するブロッホ関数
 1    e  2 ∂i unk (r , t ) =  p + k − A(t ) + V (r , t )unk (r , t )c∂t
 2m  unk (r , t ) = unk (r + a , t )
1 ∂A(t )⇔E (t ) = −c ∂t与えられた電場 E(t) に対して、分極 P(t) が得られる。

tA(t ) = −c ∫ dt 'E (t ')zˆ
 e 
 1*.
J (t ) = ∫ d r −i∇+A−ccψψ.
∑i i i2imctP(t ) = ∫ dt ' J (t ').

与えられた電場に対して、結晶中の時間依存Kohn-Sham方程式を解くことにより、摂動展開に依らない構成方程式を得ることができる。

Maxwell + TDDFT第一原理マルチスケール計算
パルス光が物質表面に垂直入射巨視的格子点（μm）上で離散化されたベクトルポテンシャルの時間発展
1 ∂24π
∂2()
()
AZ,tA
Z,tJ (Z , t )−=c 2 ∂t 2c∂Z 2J (Z , t )巨視的ベクトルポテンシャルと巨視的電流を介して結合各巨視的格子点上で独立に、電子ダイナミクス（nmスケール）を記述する時間依存Kohn-Sham方程式
A(Z , t )2 e 
δE xc1
∂

()
ψ,
ψ−
φψ+
ψ i ,Z
=−∇+iiAZtei ,Zi ,ZZ i ,ZJ ( Z , t ) = ∫ dr je,Z2δ∂
cntmΩ2eje ,Z =∑ (ψ i*,Z ∇ψ i ,Z −ψ i ,Z ∇ψ i*,Z ) − ne,Z A ∇ φZ = −4π {enion − ene,Z }

2mi i

4πc

バルクSiに、強いパルス光が垂直入射する場合

I=5 x 1012W/cm2Vacuum

光の振動数：1.55eVSi

直接バンドギャップ：
2.4eV(LDA計算）[µm]

Z=0 [µm]

Z=0.8 [µm]Z=1.6 [µm]

数サイクルパルスの薄膜透過計算： 10µm SiO2 （αクォーツ）
Macro-grid 1000 points
Red: 5x1013 W/cm2
Blue: 1x1011 W/cm2

9

10 μm SiO2 薄膜透過波形

Blue: 1x1011 W/cm2

Red: 1x1013 W/cm2

Red: 2x1013 W/cm2

Red: 5x1013 W/cm2
10

計算では、5x1013W/cm2の強度で著しい吸収が起こる。
実験では、1x1014W/cm2を超える強度まで変化していないという報告がある。

より定量性のある記述を目指して：誘電関数のTDDFT計算（バルクSi)
メタGGAを用いた誘電関数

ハイブリッド汎関数を用いた誘電関数

V.U. Nazarov, G. Vignale, PRL107, 216402 (2011)[  ]ε xc ρ , ∇ρ ,τ

J. Paier, M. Marsman, G. Kresse, PRB78,121201(2008)

HSE：非局所交換ポテンシャルを25%混合 1 2τ (r ) = ∑ ∇ψ i (r )2 iVS98: T.V. Voorhis, G.E. Scuseria, JCP109,400 (1998)
exp
TDDFT
(meta-GGA)

345

6HSE03, J. Heyd, G. E. Scuseria, and M. Ernzerhof,
J. Chem. Phys. 118, 8207 (2003)

メタGGA汎関数（TBmBJ）
F. Tran, P. Blaha, PRL102, 226401 (2009)を用いたバルクSiの誘電関数の実時間計算

LDA

TB-mBJ

計算のコストは、LDAのおよそ2－3倍

Exp

TB-mBJでは、ギャップは改善するが、エネルギーをシフトしているだけのよう。

（VS98では、誘電関数の振動数依存性も変わっていた）

TB-mBJ

VS98: T.V. Voorhis, G.E. Scuseria, JCP109,400 (1998)
exp
TDDFT
(meta-GGA)

LDAExpSiO2の光励起： 汎関数による差単一セルでの電子ダイナミクス計算
パルス光照射による電子励起エネルギー
LDA vs meta-GGA

LDA
4光子吸収

励起電子数
（#/atom)meta-GGA
5光子吸収

パルス光の強度汎関数の差（ギャップの差）は、パルス光強度を2－3倍変化に相当。

ハイブリッド汎関数を用いた誘電関数（バルクSi)
長距離部分をスクリーニングした、ハイブリッド汎関数
HSE03, J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 118, 8207 (2003)

我々の計算（実時間・実空間・擬ポテンシャル）
VASPによる計算（平面波・PAW)
J. Paier, M. Marsman, G. Kresse, PRB78,121201(2008)

ハイブリッド汎関数計算： GPUを用いた加速非局所交換ポテンシャルの計算： 高速フーリェ変換を用いる2 e2
πe4iKrK*u  (r )un 'k ' (r ')wnk (r ') = − ∑u  (r ) e f n 'k ',nk    2∫ dr ' r − r' ∑ n 'k '  n 'k 'n 'k 'n 'k ' K
K +k −k'

f nK'k ',nk = −iKr *    1dr
e un 'k ' (r )wnk (r )∫ΩΩ

HA-PACS（筑波大）：
単位胞： cubic, 8 atoms/cell
k点：12^3 ダウンサンプリングk点：4^3 空間格子：16^3108 GPU(K2090), cuFFT ⇒ 15時間程度
CPU計算（FFTW)との比較：
１ノード（4GPU、または16CPUコア）で、GPUの方が4倍以上速い。

まとめ
光電磁波と電子のダイナミクスを同時に記述するマルチスケール・シミュレーション- 高強度パルス光の固体伝播、光から電子へのエネルギー移行を第一原理から記述することが可能に- 数千コアを要する超大規模計算。
- 実験との比較のために、バンドギャップが正しく再現できる汎関数を用いることが極めて大切。
- メタGGA汎関数（運動エネルギー密度を含む）は、LDAとそれほど変わらないコストで計算可能。定量性はどの程度か？
- ハイブリッド汎関数は、より定量性に優れていそうだが、計算コストが著しく高い。