時間に依存する種々の外場中に置かれたナノ構造の電気伝導及びスピン伝導現象のシミュレーション

神戸大学大学院 工学研究科 電気電子工学専攻
相馬聡文 （渡邉班）
共同研究者：貝野昭造，坂本理朗，小川真人

Outline
１．フェルミ行列展開法を用いた強束縛分子動力学法の高速化
２．波束ダイナミクスを用いた電気伝導度計算と，光照射したグラフェンの電気伝導への応用
３．量子ポンピング現象のシミュレーション手法の開発と，グラフェンナノリボンにおける量子ポンピング効果への適用．

フェルミ行列展開法を用いた強束縛分子動力学法の高速化

研究目的
9グラフェンは種々の優れた物性を持つが，フェルミ面においてバンドギャップが存在しないためスイッチング特性に問題がある

9スイッチング素子への応用を考えるためにはフェルミ面でのバンドギャップを発生させる必要がある
example・・・
・ グラフェンナノリボン(GNR)・ バイレイヤーグラフェン
GNRに対して歪みを印加する事により，更なる電気伝導性の制御の可能性を探る．特に，フェルミ面におけるバンドギャップを制御する可能性について考察。
引っ張り歪みを加えたＧＮＲの安定構造→ 外力の下での構造緩和計算が必要．
→強束縛分子動力学法（Tight-Binding Molecular Dynamics (TBMD) ）4C. H. Xu, et al, J. Phys. Condens. Matter . 4, 28 (1992).

解析モデルfree-standingの状態で構造緩和済みの7×6 AGNR

lx歪み率(Stretching ratio: SR)

unit cellzyz軸(輸送)方向への引っ張り歪み

l '−lSR =× 100[%]
ll’

引っ張り歪みを加えると電子状態・電気伝導特性はどう変化するか？
5

TBハミルトニアン行列Hを求める①① Slater-Kosterテーブルと構造緩和の計算フロー
ハミルトニアンのパラメータ
全エネルギーEtot(=Ebs+Erep)を求める
各原子が受ける力Ftot(tn)を求める
位置ri(tn+1)を求める
力Fi(tn+1)、速度vi(tn+1)、を求める②(スケーリング)関数を用いる③② 結合エネルギーEbs
反発エネルギーErepを求める④③ Ftot=-∇Etot の関係からFtot(tn)を求める⑤④，⑤ 速度Verlet法により求める

系を安定位置に収束させるために、速度に比例する摩擦係数μを用いた計算①～⑤の操作を収束するまで繰り返す
速度Verlet法 次の時刻での位置と速度を現在の時刻の情報(位置，速度，力)から計算する方法

6

TBMD法において、結合エネルギーE はTBハミルトニアン
従来の対角化法
行列Hの固有値の総和である
bs

温度による影響を取り入れると従来の対角化法では固有値εiはハミルトニアン行列Hを対角化して求めていた問題点9εiの計算量が原子数Nの3乗に比例するため、MD計算において，大きな系を扱う場合に計算時間が膨大となる。
9並列計算に不向き計算時間を減らすためオーダーN法である
Fermi operator expansion(FOE）を導入

7

FOEにおける結合エネルギーE の計算はFOEの導入
bs

ハミルトニアン行列H ×フェルミ行列Fの演算で求められる。

フェルミ行列Fはフェルミ分布関数fにハミルトニアン行列Hが入った形で与えられる。
フェルミ分布関数関数の中に行列が入っているためこのままではフェルミ行列Fを求めることが出来ない・・・

FOEではフェルミ分布関数をChebyshev展開することで行列多項式p(H)に持ち込む

8

Chebyshev近似

近似の精度を決定

Chebyshev展開において多項式p(H)は以下の漸化式を満たす

上式でハミルトニアン行列H×Tｊ(H)の計算量はn
⎡ H11
⎢H
⎢ 21
T j +1 ⇒ ⎢ H 31⎢
⎢ M⎢⎣ H n1

H12HnTj

H13 L H1n ⎤ ⎡T11j T12j T13j L T 1nj ⎤ ⎫⎥ ⎪⎥ ⎢ j⎥ ⎪⎪⎥ ⎢T21j⎢
⎥ ⎬n
⎥ T31
⎥ ⎪⎥ ⎢⎥ ⎪⎥ ⎢M⎥ ⎪⎭⎥⎦ ⎢T j
⎣ n1⎦

T j +1⎡T11j +1 T12j +1 T13j +1 L T j1+n1 ⎤ ⎫⎢ j +1⎥ ⎪⎢T21⎥ ⎪⎪⎥ ⎬n= ⎢T31j +1⎢⎥ ⎪⎢ M⎥ ⎪⎢T j +1⎥ ⎭⎪
⎣ n1
⎦

n×n

Tｊ+1(H)の1列あたりn×n、それをn個の列について求めるので全計算量はn3となる。

ハミルトニアンのスパース性を利用 + cutoff半径を導入することで全計算量をnに減らすことができる

9例として…
対角化法とFOEの比較
1 unit cell あたりの原子数が14個のAGNRに

CPU times (s)

4000
3500

diagonalization

3000

FOE serial

ついて、unit cell の数を増やしていったときのTBMD(free-standing)にかかる時間を計算した。
なお、並列化にはopenMPを用いている。

2500
2000
1500
1000

H-F行列積における並列化

500
0020406080100120
number of atoms

140160180

F1700

F3F4FOE serial

600

FOE 2 threads

500
CPU times (s)F21/10

FOE 4 threads

400

FOE 8 threads

300

FOE 12 threads

CPU 1 CPU 2 CPU 3 CPU 4200

H x F1

1000020406080100120
number of atoms

140160180

H x F2

H x F3

H x F4

大幅な時間短縮の実現 10

波束ダイナミクスを用いた電気伝導度計算と，光照射したグラフェンの電気伝導への応用

グラフェンへの光照射
本研究

光を照射してエネルギーギャップを発生させる手法に着目
グラフェンに対して面直に円偏光を照射したときの時間依存シュレディンガー方程式⎡ ⎛eEeE⎞⎤⎛⎞H (t ) ≡ hvF ⎢σ x ⎜ k x − 0 cos ωt ⎟ + σ y ⎜ k y + 0 sin ωt ⎟⎥hωhω⎠⎦⎝⎠⎣ ⎝σ x ( y ) :パウリのスピン行列 E 0 : 照射光の電界強度∂ih φ (t ) = H (t )φ (t )∂t

適当なkx,ky,ω,E0を与えて波動関数の時間発展を数値的に求めると、与えたkx,ky毎に波動関数の時系列が得られる。

φk ,k , A (t )xyReIm
フーリエ変換t

この波の中にどのような角振動数（=エネルギー）の波がどの程度の強度で含まれているかがわかるグラフェンの動的バンド構造
E0 = 0, ω = 2π [rad/fs]

E0 = 7.0 [V/nm ], ω = 2π [rad/fs]

4422

10-1-210-1-2-3-3-4-6.0-4.0-2.00k x [rad / nm ]2.04.0-46.0-6.0

E0 = 10.0 [V/nm ], ω = 2π [rad/fs]

た。

30-1-2-3-4-6.0-4.0-2.00k x [rad / nm ]

2.0

4.0

6.0

Energy [eV]1-4.0-2.00

2.0

k x [rad / nm ]4.06.0k y = 0 としてスペクトル強度を明るさでプロットし

42Energy [eV]3Energy [eV]

3kx = 0電磁波を照射しない場合、でエネルギーバンドが交差していることがわかる(元のバンド構造を再現)。
円偏光の光を照射した場合、元の交差点においてバンドギャップ(動的バンドギャップ)が生じていることがわかる。
また、照射光の強度を大きくするとより大きなバンドギャ プが発生していることがわかる

研究目的
半導体デバイスへの利用に向けた、光によるグラフェンの電気伝導特性制御のシミュレーション
チャネル部分に光を照射することで電子の透過率がどのように変化するかを、ガウス波束の時間発展計算によって求める
時間発展計算ガウス波束
ψ(x)t=t0 t=t1μσx・・・・・・ x-2様々な波長の重ね合わせを表現できる

透過率計算モデルx-1

x0x1x2・・・・・・

x

時間依存である電磁波を考慮した計算をするためd

入射波

幅d,高さVのチャネル部分に侵入した入射波は、

V

十分時間経過後は透過波と反射波とに分かれる。
十分時間経過反射波

透過波x

透過波、反射波を積分することで透過率・反射率を算出することができる。

透過率計算シミュレータの妥当性構築した時間発展による透過率計算シミュレータの妥当性を定常状態において確認する
kｘkθ

VEkｙθ

yz
グラフェン

dxx,yk電子が入射角θで幅d,高さVのバリアに衝突するときの透過率の角度依存性を時間発展計算で求め、解析解と比較した。
M. I. Katsnelson et al Nature phisics vol.2 (2006)

解析解
時間発展解

1

d=100[nm],V=0.2[eV],E=0.08[eV],k=0.125[rad/nm]

0.8Transmittanceとして比較したものが左図である。
0.6

解析解のような大きな振幅は得られていないものの、振動の周期は解析解と時間発展解でほぼ同じになっていることがわかる。

0.4

解析解と大きな差異は見られないことから、構築した時間発展による透過率計算シミュレータの妥当性を確認できた。

0.20-80-40040

Incident angle [degree]80グラフェンの電気伝導特性のモデル
hω, E0d

チャネル：長さd,高さVy

照射光：強度E0,振動数ωチャネルxEFEgEΔEΔkk0V (= EF )E=V=EF(=1.0 [eV])となるように波束の波数k0をk0=1.53 [rad/nm]として計算した。

2⎛ ev E ⎞2E g = 4⎜ F 0 ⎟ + (hω ) − hω⎝ ω ⎠T. Oka and H.Aoki, Phys. Rev. B 79, 081406(R) (2010)

E0=10.0 [V/nm], ω=2π=6.28 [rad/fs]とすると、動的バンドギャップEg≒1.08 [eV]
よって、E0=10.0 [V/nm]のとき、V≒EF±0.5[eV]=1.5, 0.5 [eV]で透過率が増加すると考えられる

時間発展計算結果
Existence probability

Position [nm]Time [fs]E0=10.0 [V/nm]

Existence probability

Position [nm]Position [nm]E0=5.0 [V/nm]
Time [fs]Time [fs]E0=0.0 [V/nm]

Existence probability

d=10 [nm],V=1.0 [eV],k0=1.53 [rad/nm], ω=6.28 [rad/fs]のときの時間発展計算結果の一例を示す

横軸にx方向位置、縦軸に時間をとり電子の存在確率を明るさでプロットした。
初期状態で50 nmの位置にあるピークが時間発展とともに正方向に進み、100 nmにの位置にあるチャネルに進入している。
照射光の強度を大きくすることで、チャネルを透過する電子の存在確率が小さくなることがわかる。

光の照射によってグラフェンの電気伝導特性を制御できる

透過率の照射光強度依存性
1d=2 [nm]d=5 [nm]d=10[nm]TransmittanceTransmittance

0.80.60.4

0.3

0.2

0.10.20d=2 [nm]d=5 [nm]d=10 [nm]0.4024681000.30.50.70.91.11.31.51.7

E 0 [V/nm]

V [eV]

チャネル長d=2,5,10 [nm]において、チャネル高さV=1.0 [eV]で固定して横軸に照射光強度、縦軸に透過率をとったもの

チャネル長d=2,5,10 [nm]において、照射光強度E0=10.0 [V/nm]で固定して横軸にチャネル高さ、縦軸に透過率をとったもの

チャネル長によらず照射光強度を大きくすることで透過率を小さくできることがわかる。

dの値によらず動的バンドギャップEg≒1.08 [eV]に対応するV=0.5, 1.5 [eV]付近で透過率が大きくなっていることがわかる。

チャネル長が長いほど小さい照射光強度で透過率を0に近づけることができ、鋭敏に透過率が変化することがわかる。

これらのことより、前述の電気伝導特性モデルのようなエネルギーバンドの変化が起こっていると考えられる。

量子ポンピング現象のシミュレーション手法の開発と，グラフェンナノリボンにおける量子ポンピング効果への適用．

Spin Density Functional Tight-Binding Method

Hσμν= H μνDFTBr σ rr+ ∫ dr Φ μ (r )Veff (r )Φν (r )3Pseudo Atomic Orbital*(1DFTB≅ H μν + S μν δVασ + δVβσ
2S μν = Φ μ Φν)

DFTB based LCAO Hamiltonian
rrr− h 2ΔDFTB
H μν = Φ μ+ ∑ VNA,γ (r ) + Vxc nαisolated + nβisolated (r ) + ∑ Vnl,γ (r ) Φν ,2m γ =α , βγ =α , β[r atom r
δV = ∫ drδVeff (r )nα (r )σα]μ ∈ α ,ν ∈ βσrrrσδVeff (r ) = VδH (r ) + VLSDA [ζ ](r )Poisson’s equation is solved in 3D gridwith appropriate boundary conditions.

σr δEVLSDA [ζ ](r ) =σrr[n(r ), ζ (r )]rδζ (r )LSDA
xc

D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 45 (1980) 566.
J. Perdew and A. Zunger, Phys. Rev. B23 (1981) 5048.

r↑ r
↓ r
ζ (r ) = n (r ) − n (r )20Non-equilibrium Green function MethodSpin dependent charge density

rrnσ (r ) = ∑ qασ nαatom (r )Spin dependent density matrixσρ μν=αSpin dependent Mulliken population[σqασ = ∑∑ Re ρ μνS μνμ∈α ν]

1∞[()(S ,σD ,σ()(E ) f E − μ D
μdE
AEfE−+AμνSμν∫π −∞S ( D ),σ(E ) = ∑ Gμδσ (E )ΓδγS ( D ),σ (E )Gνγσ (E )Aμν[δγ()](E )]G (E ) = δΓδγS ( D ),σ = i ΣδγS ( D ),σ − ΣγδS ( D ),σ∑ [ES − H − Σ (E ) − ΣνσμνμνS ,σμνD ,σμν*σνξμξμ ,ν , ξ ∈ Device Region
Spin dependent current

Total energy in the device region[() (e ∞I D = ∫ dET σ (E ) f E − μS − f E − μ Dh −∞σT (E ) = ∑ Γμν (E )Gσ

D ,σμμS − μ D = e VDSσμν)](E )Γνξ (E )G (E )S ,σσ*μξ

Etot = Eband + Erep
Eband = ∑σσσH μνρ μν∑μν∈device21)]

Edge states in Zigzag Edged Graphene Nanoribbon (ZGNR)
Spin-unpolarized case

N=4EFk x a0 3 / π
Partly flat band at the Fermi level(its wavefunction is mainly localized at the “zigzag” edges)
High Density of States at the edges around the Fermi level!

LSDA Results(Spin-polarized case)Spin density profile

3Potential profile

Energy [eV]210

Spin Polarized CaseSpin Unpolarized Case

-1

Large local density of statesat the edges around the Fermi level.

-2-30.00.5

kza/π1.0

Spin-dependence of the electron-electroninteraction
Locally spin-polarized (edge-spin polarized)
case is energetically stable.
(Energy gap appears)

Theoretical Formulation of Adiabatic Pumping
Particle transport at zero bias voltage bychanging some properties of a phase-coherent
mesoscopic conductor periodically in time.
D. Thouless, PRB 27, 6083 (1983)
Pumped charge by small deviation of two parameters

⎡ dnα
⎤dnα
δQα = e ⎢
δX 1 +
δX 2 ⎥dX 2⎣ dX 1⎦

Emissivity→Scattering Matrix
orGreen’s functionX 2

Pumped charge during the cycle:⎡ dnα dX 1 dnα dX 2 ⎤+Qα = e ∫ dt ⎢⎥dtdXdtdX2⎣ 1⎦0⎡ ∂ dnα∂ dnα ⎤−= e ∫∫ dX 1dX 2 ⎢⎥∂∂XdXXdX221⎦⎣ 1STTω TQα
Pumping current: I α =2π1

ω>> τ dwellSX 1 (t ) = X 10 + δX 1 cos(ω t )X1

X 2 (t ) = X 20 + δX 2 cos(ω t − ϕ )Pumped Charge Current in ZGNR
Time-dependentgate Voltage

X1

X2X 1 (t ) = δX 1 cos(ω t )X 2 (t ) = δX 2 cos(ω t − ϕ )ϕEnergy [eV]

21E0-1-2-3
0.00.5kza/π

1.00.1

Pumped Charge/ePhase
difference3

δX 1 = δX 2 = 0.1 [eV]
up0.0

down
-0.1
0.0

Pumped charge (pumping current) is inverted by thechange of relative phase by pi.

0.5

1.0

φ/π

1.5

2.0

Pumped Spin Current in ZGNRcos(ω t )X1cos(ω t − ϕ1 )ϕ1 = π / 2X2

X4

X3

0.4

0.4

0.2

0.2

0.0

0.0

-0.2

-0.2

-0.4

-0.4

-1.0 -0.5 0.0 0.5 1.0-0.05 0.00 0.05

ka/πcos(ω t − ϕ 2 )ϕ 2 = −π / 2Charge/edown

Energy [eV]

Energy [eV]cos(ω t )Charge current is canceledÆ Pure spin currentup

Summary
１．強束縛分子動力学法（Tight binding molecular dynamics)を用いた構造緩和計
算のFermi Operator Expansion 法を用いた高速化の実装→
今後の計算への活用，強束縛密度汎関数法（Density-functional tight-binding ：DFTB)法のフレームワークへの拡張

２．波束の時間発展を用いた時間依存電気伝導計算プログラムの構築と，光照射されたグラフェンの電気伝導への応用．
→Floquet scattering theoryを用いた解析との比較など．

３．強束縛密度汎関数法と非平衡グリーン関数法を用いたナノ構造の電気/スピン伝導プログラム→Adiabatic Pumping現象を扱うように拡張．
→ジグザグ端グラフェンナノリボンにおける電荷，スピンポンピングの解析．
Pure Spin Current 生成を示唆