Classical Molecular Dynamics of Phonons and Electrons

Classical molecular dynamics of phonons and electrons with quantum baths Jian-Sheng WANG 1

Outline A quick introduction to molecular dynamics Langevin equation Generalized Langevin equation with quantum heat baths Application examples Quantum heat transport in nanostructure Electron transport Electron-magnet coupled system 2

Molecular Dynamics 3

Basic Idea Solve Newton s equations of motion 2 r d dt V r j = = = F F , , 1,2, , m j N j j j 2 j Choose a force field (specified by a potential V) appropriate for the given system under study Decide a statistical ensemble to use, choice of boundary conditions; collect statistics of observables 4

Commonly Use Force Fields Lennard-Jones Potential For noble gas and generic fluids Tersoff, Brenner, Stillinger-Weber, 3-, 4-body potentials For C, Si, Ge, AMBER, CHARMM, GROMOS, MM4, etc For biomolecules GULP, DFT codes, etc For solids 5

Ensembles Micro-canonical Ensemble Energy is fixed Canonical ensemble Need to use thermostat to fix temperature Langevin dynamics Nos -Hoover Generalized Langevin 6

Langevin theory of Brownian motion (1908) dv dt = m v + ( ) m F t = 2 2 x Dt 1 2 k T = 2 mv B 2 k T m k T = = D B B = ( ) ( ) F t x t 0 6 a dx dt = v P. Langevin 1872- 1946 7

Langevin Dynamics 2 r F d dt j j = + = v , 1,2, , j N j 2 m = ( ) t 0, 2 k T m = T B ( ) ( ') t ( ') t t t I 8

Observables, Statistics Equilibrium temperature (in micro-canonical ensemble) by equipartition theorem. 1 1 2 2 = = 2 j , , , x y z k T m v , B j Pressure of a fluid (for pair potential) 1 d r = + r F ( ) PV Nk T B i j ij i j Where d is dimension, Fij is the force acting on particle i from particle j. 9

Transport Coefficients The diffusion constant can be computed through velocity correlation function 2 r r ( ) t (0) = = v v ( ) t (0) lim t 3 dt D 2 t 0 10

Transport Coefficients Thermal conductivity can be computed through energy-current correlation using Green-Kubo formula; or nonequilibrium simulation by directly computing the energy current = Fourier's law: J T 11

Generalized Langevin t = + r ( ') ( ') t u t dt ' u F t = ( ) t 0, = ( ) t ( ') t ( ') i t t 1 2 ( ) is known as self-energy = + 12

Quantum noise with Langevin dynamics For thermal transport, Dhar & Roy, J Stat Phys (2006); JSW, PRL 99, 160601 (2007). Thermostat to solids, Ceriotti, Bussi, Parinello, PRL 103, 030603 (2009); Dammak, Chalopin, Laroche, Mayoun, Greffet, PRL 103, 190601 (2009). Kantorovich, et al, PRB 78, 094303 (2008); 94, 184305 (2016). 13

Thermal conduction at a junction Left Lead, TL Right Lead, TR semi-infinite Junction Part 14

Three regions 1 u u L L 2 = = = , , u u u u u C L L C u R i = = T ( , ) ' ( ) ( ) ' , , , , G T u u L C R C Tc is contour order super-operator; u is m x 15

Quantum heat-bath & MD Consider a junction system with left and right harmonic leads at equilibrium temperatures TL & TR, the Heisenberg equations of motion are = = = L LC , u K u V u u u L L C 1 C CL CR , u F V u V u C L R 2 R RC = u K u V u u R R C u The equations for leads can be solved, given j t = + 0 L LC ( ) ( ) ( ') ( ') ', u t u t g t t V u t dt L L C where 2 2 d dt d dt + = + = 0 L L L ( ) 0, ( ) ( ) t I K u t K g t L 2 2 16

Quantum Langevin equation for center Eliminating the lead variables, we get t = + + C ( ') ( ') ' u F t t u t dt C C L R where retarded self-energy and random noise terms are given as = + = = C C , , V g V = L R u 0 C , , V L R See JSW, Phys. Rev. Lett. 99, 160601 (2007); JSW, et al, Phys. Rev. B 80, 224302 (2009). (such idea for classical ensemble appears as early as in 1976, Adelman & Doll, JCP. But see Kantorovich PRB 2008) 17

Properties of the quantum noise = ( ) t 0, = CL T LC L 0 L 0 L ( ) t ( ') t ( ) ( ') V u t u t V L = = CL LC ( ') ( '), V i g t t V i t t L L T = L ( ') t ( ) t ( '), i t t L L + T i t = [ ] = ( )Im [ ] L ( ') t ( ) t 2 e dt i f L L L For nonequilibrium Green s function (NEGF) notations, see JSW, Wang, & L , Eur. Phys. J. B, 62, 381 (2008); or JSW, Agarwalla, Li, & hingna, Front. Phys. 9, 673 (2014). 18

Quasi-classical approximation (Schmid 1982) Replace operators uC & by ordinary numbers Using the symmetrized quantum correlation, i ( > + <)/2 for the correlation matrix of . For linear systems, quasi-classical approximation turns out exact! See, e.g., Dhar & Roy, J. Stat. Phys. 125, 805 (2006). 19

Implementation Generate noise using fast Fourier transform 1 ( ) k k hM lk M = = 2 / i t lh e Solve the differential equation using velocity Verlet Perform the kernel integration using a simple rectangular rule Compute energy current by t dH dt = = + T C ( ') ( ') c t u t dt ' I u t L L L L 20

Equilibrium simulation 1D linear chain (red lines exact, open circles QMD) and nonlinear quartic onsite (crosses, QMD) of 128 atoms. From Eur. Phys. J. B, 62, 381 (2008). 21

Comparison of QMD with NEGF Three-atom junction with cubic nonlinearity (FPU- ). From JSW, Wang, Zeng, PRB 74, 033408 (2006) & JSW, Wang, L , Eur. Phys. J. B, 62, 381 (2008). QMD ballistic QMD nonlinear kL=1.56 kC=1.38, t=1.8 kR=1.44 22

From ballistic to diffusive transport 1D chain with quartic onsite nonlinearity ( 4 model). The numbers indicate the length of the chains. From JSW, PRL 99, 160601 (2007). Classical, 0 4 16 NEGF, N=4 & 32 64 256 1024 S L 4096 = = , J T 23

FPU- model It is not clear whether QMD over or under estimates the nonlinear effect. From Xu, JSW, Duan, Gu, & Li, PRB 78, 224303 (2008). }relaxing rate NEGF 24

Conductance of graphene strips Sites 0 to 7 are fixed left lead and sites 28 to 35 are fixed right lead. Heat bath is applied to sites 8 to 15 at temperature TL and site 20 to 27 at TR. JSW, Ni, & Jiang, PRB 2009. 25

Zigzag (5,0) carbon nanotubes Temperature dependence of the thermal conductance (and conductivity) for lengths 4.26 (green), 12.8 (red) and 25.6 (blue) nm, respectively. JSW, Ni and Jiang, PRB 2009. 26

Electron transport & phonons For electrons in the tight-binding form interacting with phonons, the quantum Langevin equations are (set = 1) t ic Hc t t c t dt = + + + k ( ') ( ') ' , M u c k t = + ( ') ( ') t u t dt ' u Ku t c Mc T ( ) ( ') t = ( ') ( ) t = ( '), ( '), t i t t t i t t T ( ) ( ') t = = T ( '), ( ') t ( ) t ( ') t i t t i t t 27

Quasi-classical approximation & NEGF D< <n= i G< G< } D< Dr Dr rn= i { + Gr G> To lowest order in coupling M, quasi- classical approximation is to replace all G> by G<. G< <n= i G> Gr G< rn= i { } + G< Ga 28

Ballistic electron transport, NEGF vs QMD Nearest neighbor hopping model with two sites in the center, lead hopping -hl = 0.1 eV, varying the center part hopping term. From L & JSW, J. Phys.:Condens. Matter 21, 025503 (2009). For NEGF method of electron-phonon interaction, see L & JSW, PRB 76, 165418 (2007). 29

Strong electron-phonon interactions Two-center-atom model with Su, Schrieffer & Heeger electron-phonon interaction. Lines are NEGF, dots are QMD. From L & JSW, J. Phys.: Condens. Matter, 21, 025503 (2009). ballistic 30

QMD is exact in low electron density limit Self-consistent Born approximation NEGF QMD Low electron density 31

Ballistic to diffusive Electronic conductance vs center junction size L. Electron-phonon interaction strength is m=0.1 eV. From L & JSW, J. Phys.: Condens. Matter, 21, 025503 (2009). 32

Electrons as a heat bath (L, Brandbyge, et al, based on path-integral formulism) t = + C ( ') ( ') t u t dt ' u F F t e t + + ( ') ( ') t u t dt ' t e e 33

Joule heating from QMD Heat per unit time generated due to electric current for a 4 2 graphene armchair configuration. From L , et al, AIP Advances, 5, 053204 (2015). 34

A simplified QMD? Consider the following Langevin equation for lattice vibration [Keblinski & JSW, unpublished; see also Buyukdagli, et al, PRE 78, 066702 (2008); Ceriotti, et al, JCTC (2010) ]: = + ( ), t u Ku = u ( ) ( ') t = T ( ) t 0, ( ) , c t I t + 1 2 i t [ ] c = = ( ) + ( ) c t e dt 2 f then in the limit of small damping ( 0), the energy of the vibrational modes is given exactly as that of the corresponding quantum system. 35

df coupling model = + l l l S H t c c J c c B c c ' ' ' l H j l B l l l = S S , , , ' j l | J '| 1, , l + z j S S J g B S 1 2 i j i j B nn nnn j Spins sit on the corners of a cubic lattice, electron at the center of the lattice. Parameters: t 0.2 eV, JH 0.1 eV or less, lattice constant a = 3.9 , J1/kB = - 0.037 K, J2/kB = 0.069 K, g = 2. S j lc Approaches: 1) Boltzmann + RPA 2) Dynamic mean-field 3) Exact diagonalisation 36

Magnetic polaron physics? A single electron within the background sea of Heisenberg spins. We can only talk about the mobility, v= E; then = en is correct at low carrier density, say n 10-4 per unit cell. Our approach: 1) Derive Heisenberg equations of motion for cl , and Sj from the df coupling model. 2) replace Sj by a classical vector, replace cl by complex number (i.e., single particle wavefunction). 3) Simulating a single electron diffuses on lattice under the influence of fluctuating magnetic moments. 4) Determine the diffusion constant D, and use Einstein relation to get = D/(kBT). 37

Operator-splitting method , ) ( , , ) { , } H A A H L A = ( , H c c S dA dt = * H S Trotter-Suzuki formulas + = + ( ) 2 t A B tA tB e e ( ) e O t = t t A A + = + = + ( ) 3 3 t A B tB ( ) ( ) ( ) e e e e O t S t O t 2 2 2 , } = tL { t H ( ) (0) (0) A t e A e A ( ) H + = (1 2 ) + ( ) 5 t A B ( ) ( ) ( ) e S t S t S t O t 2 2 2 = 1.35120 Poisson bracket definition: See Natano & Suzuki, arXiv:math-ph/0506007 1 i A B B A { , } A B * l * l l l l 1 A S B S S j j j j 38

Electron diffusion A typical MD run result of the mean-square displacement for a 2D version of the model. Temperature T=2K, JH=0.005 eV, t=0.05 eV. Unpublished, from Ruofan Chen. 39

Conclusion & Outlook QMD for phonon is correct in the ballistic limit and high-temperature classical limit Produce correct physics for single electron Much large systems can be simulated (comparing to NEGF) Quantum corrections in dynamics? Solve a hierarchical set of Heisenberg equations (Prezhdo et al, JCP) Can we treat the noises & as operators (i.e. matrices) thus restore > <? 40

Collaborators/students Saikong Chin Chee Kwan Gan Jinghua Lan Yong Xu Baowen Li Pawel Keblinski Jian Wang Jingtao L Jin-Wu Jiang Eduardo Cuansing Nan Zeng Bijay Kumar Agarwalla Thingna Juzar Yahya Ruofan Chen Lifa Zhang Xiaoxi Ni 41