Role of direct electron-phonon coupling across metal-semiconductor interfaces in thermal transport via molecular dynamics Keng-Hua Lin and Alejandro Strachan

Citation: The Journal of Chemical Physics 143, 034703 (2015); doi: 10.1063/1.4922893 View online: http://dx.doi.org/10.1063/1.4922893 View Table of Contents: http://aip.scitation.org/toc/jcp/143/3 Published by the American Institute of Physics

Articles you may be interested in Role of electron–phonon coupling in thermal conductance of metal–nonmetal interfaces The Journal of Chemical Physics 84, (2004); 10.1063/1.1758301 The effect of the electron-phonon coupling on the effective thermal conductivity of metal-nonmetal multilayers The Journal of Chemical Physics 109, 094310094310 (2011); 10.1063/1.3585824

THE JOURNAL OF CHEMICAL PHYSICS 143, 034703 (2015)

Role of direct electron-phonon coupling across metal-semiconductor interfaces in thermal transport via molecular dynamics Keng-Hua Lin and Alejandro Strachan School of Materials Engineering and Birck Nanotechnology Center, Purdue University, West Lafayette, Indiana 47907, USA

(Received 27 January 2015; accepted 7 June 2015; published online 16 July 2015) Motivated by significant interest in metal-semiconductor and metal-insulator interfaces and superlattices for energy conversion applications, we developed a molecular dynamics-based model that captures the thermal transport role of conduction electrons in metals and heat transport across these types of interface. Key features of our model, denoted eleDID (electronic version of dynamics with implicit degrees of freedom), are the natural description of interfaces and free surfaces and the ability to control the spatial extent of electron-phonon (e-ph) coupling. Non-local e-ph coupling enables the energy of conduction electrons to be transferred directly to the semiconductor/insulator phonons (as opposed to having to first couple to the phonons in the metal). We characterize the effect of the spatial e-ph coupling range on interface resistance by simulating heat transport through a metal-semiconductor interface to mimic the conditions of ultrafast laser heating experiments. Direct energy transfer from the conduction electrons to the semiconductor phonons not only decreases interfacial resistance but also increases the ballistic transport behavior in the semiconductor layer. These results provide new insight for experiments designed to characterize e-ph coupling and thermal transport at the metal-semiconductor/insulator interfaces. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4922893]

I. INTRODUCTION

Thermal transport across metal-semiconductor and metalinsulator interfaces is of great interest due to the ubiquitous nature of these structures in solid-state devices,1 microelectromechanical systems (MEMS),2 and thermoelectrics.3,4 For example, metal-semiconductor superlattices have been proposed for thermoelectric applications since a high figure of merit, ZT, could be achieved by varying superlattice periods to tune phonon thermal conductivity without affecting the electronic properties.3,4 Furthermore, a recent device-level sensitivity analysis shows that the effect of thermal interface resistance between the semiconducting thermoelectric materials and the metal contacts on ZT increases more than 50% as the heat flux increases from 100 to 500 W/cm2.5 More generally, interfacial properties have become increasingly important with device miniaturization due to increasing surface to volume ratios; therefore, understanding the fundamentals of thermal transport across metal-semiconductor/insulator interfaces is necessary to advance the design of nanostructured materials and devices. An interesting and important aspect of thermal transport across the metal-semiconductor/insulator interfaces is the change in the dominant heat carriers. Phonons govern transport in undoped semiconductors and insulators, while electrons are the main heat carriers in metals; therefore, energy exchange between electrons and phonons must occur at the metal-semiconductor/insulator interfaces. There are two possible pathways for energy to transfer from the metal to the semiconductor/insulator layer across the interface:6 (1) the conduction electrons couple to the metallic atoms, and 0021-9606/2015/143(3)/034703/9/$30.00

then phonons carry the energy across the interface or (2) the electrons in the metal couple directly with the phonons in the semiconductor/insulator, providing an additional path across the interface. Evidence of heat transfer via pathway (2) and estimates of the resistance associated with direct electron-phonon (e-ph) heat transfer across the interface have been obtained in the recent experimental studies.4,7 In addition, theoretical studies have shown that this pathway is non-negligible when e-ph coupling is strong or when the phonon interface resistance is high by calculating the thermal conductance from the phonon emission/absorption rates and the e-ph scattering matrices.8–10 These studies suggest that the effect of spatially extended eph coupling on thermal transport becomes more significant when the electrons and the phonons in the metals are under non-equilibrium conditions. This occurs during the initial equilibration stage of the ultrafast laser heating in the transient thermoreflectance measurements7 in which electrons obtain a significant amount of energy from the laser pulses or during heat transport in the metal-semiconductor superlattices with periods of a few nanometers,4 shorter than the distance required for electrons and phonons to achieve thermal equilibrium. Energy transport in metal-semiconductor systems has been studied by analytically solving the non-equilibrium Green’s function (NEGF),11–13 Boltzmann transport equation (BTE),14 two-temperature model (TTM),6,15–18 and by molecular dynamics (MD) combined with TTM (MD-TTM).19 NEGF and BTE methods are often combined with density functional theory calculations of force constant for interfaces or lowdimensional materials and require a model to treat scattering.14 TTM approaches use two coupled heat diffusion equations

143, 034703-1

© 2015 AIP Publishing LLC

034703-2

K.-H. Lin and A. Strachan

to describe the energies and temperatures of the electronic and the atomic subsystems and consider the energy exchange between them. Since the heat diffusion equations are solved in TTM, ballistic transport, which is important in nanoscale specimens, is not captured. When combined with MD, the heat diffusion equation of the atomic subsystem is replaced by explicit atomistic dynamics so that the real atomistic structure and phonon properties are captured together with any possible ballistic or coherent transport phenomena. However, prior MD-TTM efforts have assumed that heat transfer follows only pathway (1). In summary, despite significant experimental and theoretical progress in the understanding of thermal interface resistance of metal-semiconductor/insulator interfaces, important questions remain to be answered. Specifically, how does direct e-ph coupling across the interface affect transport? No model to date both captures direct e-ph coupling across the interface and includes the non-diffusive characteristics and possible lack of local equilibrium of ionic thermal transport, and the goal of this study is to fill this gap and characterize the nature of thermal transport across metalsemiconductor interfaces. We have previously extended the model of dynamics with implicit degrees of freedom (DOFs) (DID)20 to describe electronic heat transport in metals via MD simulations coupled to the TTM.21 In this study, we introduce a key extension to our model, called eleDID (electronic version of DID), to account for the effect of spatially extended electron-phonon coupling. EleDID differs from the prior MD-TTM efforts in two aspects: (i) the temperatures of the electronic DoFs are described using the positions of the metallic atoms as compared to the separate, rigid spatial grid used by MD-TTM; therefore complex atomistic structures can be captured naturally; (ii) the conduction electrons in the metal can exchange energy with the atoms in a non-local manner, i.e., the electrons around an atom exchange energy with the atoms in their vicinity regardless of their nature (metallic or non-metallic). The spatial extent of this non-local e-ph coupling is the only new adjustable parameter in our model. In this paper, we focus on the effect of the spatial extent of e-ph coupling on thermal interface resistances between a metal and a semiconductor, mimicking ultrafast laser heating in the transient thermoreflectance measurements.7 Interestingly, we find that an e-ph coupling extent of approximately two lattice parameters results in an interfacial thermal resistance between the electrons in the metal and the phonons in the semiconductor comparable to the recent experimental measurements. Also, our results show that direct e-ph coupling accounts for about 50% of the total energy transferred across the interface in the early stages of pulsed laser heating experiments when electrons and phonons in the metal are away from local equilibrium. In addition, we find that the rapid energy transfer across the interface due to direct e-ph coupling enhances ballistic thermal transport of the high-frequency phonons in the semiconductor. II. TWO TEMPERATURE MODEL MD SIMULATIONS WITH NON-LOCAL ELECTRON-PHONON COUPLING

DID describes particle dynamics using classical equations of motion and describes the additional DoFs associated with the particles in an averaged way via an additional

J. Chem. Phys. 143, 034703 (2015)

dynamical variable(s), which couple to particle dynamics to enable the interaction between the explicit and implicit DoFs. In coarse-grained or mesodynamics simulations, particles represent groups of atoms, and the DoFs internal to the particles are described implicitly. Such mesoDID approach has been applied to shockwave propagation,20,21 dynamical failure,22 and thermal transport23 and was recently extended to describe chemical reactions associated with implicit DoFs.24 DID can also be generalized to describe the thermal transport role of conduction electrons in metals (eleDID) using the twotemperature model; in this case, particles represent atoms, and an additional dynamical variable, Tielec, is assigned to each atom to represent the electronic temperature of the electrons in the vicinity of atom i 21 as pioneered by Hakkinen and Landman.25 In this paper, we generalize eleDID to enable electrons and atoms to couple in a spatially delocalized manner and apply the method to characterize transient transport across the metal-semiconductor interface. The remainder of the section describes eleDID focusing on these new features. As discussed above, the thermal state of the electronic DoFs is described with a temperature associated with each atom in the system. In order to couple the electronic and the atomic subsystems, we define the local atomic temperature in the vicinity of every atom in the system, Tiatom. Tiatom is obtained from the velocity fluctuations around the average velocity of each atom and its neighboring atoms. As in Refs. 20 and 21, we define the average velocity in the vicinity of atom i as N 

⟨u⟩i =

j=1

w(r i j )m j u j

N  j=1

,

(1)

w(r i j )m j

where m j is the mass of atom j and w(r i j ) is a properly normalized weighting function describing the contributions of atom j to the local average velocity around atom i based on their separation distance, r i j . The local atomic temperature around atom i is dk BTiatom =

N 

 2 w(r i j )m j u j − ⟨u⟩i ,

(2)

j=1

where d is the dimension of the simulation and k B is the Boltzmann’s constant. As in prior work,21 we use a function, w(r i j ), which smoothly decays to zero at a cutoff distance, R, as the weighting function, which has the following form:  f (r i j )4 1    · , if r i j ≤ R,  W 1 + f (r i j )4 (3) w(r i j ) =      0, if r i j > R,  where f (r i j ) = (r i j − R)/σ and r i j is the interatomic distance between the atoms i and j. W is the normalization factor for the weighting function. In the simulations below, we will use R = 8 Å and σ = 2.67 Å. The weighting function and its parameters used here are similar to those used in the past21 and the results are not sensitive to them. Energy exchange between the atoms and the conduction electrons is driven by the differences between the electronic and the atomic temperature fields. Following DID, energy

034703-3

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

exchange is done through a viscous term in the position update equation which naturally leads to the Galilean invariance.20,21 An additional displacement is imposed on each atom in the direction of its total force with a magnitude and a sign depending on the difference in local temperatures,    T jatom − T jelec  ∗ + * · ν · w (r i j ) r˙i = ui + Fi ·  θ0 i ∈all   j ∈metal , 1

, (4) · mi · ω2i Fi u˙i = , (5) mi where r i and ui are the position and velocity of atom i, Fi is the

force  exerted on atom i by its neighboring atoms, mi and ω2i are, respectively, the atomic mass of atom i and the mean-square frequency of this atom type (obtained from the vibrational density of states). θ 0 is the reference temperature. w ∗(r i j ) is a normalized, local weighting function that describes the spatial extent of e-ph coupling; the function depends on the atomic separation distance, r i j , and monotonically decays to zero at a cutoff distance. Finally, ν is the coupling rate between the electronic and the atomic subsystems. The energy lost or gained by the atomic subsystem needs to be transferred to the electrons, and energy conservation requires the time derivative of the electronic energy as follows: T atom − T jelec elec ˙ elec +·ν * j E˙ elec = C · T = j j j 0 θ j ∈metal ,  |Fi |2

2 · w ∗(r i j ) + κ elec · ∇2T jelec, · m · ωi i i ∈all

(6)

where the first term on the right hand side represents the energy exchanged with the atoms (and is needed for energy conservation) and the second term describes diffusive electronic thermal is the electronic specific heat and κ elec is the transport. C elec j electronic thermal conductivity. Equation (6) in conjunction with the equations of motion (4-5) leads to energy conservation as demonstrated in the supplementary material.26 We note that eleDID and similar methods21,25 are based on the definition of local thermodynamic properties; this is a common practice in atomistic simulations and is valid as long as the properties of interest vary slowly within atomic distances. We now explain the key features of eleDID. The first term on the right hand side in Eq. (6) describes the energy transferred between the conduction electrons and the atoms. Note that all atoms (either metallic or non-metallic) are affected as long as they are near a metallic one. The electrons in the vicinity of metallic atom j, described by T jelec, exchange an amount of energy with the atoms that depend on the local temperature difference between the electrons and the atoms. If T jatom is larger than T jelec, the second term on the right hand side in Eq. (4) will add an extra displacement to atom i in the direction of its force so that energy will be transferred from the atomic subsystem to the electronic subsystem. The reverse occurs if the electrons are locally hotter. We set the reference temperature θ 0 as T jatom. Since T jatom is proportional to the long |Fi |2 time average of 2 as shown in Ref. 21, this choice i ∈all m i ·⟨ω i ⟩

makes the amount of energy exchange proportional to the nonequilibrium degree between the subsystems, (T jatom-T jelec), and an temperature-independent coupling constant (the e-ph coupling factor, G, used by MD-TTM19). In reality, G is typically temperature dependent;27 this effect can be considered by our model by making ν a temperature dependent variable but we will ignore this in the present simulations. The Laplacian term on the right hand side in Eq. (6) accounts for thermal diffusion among the conduction electrons and depends on the electronic thermal conductivity, κ elec. Since the time scale for electron relaxation is short as compared to that of atomic vibrations, an electronic time step smaller than a MD time step should be used to evolve the electronic temperature. In addition to ν and (T jatom-T jelec), the energy exchanged between the subsystems also depends on w ∗(r i j ), which controls the spatial extent of the coupling. Local e-ph coupling is achieved by setting w ∗(r i j ) as   1, if i = j, w ∗(r i j ) =   0, if i , j. 

(7)

To capture non-local e-ph coupling and, consequently, direct e-ph coupling across the interfaces of interest here,4,7 the energy exchange range between the electronic and the atomic subsystems among the atoms is extended by using the same functional form as Eq. (3) for the local atomic properties but with different parameters R∗ and σ∗. The details of e-ph coupling across metal-semiconductor interfaces, including the coupling rate and the spatial extent, could be calculated using ab initio electronic structure calculations28–30 for specific materials combinations. In the initial studies presented here, we assume the coupling occurs within two unit cells into the semiconductor; we take R∗ = 8 Å and σ∗ = 2.67 Å. We will show that this value leads to e-ph resistance across the interface consistent with the experimental measurements. We note that the extent of the electron-phonon coupling in eleDID (local or non-local) does not affect the overall time scales for equilibration between conduction electrons and metallic ions;26 the effect of non-local coupling manifests itself strongly in the case of metal-insulator interfaces. III. SIMULATION DETAILS AND ANALYSIS A. Model metal/semiconductor system

To study the effect of e-ph coupling range on thermal resistance of metal-semiconductor interfaces, we simulate ultrafast laser heating of a system consisting of a thin Al layer over a fictitious semiconductor material that we will denote Al∗. The metallic layer consists of 1000 atoms (5 × 5 × 10 unit cells oriented along the [001] direction) and the Al∗ substrate contains 3000 atoms (5 × 5 × 30 unit cells); see Fig. 1(a). In addition to the properties mentioned above, the Al∗ system does not have electronic DoFs. The atomic structures of Al and Al∗ are both fcc and their interatomic interactions are described with an embedded atom method interatomic potential.31 To generate phonon thermal resistance between the two materials, Al∗ is assigned with a lighter mass, 5.611 g/mol, as compared to 28.0855 g/mol for Al, which creates acoustic mismatch between the two materials similar to that between Au and Si.

034703-4

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

FIG. 1. (a) Structure of the metal-semiconductor system. The blue lines represent the simulation box. (b) Temporal evolution of the atomic and the electronic temperatures in the 3rd, 2nd, and 1st bins in the metal and the semiconductor layers from the interface. The solid red and blue lines represent the atomic temperature for the local and non-local coupling cases, respectively; the dashed red and blue lines represent the electronic temperatures in the metal layers for the local and the non-local coupling cases, respectively.

This choice is motivated by the availability of the experimental data for Au/Si in Refs. 7 and 32. We stress that our aim is not to model a specific metal/semiconductor interface but to study the physics of energy transport at such interfaces that necessarily involve a change in heat carriers. Modeling realistic material interfaces would require force fields capable of describing both materials and their interactions; in addition, different crystal structures and/or lattice parameters would lead to interfacial defects that affect phonon heat transport.

ν = 0.0057 ps−1, which leads to G ∼ 4.3 × 1016 W(m3 K)−1 as shown in the supplementary material.26 For simulations with local e-ph coupling, Eq. (7) is used as the weighting function, and for the non-local e-ph coupling case, we use Eq. (3) with R∗ = 8 Å (just below two lattice parameters) and σ ∗ = 2.67 Å. The atomic time step for MD integration is 0.1 fs, and the electronic time step for solving the electronic temperatures by diffusive heat transport is 0.0025 fs. C. Simulating transient laser heating

B. eleDID parameterization

The remaining parameters for the eleDID simulations are obtained from experiments. We use a temperature dependent electronic specific heat C elec = 132.63T elec J(m3 K)−1 33 and thermal conductivity κ elec = 221.752 W(m K)−1, which is the total thermal conductivity of pure Al,34 since the phonon contribution to thermal conductivity

 at 300 K in Al is less than 10 W(m K)−1.35 The frequency ω2 = 65.30 ps−2 is calculated from the vibrational density of states obtained from a MD simulation using our force field. The coupling constant ν is obtained from the e-ph coupling factor, G, which for metals varies in the range of 1–10 × 1016 W(m3 K)−1.36,37 Therefore, we take

Before simulating the heating due to a laser pulse, the metal-semiconductor material is equilibrated via a three-step procedure: (i) the system is relaxed under isobaric, isothermal (NPT) conditions at 300 K for 50 ps under 3-D periodic boundary conditions and (ii) the system is relaxed under isochoric, isothermal conditions (NVT) at 300 K for 60 ps with the equilibrium lattice parameter calculated from 40 to 50 ps from the previous NPT run. During NVT runs, the velocities of the atoms are reassigned to 300 K every 10 ps, and periodic boundary conditions are applied only along the x and y directions, and open boundary condition is imposed along z. Re-assigning atomic velocities periodically helps dissipate

034703-5

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

the long wavelength oscillations of the structure caused by the relaxation at the free surfaces. (iii) Finally, we perform a 50 ps-long NVE simulation to equilibrate the atomic and the electronic subsystems at 300 K. After the structural relaxation, ultrafast laser heating7 is simulated by assigning an electronic temperature of 5000 K to all the Al atoms and using the eleDID equations to simulate the dynamics of the system. IV. ENERGY TRANSPORT ACROSS THE INTERFACE AND INTERFACE RESISTANCE

To understand the details of thermal transport in metalsemiconductor systems, we compute the temporal evolution of the local atomic and electronic temperatures and the associated energies. As discussed in Subsections IV A–IV C , the thermal interface resistance between the electrons in the metal and the phonons in the semiconductor and that between the phonons in the metal and in the semiconductor are calculated from the amounts of the energy transferred across the interface and the temperature jumps at the interface. A. Temporal evolution of the local temperatures

We divide the material into 20 bins along the z direction to investigate the evolution of the local temperatures as a function of position and time. Figure 1(b) shows the temporal evolution of the local temperatures in the first three bins of the metal and the semiconductor layers near the interface. We observe a fast equilibration between the electronic and the atomic temperatures in the metal and a progressive heating of the semiconductor layer. The equilibration time between the electronic and

∆Ee−ph (t) =

the atomic temperatures for the metal bins other than the one adjacent to the interface is about 1–2 ps, which is comparable to the experimental data.7,37 For the metallic bin adjacent to the interface, energy transfer to the neighboring semiconductor bin causes a longer equilibration time between the electronic and the atomic subsystems within the bin. Clearly, direct coupling between the conduction electrons in the metal and the phonons in the semiconductor (non-local e-ph coupling) leads to a faster increase in atomic temperature in the semiconductor region adjacent to the interface. Interestingly, the atomic temperature in the first semiconductor bin overshoots the final temperature for both the local and the non-local coupling cases and then cools down as the heat flows into the rest of the semiconductor. B. Energy transfer across the interface

The total energy increase in the semiconductor, ∆ES, has two origins: (i) the direct coupling between the electrons in the metal and the phonons in the semiconductor, ∆Ee-ph, and (ii) the phonon-phonon coupling between the metal and the semiconductor atoms across the interface, ∆Eph-ph, ∆ES = ∆Ee−ph + ∆Eph−ph.

∆ES at time t can be obtained from the increase in potential energy, Ep, and kinetic energy, Ek, of the semiconductor layer, ∆ES (t) = E p (t) − E p (t = 0) + Ek (t) − Ek (t = 0). (9) ∆Ee-ph at time t can be obtained as the sum of the energy exchange between the electrons in the metal and the phonons in the semiconductor at each MD time step as

 t     T jatom − T jelec |Fi |2 ∗  · dt, +·ν· *

 · w (r ) i j 0 2 θ   m · ω i i i ∈semiconductor t=0 j ∈metal , -

where dt is the MD atomic time step. The temporal evolutions of ∆ES and ∆Ee-ph using local and non-local couplings are shown in Fig. 2(a). Direct energy transfer between electrons and phonons across the interface leads to a faster rise in the energy of the semiconductor layer; furthermore, the results show that for short time scales, ∆ES is dominated by the direct e-ph coupling process (∆Ee-ph) across the interface when non-local coupling is enabled. Even after 4 ps, direct e-ph coupling across the interface accounts for about half of the total energy transferred to the semiconductor layer for the non-local coupling case. These results lend credence to the recent experimental efforts suggesting that ultrafast laser heating experiments are ideally suited to measure e-ph coupling across the metal-semiconductor interfaces. We clarify that, as mentioned above, the total energy of the system is conserved during this energy transfer process as shown in the supplementary material.26 C. Thermal interface resistance

Figure 2(b) shows the temporal evolution of e-ph and phonon-phonon thermal interface resistance with local and

(8)

(10)

non-local couplings. The thermal interface resistance between the electrons in the metal layer and the atoms in the semiconductor layer is calculated as Re−ph =

elec TM − TS , e→ J ph

(11)

where TM elec and TS are the electronic temperature of the metallic bin and the atomic temperature of the semiconductor bin adjacent to the interface. J e→ ph is the heat flux from the electronic subsystem in the metal to the atomic subsystem in the semiconductor and is equal to dEe−ph/(A · dt). A is the cross sectional area of the material. This definition of thermal interface resistance between the electrons in the metal layer and the atoms in the semiconductor is commonly applied to study interfacial phenomena.6,17,19 Thermal interface resistance between the atomic subsystems in the metal layer and in the semiconductor layer is calculated as atom TM − TS , (12) ph→ J ph where TM atom is the atomic temperature of the metal bin adjacent to the interface. J ph→ ph is the heat flux transferring from

Rph−ph =

034703-6

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

FIG. 2. (a) Increase of total energy in the semiconductor layer, ∆Es, and the energy obtained from direct electron-phonon coupling across the interface, ∆Ee-ph, for the local and non-local coupling cases. (b) Temporal evolution of the electron-phonon and the phonon-phonon thermal interface resistances, Re-ph and Rph-ph, with electron-phonon coupling rate of 0.0057 ps−1 for the local and the non-local coupling cases.

the atomic subsystem in the metal layer to the atomic subsystem in the semiconductor and is equal to dE ph-ph/(A · dt). ∆Eph-ph is equal to (∆ES − ∆Ee-ph) as stated in Eq. (8). Since ∆Eph-ph fluctuates over time, J ph→ ph at time t is calculated from a linear fit of the instantaneous ∆Eph-ph versus t curve in the time interval of ±0.1 ps around t. As shown in Fig. 2(b), Rph-ph obtained from the local and non-local coupling simulations are approximately 1 × 10−9 m2 K/W. When non-local coupling is allowed, Re-ph increases from 0.25 × 10−8 m2 K/W during the early stages of heat transfer process when the metallic electrons are very hot to a steady-state value of ∼1.2 × 10−8 m2 K/W after the atomic and electronic subsystems in the metal achieve equilibrium. Interestingly, using a nominal parameter for e-ph coupling and assuming the interaction between the electrons and the semiconductor atoms extends over approximately two lattice parameters, the values of Re-ph obtained from our simulation are comparable to the results from the transient thermoreflectance measurements of Au/Si by Guo et al. (Re-ph = 0.1–1 × 10−9 m2 K/W)7 and by Hopkins and Norris (Re-ph = 0.5–10 × 10−9 m2 K/W).32 Previous studies also show that metal-semiconductor systems with rather different acoustic mismatches possess Re-ph values within the same range.9,32 An interesting result from our simulation is that the smaller the non-equilibrium degree between electrons and phonons in the metal during the early stage of the equilibration process, the smaller the heat flux J e→ ph, and therefore the larger the value of Re-ph; this result is in agreement with the observations by Hopkins and Norris,32 but different from the assumption by Guo et al.,7 who has assumed Re-ph to be a constant in their description of the experimental data. During the initial stage of the process, when the electrons and phonons in the metal are away from equilibrium, the value of Re-ph varies between 1 and 10 times Rph-ph; this small difference between Re-ph and Rph-ph indicates that direct e-ph coupling across the metal-semiconductor interface provides a non-negligible energy transfer pathway in this stage. At time intervals between 2 and 5 ps, the electronic and the atomic subsystems in the metal layer have achieved local equilibrium except the bins

right next to the interface. The observed Re-ph in this period shows a steady-state value, which indicates that J e→ ph and (TM atom-TS ) are in a linear response region. After 5 ps, Re-ph can no longer be defined by Eq. (11) since both the heat flux and temperature difference are close to zero.

V. MODE-RESOLVED ENERGY TRANSFER ACROSS INTERFACE AND PHONON THERMAL TRANSPORT

To understand how energy in the semiconductor is distributed among phonon modes during the energy transfer process, we calculate the velocity power spectrum, which quantifies the spectral content of the kinetic energy, from the Fourier transform of atomic velocities by the method stated in Ref. 38. The power spectrum is obtained as  3N  2 1 τ m j · lim P(ν) = dt · exp(−i2π f t) · u j (t) , (13) τ→ ∞ 2τ −τ j=1 where f is the frequency of the phonon mode and the power spectrum is calculated within a time interval of ±τ. Under equilibrium conditions and for a classical (i.e., not quantum) system, every phonon mode possesses the same amount of kinetic energy: 21 k BT; therefore, the vibrational power spectrum is proportional to the vibrational density of states, DoS. In contrast, the kinetic energy may distribute unevenly among the phonon modes under non-equilibrium conditions. To understand how energy is deposited into the semiconductor and then carried through it, the systems are divided into 20 bins along the z direction and the vibrational power spectrum of each bin is calculated as a function of time using a time interval of 5 ps to compute the Fourier transform. Due to the short time interval required to achieve the desired time resolution, the spectra are averaged over 10 independent runs in order to obtain better statistical results. To quantify where the extra energy in the semiconductor goes, we track the increase in kinetic energy in the semiconductor bins by computing (E non−eq(ω,t) − E eq(300K)(ω)), where E non−eq(ω,t) is the power spectrum during the ultrafast laser

034703-7

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

simulations and depends on phonon frequency, ω, and time, t, and E eq(300K)(ω) is the equilibrium power spectrum at 300 K. The contour plots in Fig. 3 show the increase in kinetic energy of the phonon modes obtained from the simulations with local and non-local e-ph coupling. We find that at the early stages, transport is dominated by the frequency domains with a high DoS (in the range of 300–800 cm−1). It is also clear that thermal transport is frequency-dependent that the modes with the frequencies around 500 cm−1 have higher diffusivity. As expected, non-local coupling leads to faster energy transfer to the semiconductor as can be clearly seen in the energy maps at 2.5 ps and 7.5 ps. This phenomenon again indicates that direct e-ph coupling is a non-negligible energy transport pathway across the metal-semiconductor interface. We now focus our attention on how fast the energy propagates across the semiconductor as a function of frequency. The symbols in Figs. 4(a) and 4(e) show total kinetic energy profiles in the semiconductor at different time intervals for the local and non-local coupling simulations, respectively. Symbols in Figs. 4(b)-4(d) and 4(f)-4(h) show kinetic energy profiles divided into three frequency ranges at different time intervals for local and non-local couplings, respectively. In order to assess the nature of transport in the semiconductor, we compare the MD results with energy profiles obtained by numerically solving the 1-D heat diffusion equation, κ E˙ = · ∇2 E + q, C

(14)

where E is the kinetic energy over a range of phonon frequencies, κ is the thermal conductivity of these phonon modes, and q is the heat flux from the metal layer as a boundary condition imposed on the left side of the simulation domain. C is the phonon specific heat and is equal to 3NkB, N being the number of phonon modes. q is obtained from fitting the MD data from 0 to 15 ps to match the kinetic energy of the semiconductor

bin near the interface. κ is adjusted to fit the MD profiles. The numerical solutions to the 1-D diffusion equation are shown as solid curves in Figure 4. Interestingly, the extracted values of the effective thermal conductivity on the thin film are different for the local and the non-local coupling simulations: κ = 36.8 W(mK)−1 and κ = 24.5 W(mK)−1, respectively. This indicates that the manner in which the energy is inputted into the semiconductor affects its response; this has also been observed in other heterogeneous systems involving interfaces.39,40 Another noteworthy result from the fittings is that the actual MD results cannot be described by a diffusive model; we observe that the diffusive model underestimates the energy far away from the interface, indicating the ballistic transport characteristic in the MD simulations.41 Interestingly, the discrepancy is larger for the nonlocal coupling case, indicating that direct, fast input of energy into the semiconductor by the electrons in the metal leads to enhanced ballistic transport. We now focus on understanding whether different phonon modes contribute to the overall thermal transport differently. Phonons are divided into three frequency ranges based on the difference in DoS between the metal and the semiconductor bins near the interface at 300 K as shown in Fig. 5: (1) the lowfrequency range (0–300 cm−1), where the number of phonons in the metal significantly exceeds that in the semiconductor, (2) the medium-frequency range (300–600 cm−1), where both the metal and the semiconductor have substantial numbers of phonons but with very different frequency dependence, and (3) the high-frequency range (600–1200 cm−1), where the metal has few phonons as compared to the semiconductor. The MD results and the numerical solutions of the 1-D diffusion equation for the three phonon ranges for the local coupling case are shown in Figs. 4(b)-4(d), and the results for the nonlocal coupling case are shown in Figs. 4(f)-4(h). As shown in Figs. 4(b)-4(d), the 1-D diffusion equation accurately describes

FIG. 3. Energy distributions among the phonon modes in the semiconductor layer for the (a) local and (b) non-local coupling cases at 2.5, 7.5, 12.5, and 37.5 ps.

034703-8

K.-H. Lin and A. Strachan

J. Chem. Phys. 143, 034703 (2015)

FIG. 4. ((a)-(d)) Kinetic energy profiles along the semiconductor layer at various times for the simulations with local coupling; (a) total kinetic energy; ((b)-(d)) frequency resolved kinetic energy in the 0-300, 300-600, and 600-1200 cm−1 ranges, respectively. ((e)-(h)) Kinetic energy profiles along the semiconductor layer at various times for the simulations with non-local coupling; (e) total kinetic energy; ((f)-(h)) frequency resolved kinetic energy in the 0-300, 300-600, and 600-1200 cm−1 ranges, respectively.

the local coupling MD results for the medium-to the highfrequency phonons, and this indicates the diffusive nature of thermal transport via these phonons. The low-frequency phonons exhibit more ballistic characteristic and deviate from the diffusive model. In contrast, when direct e-ph energy transfer is enabled across the interface, all phonon frequency ranges show some degrees of ballistic behavior and the numerical solutions to the diffusion equation are unable to fit the MD results, as shown in Figs. 4(f)-4(h). We attribute the difference in thermal transport characteristics in the semiconductor to the fact that the high-frequency phonons can get energy directly from the electrons for the nonlocal coupling case. As compared to the local coupling case,

the high-frequency phonons for the non-local coupling case obtain energy in a shorter time scale since the energy can be transferred from the electrons in the metal to the phonons in the semiconductor directly. The high-frequency phonons for the local coupling case can only obtain energy from the phononphonon scattering process among the low- and the mediumfrequency phonons, and this process leads to the overall diffusive transport behavior. This observation indicates that simple diffusive transport may be inappropriate to describe thermal transport in the semiconductor layer during the equilibration process.

VI. CONCLUSIONS

FIG. 5. Vibrational density of states of the first metallic and the first semiconductor bin on either sides of the interface at 300 K.

We extended the eleDID method by making the coupling between electrons and phonons non-local; this is key to enable direct energy transport between the electrons and phonons across metal-semiconductor/insulator interfaces. The method describes phonons via MD and, thus, makes no approximations regarding phonon transport except those intrinsic in the description of atomic interactions and the use of classical, not quantum, dynamics. Our simulations in a model system show that non-local electron-phonon coupling accounts for about 50% of the energy transferred across the metalsemiconductor interface in the early stages of pulsed laser heating experiments, greatly reducing the interface resistance. This result is in good agreement with the recent experimental transient thermoreflectance measurements. Moreover, instead of being a constant, the effective thermal interface resistance is lower at early stages when the electrons in the metal are hotter than the phonons. Most importantly, our simulations suggest that non-local electron-phonon coupling across the

034703-9

K.-H. Lin and A. Strachan

metal-semiconductor interface enhances ballistic transport of the high-frequency phonons in the semiconductor. Therefore, the details of electron-phonon coupling, including the coupling rate, coupling distance, and coupling strength between different electron and phonon modes, across the metalsemiconductor interface are worth understanding. ACKNOWLEDGMENTS

This work was supported by the US National Science Foundation through Grant Nos. EEC-0634750 (Network for Computational Nanotechnology) and EECS-1028667. Computational resources from nanoHUB.org are gratefully acknowledged. 1D.

G. Cahill, K. Goodson, and A. Majumdar, J. Heat Transfer 124, 223 (2002). 2S. M. Spearing, Acta Mater. 48, 179 (2000). 3V. Rawat, Y. K. Koh, D. G. Cahill, and T. D. Sands, J. Appl. Phys. 105, 024909 (2009). 4Z. Li, S. Tan, E. Bozorg-Grayeli, T. Kodama, M. Asheghi, G. Delgado, M. Panzer, A. Pokrovsky, D. Wack, and K. E. Goodson, Nano Lett. 12, 3121 (2012). 5S. Harsha Choday and K. Roy, J. Appl. Phys. 113, 214906 (2013). 6A. Majumdar and P. Reddy, Appl. Phys. Lett. 84, 4768 (2004). 7L. Guo, S. L. Hodson, T. S. Fisher, and X. Xu, J. Heat Transfer 134, 042402 (2012). 8J. Ren and J.-X. Zhu, Phys. Rev. B 87, 241412 (2013). 9A. V. Sergeev, Phys. Rev. B 58, R10199 (1998). 10M. L. Huberman and A. W. Overhauser, Phys. Rev. B 50, 2865 (1994). 11L. Zhang, J.-T. Lü, J.-S. Wang, and B. Li, J. Phys.: Condens. Matter 25, 445801 (2013). 12B. K. Nikoli´ c, K. K. Saha, T. Markussen, and K. S. Thygesen, J. Comput. Electron. 11, 78 (2012). 13Z. Huang, T. Fisher, and J. Murthy, J. Appl. Phys. 109, 074305 (2011). 14B. Saha, T. D. Sands, and U. V. Waghmare, J. Appl. Phys. 109, 083717 (2011).

J. Chem. Phys. 143, 034703 (2015) 15Y.

S. Ju, M.-T. Hung, and T. Usui, J. Heat Transfer 128, 919 (2006). S. Ju, J. Heat Transfer 127, 1400 (2005). 17J. Ordonez-Miranda, J. J. Alvarado-Gil, and R. Yang, J. Appl. Phys. 109, 094310 (2011). 18K.-H. Yoo and A. C. Anderson, J. Low Temp. Phys. 63, 269 (1986). 19Y. Wang, X. Ruan, and A. K. Roy, Phys. Rev. B 85, 205311 (2012). 20A. Strachan and B. Holian, Phys. Rev. Lett. 94, 014301 (2005). 21K.-H. Lin, B. L. Holian, T. C. Germann, and A. Strachan, J. Chem. Phys. 141, 064107 (2014). 22K. Lynch, A. Thompson, and A. Strachan, Modell. Simul. Mater. Sci. Eng. 17, 015007 (2009). 23Y. Zhou and A. Strachan, J. Chem. Phys. 131, 234113 (2009). 24E. Antillon, K. Banlusan, and A. Strachan, Modell. Simul. Mater. Sci. Eng. 22, 025027 (2014). 25H. Hakkinen and U. Landman, Phys. Rev. Lett. 71, 1023 (1993). 26See supplementary material at http://dx.doi.org/10.1063/1.4922893 for a definition of the electron-phonon coupling factor. 27Z. Lin, L. Zhigilei, and V. Celli, Phys. Rev. B 77, 075133 (2008). 28G. Q. Huang, New J. Phys. 13, 093023 (2011). 29P. Hofmann, I. Y. Sklyadneva, E. D. L. Rienks, and E. V. Chulkov, New J. Phys. 11, 125005 (2009). 30K. Fukutani, H. Hayashi, I. N. Yakovkin, T. Habuchi, D. Hirayama, J. Jiang, H. Iwasawa, K. Shimada, Y. B. Losovyj, and P. A. Dowben, Phys. Rev. B 86, 205432 (2012). 31Y. Mishin, M. J. Mehl, and D. A. Papaconstantopoulos, Phys. Rev. B 65, 224114 (2002). 32P. E. Hopkins and P. M. Norris, Appl. Surf. Sci. 253, 6289 (2007). 33D. A. Dicke and B. A. Green, Phys. Rev. 153, 800 (1967). 34W. D. Callister, Materials Science and Engineering–An Introduction, 6th ed. (John Wiley and Sons, Inc., New York, 2003), p. 755. 35H. N. De Lang, H. van Kempen, and P. Wyder, J. Phys. F: Met. Phys. 8, L39 (1978). 36P. E. Hopkins, J. L. Kassebaum, and P. M. Norris, J. Appl. Phys. 105, 023710 (2009). 37J. L. Hostetler, A. N. Smith, D. M. Czajkowsky, and P. M. Norris, Appl. Opt. 38, 3614 (1999). 38A. Strachan, J. Chem. Phys. 120, 1 (2004). 39Y. Zhou and A. Strachan, J. Chem. Phys. 138, 124704 (2013). 40L. Hu, T. Desai, and P. Keblinski, Phys. Rev. B 83, 195423 (2011). 41B. Vermeersch, A. M. S. Mohammed, G. Pernot, Y. R. Koh, and A. Shakouri, Phys. Rev. B 90, 014306 (2014). 16Y.