新学術領域研究「コンピューティクスによる物質デザイン：複合相関と非平衡ダイナミクス」研究会
東京大学、H23.10.7-8

パルス光の伝播を記述する
Maxwell-TDDFTマルチスケール・シミュレータの開発
筑波大学数理物質科学研究科
計算科学研究センター
矢花一浩

TDDFTを用いた光と物質の相互作用の第一原理計算
光電磁場を記述する
マクスウェル方程式

電子ダイナミクスを記述する
時間依存Kohn-Sham方程式

ψ i (r , t )E (r , t ), B (r , t )基本となる方程式は、
時間依存Kohn-Sham方程式

 e  2
∂δE1 i ψ i = − i∇ + A  ψ i − eφψ i + xc ψ i∂tδn2m c n = ∑ψi2i


1 ∂ A  2  4π −∇ A=j22c ∂tc2

マクスウェル方程式

 e 1ψ i*  − i∇ + A ψ i − c.c.
j=∑2im ic 2∇ φ = −4π {enion − ene }通常の（弱い）光と物質の相互作用では、マクスウェル方程式とシュレディンガー方程式は分離する。
- 電磁場を摂動と扱い、シュレディンガー方程式から線形応答係数（誘電率、導電率）が計算される。
- 光の波長（μm）と、電子ダイナミクス（nm以下）の空間スケールの乖離から、巨視的マクスウェル方程式を解く。
（仮定 ： 電子の平均自由行程 << 巨視的空間スケール～光の波長）

今日の先端の光科学では、様々な局面で上記のアプローチが破綻する。
- ナノサイズの物質の光応答では、空間スケールが分離しなくなる。
t
  

D (r , t ) = ∫ dt ' dr 'ε (r , r ' , t − t ' )E (r ' , t ' )- 超短パルス光では、振動数ではなく時間を変数とした記述が必要となる。
tD (t ) = ∫ dt 'ε (t − t ' )E (t ' )- 高強度パルス光では、光電場による電子の非線形ダイナミクスが顕著になり線形応答理論が使えなくなる。（マクスウェル方程式とシュレディンガー方程式は分離しなくなる。）

D = D[E ] ≠ εE

高強度超短パルス光と物質の相互作用を第一原理計算により記述する試み

１．与えられたパルス光がもたらす電子・フォノンダイナミクス- 結晶中の電子ダイナミクスを記述する方法論- コヒーレントフォノン生成の第一原理計算- 電子・ホールプラズマ生成と光絶縁破壊

２．高強度超短パルス光伝播の記述- Maxwell-TDDFTマルチスケール・シミュレータ- 予備的な計算

高強度超短パルス光と物質の相互作用を第一原理計算により記述する試み

１．与えられたパルス光がもたらす電子・フォノンダイナミクス- 結晶中の電子ダイナミクスを記述する方法論- コヒーレントフォノン生成の第一原理計算- 電子・ホールプラズマ生成と光絶縁破壊

２．高強度超短パルス光伝播の記述- Maxwell-TDDFTマルチスケール・シミュレータ- 予備的な計算

無限周期系（固体）の実時間・実空間計算空間に関して一様で、時間に関して周期的な電場の扱い

周期ポテンシャルeE(t)z
周期的ポテンシャル：ブロッホの定理

()

 

ik R
ψ nk r + R = e ψ nk (r ),() h r + R = h( r )一様な外場ポテンシャル eE(t)z があると使えない。

ゲージ変換をして、ベクトルポテンシャルを使うと、周期性が回復する。


∂A
E = −∇φ −∂t
φ = eEz ⇔ A = eEzˆt 1   e  2
 ∂i ψ (t ) =  p − A(t ) − eφ (r , t )ψ (t )c∂t 2m ベクトルポテンシャルの方程式： 分極場のダイナミクス
一様電場はレーザー電場＋分極電場A(t ) = Alaser (t ) + Apolarization (t )

誘電体
 

E = Elaser + E polarization

Elaser (t )～分極の運動方程式

d Apolarization (t )
2dt2

 
4π
=
dr j (r , t )∫V cell∂A(t )E (t ) = −∂tt A(t ) = −c ∫ E (t )dt j (r , t )A(t )

時間依存Kohn-Sham方程式
 e  2∂δE xc
1 −∇+i ψ i =iAψ−eφψ+ψi ii∂tc δn2m n = ∑ψii21 *  e  j=ψ i  p + A ψ i − c.c. ∑2m i  c 

Bertsch, Iwata, Rubio, Yabana, Phys. Rev. B62(2000)7998.

Si結晶にパルス光を入射した際に起こる電子ダイナミクス

Electric Field (a.u.)

I= 3.5×1014 W/cm2, T=50 fs, ħω=0.5 eV0102030Time (fs)Density change from the ground state （１１０）40Ground state density （１１０）Etotal vs Elaser : スクリーニング効果
誘電体 
E = Elaser + E polarization

Calculations: I=1012W/cm2, ħω=1.03eV
0.006

Bulk Si

0.004

electric field

0.002
00.002

Laser electric field (D)0.004

Total electric field (E)= Laser + Polarization

0.006

0510

D = ε (ω )E

15time [fs]20

ε ≈ 142530(TDDFT )Elaser (t )～TDDFT（断熱局所近似）のもとでの誘電関数（Si)

定量的には不十分- バンドギャップの過少評価- 電子・ホール相互作用

高強度超短パルス光と物質の相互作用を第一原理計算により記述する試み

１．与えられたパルス光がもたらす電子・フォノンダイナミクス- 結晶中の電子ダイナミクスを記述する方法論- コヒーレントフォノン生成の第一原理計算- 電子・ホールプラズマ生成と光絶縁破壊

２．高強度超短パルス光伝播の記述- Maxwell-TDDFTマルチスケール・シミュレータ- 予備的な計算

超短パルス光によるコヒーレントフォノン生成過程の記述
Y. Shinohara, K. Yabana, Y. Kawashita, J.-I. Iwata, T. Otobe, G.F. Bertsch, Phys. Rev. B82, 155110 (2010).

Raman active direction
Raman inactive direction

Si:ポンプ=プローブ分光
M. Hase, et.al,
Nature 426, 51 (2003)

E//[110]⑥- ポンプ光が電子励起を起こし、コヒーレントフォノンを生成。.
- プローブ光の反射率変化を、ポンプ光との時間差を関数にして測定。

⑤⑧⑦② ④① ③E//
[001]

T

Polarization [110]: Raman active

4

I= 1×1012 W/cm2,
T=16 fs (8 mJ/cm2),
ħω=2.06 eV

20Applied laser pulse Eext(t)2Force on ① [eV/pm]

1020304050E//[110]F1,x×105∫ F (t )dt ≅ 0⑧⑤⑦② ④① ③⑥1x1020304050

F1,y×105∫ F (t )dt ≅ 01y10203040

F1,z×4×106T∫ F (t )dt ≠ 001020

time [fs]301z40dispacement [pm]

40

2015
1050-5101520 020151050-5101520 020
151050-51015
20 0z

50

6420-2-4-650

R1,z×10401020304050

原子に働く力[eV/nm]

単位胞あたりの励起電子数

Laser pulse
Laserpulseレーザー光の振動数が直接バンドギャップ(2.4 eV in DFT)を越えると、コヒーレントフォノンの生成機構が変化する

ISRS (below gap)
DECP (above gap)

Atomic coordinates are fixedduring the calculation.

TDDFTは、コヒーレントフォノンの２つの生成機構を統一的に記述できる。
Y. Shinohara, K. Yabana, Y. Kawashita, J.-I. Iwata,
T. Otobe, G.F. Bertsch, Phys. Rev. B82, 155110 (2010).

ラマン活性な[110] 方向に光電場を照射した際の電子ダイナミクス
E//[110]
I= 1×1012 W/cm2,
T=16 fsħω=2.5 eV⑧⑤⑦② ④① ③⑥Color map changed
by factor 40

Calculated direct bandgap

コヒーレントフォノンの振幅・位相・平衡位置のずれq(t ) = −q0 cos(ω phonont + φ ) + q

振幅 q0ISRS

DECP

位相 φ

平衡位置のずれ q

Measurements:
M. Hase et.al, Nature 426, 51 (2003)
D.M. Riffe et.al, Phys. Rev. B76, 085207 (2007)
K. Kato et.al, JJAP 48, 100205 (2009)

Laser frequency

高強度超短パルス光と物質の相互作用を第一原理計算により記述する試み

１．与えられたパルス光がもたらす電子・フォノンダイナミクス- 結晶中の電子ダイナミクスを記述する方法論- コヒーレントフォノン生成の第一原理計算- 電子・ホールプラズマ生成と光絶縁破壊

２．高強度超短パルス光伝播の記述- Maxwell-TDDFTマルチスケール・シミュレータ- 予備的な計算

高強度パルス光を用いたレーザー加工
光から物質へのエネルギー移行
Femtosecond laser machining

Short pulse laser enables non-thermal,
precise operation on materials

1.4 nsM.D. Perry et.al,
J. Appl. Phys. 85, 6803 (1999)

Short pulse laser
dig a hole on enamel

350 fsさらに強度を増すと⇨光絶縁破壊
Electric field (au)T. Otobe, M. Yamagiwa, J.-I. Iwata, K.Y. T. Nakatsukasa,
G.F. Bertsch, Phys. Rev. B77, 165104 (2008)

0.2
ext5d14_400
tot5d14_400

0.1
0.0
-0.1

5 x 1014 W/cm2-0.2
0200

400

600

800

1000

time (au)Elaser(t) vs Etotal(t)ext1d15_400
tot1d15_400

0.1
0.0-0.1

Thresholdfor breakdown -0.2 0

1 x 1015 W/cm2
200

400

600

800

1000

time (au)0.6Electric field (au)ダイアモンド振動数: 3.1eV(直接ギャップ以下)
パルス長: 16fs

Electric field (au)0.2

ext2d15_400
tot2d15_400

0.4
0.2
0.0-0.2

2 x 1015 W/cm2-0.4
-0.6
0200

400

600

time (au)

800

1000

TDDFTでの光絶縁破壊のメカニズム(1x1015 W/cm2, 3.1eV, 40fs)
最初は誘電体として応答
リレーションシップ ID rId6 のイメージ パーツがファイルにありませんでした。

次第に電子励起が進む, > 15 fs- 外場 Eext(t) と内部電場 Etot(t) に位相差- 励起電子数の急激な増加⇒ 絶縁破壊
金属的な応答, > 25 fs- 励起電子数、エネルギー移行が止まる

Eext(t)Etot(t)励起電子数

励起電子数 0.4/atom でのプラズマ振動数

 4π nex  ≈ 4eV mε (0 ) 
これは照射したレーザーの振動数に近い。

ω p = レーザーから電子へのエネルギー移行

高強度超短パルス光と物質の相互作用を第一原理計算により記述する試み

１．与えられたパルス光がもたらす電子・フォノンダイナミクス- 結晶中の電子ダイナミクスを記述する方法論- コヒーレントフォノン生成の第一原理計算- 電子・ホールプラズマ生成と光絶縁破壊

２．高強度超短パルス光伝播の記述- Maxwell-TDDFTマルチスケール・シミュレータ- 予備的な計算

高強度超短パルス光の物質内伝播の第一原理計算
物質中の光伝播（マクスウェル方程式）

光電場による電子ダイナミクス（時間依存Kohn-Sham方程式）


E (r , t ), B (r , t )ψ i (r , t )波長 : µm

単位胞: nm

２つの異なる空間スケール

弱い電磁場： 摂動論により誘電関数を用いた分離が可能。
強い電磁場： マクスウェル方程式とシュレディンガー方程式が分離しなくなる。

D = D[E ] ≠ ε (ω )E⇒ マクスウェル=TDDFTマルチスケール・シミュレータ

光電磁場の１次元的な伝播と反射ｘ

マクスウェルの方程式ｚε (z ) ∂ 2∂2A(z, t ) − 2 A(z, t ) = 022∂zc ∂t A(r , t ) = A(z, t )xˆIndex of refraction

真空誘電体
ε=20n= εVelocity of wave
cnReflectancev=λ = 800nm, ω = 1.55eV

1− n R=+n12パルス光がSi表面に垂直入射する場合のマクスウェル-TDDFTマルチスケール・シミュレーション
ｘVacuumSiZ

巨視的格子点（μm）上で離散化されたベクトルポテンシャルの時間発展A = xˆA(Z , t )1 ∂2∂24π()()
AZ,t−AZ,t=J (Z , t )∂Z 2c 2 ∂t 2c
J (Z , t )巨視的ベクトルポテンシャルと巨視的電流を介して結合

各巨視的格子点上で独立に、電子ダイナミクス（nmスケール）を記述する時間依存Kohn-Sham方程式
A(Z , t )J ( Z , t ) = ∫ dr je,ZΩ() *e*je ,Z =
∇−
∇−nA
ψψ
ψ
ψ
∑ i ,Z i ,Z i ,Z i ,Z 4πc e,Z2mi i e  2∂δE1 i ψ i , Z = − i∇ + A  ψ i ,Z − eφZψ i ,Z + xc ψ i ,Z∂tδn2m c 2∇ φZ = −4π {enion − ene,Z }弱い光パルス：線形応答領域（電磁場と電子ダイナミクスが分離する場合）
I=1010W/cm2
Vacuum

光の振動数：1.55eVSi

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

Z=0 [µm]

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

やや強い光パルス：非線形励起を伴う場合（電磁場と電子ダイナミクスは分離できない）
I=5 x 1012W/cm2
Vacuum

Si[µm]

Z=0 [µm]

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

計算サイズｘｚ
3.5µm, 256 points∆z = 13nm
Microscopic (Schroedinger)
spatial grid: 163
k-points : 83 (reduced by symmetry)

Macroscopic (Maxwell)
spatial grid: 256Time step (commom) = 16,000
Parallelize(256 points)*(4 k-points) = 1024 cores (MPI only)
11 hours at
SGI Altix ICE 3800EX (ISSP, Univ. Tokyo)
k点に関して十分な収束を得るたえにはコア数がまだ不足。

I=1d10 W/cm2

光電場

単位体積当たりの電子励起エネルギー

I=1d12 W/cm2

Si表面に誘起される電子・ホールプラズマ（Preliminary)
Si表面での光パルスの反射率

2-photonabsorption(n − 1) + k(n + 1)2 + k 22

R=

Si表面での電子励起

2e-h plasma freq.
= laser freq.

誘電応答
レーザー
強度

２光子吸収による電子励起
多電子励起による電子・ホールプラズマ生成関連する測定“Generation of dense electron-hole plasma in silicon”
K. Sokoowski-Tinten, D. von der Linde, Phys. Rev. B61, 2643 (2000)

λ=625nm, 100fs pulse

Si

Pump pulse, s-polarizatione-h density whereplasma freq. = laser freq.

Probe pulse, p-polarization

Reflection ofprobe pulse
Exp.
Exp.

Electron-hole density

まとめ
近年の光科学では、誘電率を用いた巨視的マクスウェル方程式では記述できない光-物質相互作用が主要な課題となっている。
- 空間的な非局所性- 時間的な非局所性- 電子ダイナミクスの非線形性与えられたパルス光が生じる電子・フォノンダイナミクス（時間依存密度汎関数理論による固体中の電子ダイナミクス計算）- コヒーレントフォノン生成- 電子・ホールプラズマ、光絶縁破壊
Maxwell-TDDFTマルチスケール・シミュレーション- 数千コアを要する大規模計算により可能- 長波長近似は仮定、摂動論は用いない- パルス光の１次元伝播の記述にほぼ成功
今後の展開- パルス光と物質の様々な相互作用への応用（Pump-probe, CEP, アト秒科学．．．）- 理論の精緻化（電子衝突、電子ホール相互作用； 実電荷・実電流の扱い、多次元化）- 高速化（数万コアの利用、GPUの利用）- 大規模化（多次元の光伝播、ナノ構造と光の相互作用）