Roles of Atomic Restructuring in Interfacial Phonon Transport

PHYSICAL REVIEW B 82, 081302R 2010

Roles of atomic restructuring in interfacial phonon transport

 Seungha Shin and Massoud Kaviany*
Department of Mechanical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2125, USA

Tapan Desai and Richard Bonner
Advanced Cooling Technologies, Inc., Lancaster, Pennsylvania 17601, USA
(Received 17 July 2010; published 5 August 2010)

Phonon resistance can dominate the interfacial thermal resistance between hard and soft solids. Using ab initio calculation, molecular-dynamics simulation and the diffuse mismatch model for Si/In as example, we decompose phonon interfacial resistance into boundary and interfacial-region resistances. These show that the interfacial atomic restructuring as well as the cross-boundary interactions reduce the phonon boundary resistance by providing additional transport channels altering their phonon density of states and cause extra interfacial-region resistances due to additional phonon scattering.

DOI: 10.1103/PhysRevB.82.081302                                              PACS numbers: 68.35.Ja, 02.70.Ns, 63.22.m

Heat transfer across interfaces is central in design of micro/nanodevices, and the acoustic mismatch model (AMM) (for specular surfaces) and diffuse mismatch model (DMM) have been used to estimate the phonon boundary resistance.1 The DMM assumes all phonons are diffusely scattered at the boundary and then transported in proportion to the phonon density of states Dp and the phonon speed up of each side, so the transmission probability from material 1 to 2, τ12(ω), is1


where j is the phonon mode and ω is the phonon frequency. While the AMM is limited to low temperatures (~ 30 K) and smooth surfaces, the DMM provides relatively good agreements with experiments at room or higher temperatures.2 Since the DMM is based on the Debye approximations and bulk properties of contacting materials, it only applies to interfaces with high and similar Debye temperatures TD. Otherwise, the measured bulk phonon dispersion and Dp,3,4 the joint phonon properties for the interfacial region,5 or the virtual crystal representing the disordered interfacial regionhave been used in the DMM. To verify the theoretical models, molecular-dynamics (MD) simulations with no restrictions on the boundary phonon scattering have been used.7,8 However, more realistic interfaces and theoretical models for them have not been addressed, and more accurate predictions are expected with guidance from ab initio calculations, e.g., advances in computing and methods have enabled prediction of the interfacial structures9 and phonon properties10 using the density-functional theory.

In this Rapid Communication, we address the interfacial-region phonon properties using nonequilibrium molecular dynamics (NEMD) and employ cross-boundary reaching force fields (potentials) from ab initio simulations. With interest in using solid-liquid phase change for thermal management of high-power transistors, we consider thermal transport across Si/In (phase change material) interface, which is challenging to treat with the DMM due to the strong contrast in material properties.1 NEMD allows only for the phonon transport with interphonon scattering, however electrons contribute to thermal energy transport, especially in metals. For our semiconductor/metal interface, we address electron (metal) electron (semiconductor) and electron (metal) phonon (semiconductor) interactions for boundary transport involving electrons. The electronic thermal resistance is found from the Wiedemann-Franz law,11 ARe,k=ARe,e /(NL,oT), where ARe,k is the electronic thermal resistance, ARe,e is the electrical resistance, NL,o is the Lorenz constant, and T is the temperature. Unless the semiconductor is very heavily doped, thermal transport by the electron-electron interaction is minor compared with interphonon transport due to the low density and tunneling possibility of electrons to pass through the carrier depletion region from contact.12 [In the Si/In interface with doping concentration 1020 cm−3 for n-type Si, ARe,k= 100– 1000 nK/(W/m2) depending on T, which is ten to 100 times larger than the calculated phonon resistance.] The transport by electron-phonon coupling is negligible at low temperatures but is more significant in the interfacial transport at the high temperatures.13 However, recent experiments14 and theoretical calculations15 show that this coupling has a small contribution to the interfacial transport, and this allows us to focus only on the phonon transport and use NEMD.

For the Si/In layers, the Si-Si and In-In interatomic potentials are the Stillinger-Weber potential for the Si interactions16 and the pairwise Morse potential for In interactions.  We developed the Morse potential from ab initio simulations with Vienna ab initio simulation package (VASP) with the parameters set as φo, In-In= 0.1336 eV, aIn-In= 0.961 Å−1, ro,In-In= 3.41 Å, and rc,In-In(cut-off radius)= 5.235 Å in the Morse potential model {φMorse(rij)=φo[(1−exp−a(rij−ro)])2− 1]} while fitting the experimental thermal properties specific-heat capacity cp= 25.8 J/mole K [26.74 J/mole K at 300 K (Ref. 17)], melting point Tm= 430 K [429.8 K (Ref. 17)], etc.). Other In potentials such as the modified embedded atom model (with stronger bond)18 were considered, however the developed Morse potential predicts the thermal properties most accurately. For the Si-In interactions which affect the interfacial structure, the Morse potential was developed from the force field between Si slab and In atom by ab initio simulations (with parameter set as φo,Si-In=0.3990 eV, aSi-In= 0.8965 Å−1ro,Si-In= 3.402 Å, and rc,Si-In= 5.0 Å), instead of using the combinative rules between Si and In potential models (because  of different forms of Si and In potentials).

NEMD simulations employing these potentials are performed  from 200 to 700 K in the NVT ensemble (constant number of particles, volume, and temperature). The temperature of the system is controlled by the Langevin thermostat, which stochastically adds or removes kinetic energy of the thermostat atoms.19 Periodic boundary conditions are imposed in the x and y directions, and Si layers are placed on both sides of the In layers in the z direction the Si end atoms are fixed and the next-to-end atoms are used as the thermostat atoms.) We use 10 000 In atoms and 5184 Si atoms (2592 each side) are simulated in a rectangular parallelepiped simulation cell of volume 4.84.816 nm3 (Si: 2 nm/In: 12 nm/Si: 2 nm). The interfaces are initially connected as diamond cubic (100) Si surface and face-centered cubic (100) In surface. The initial simulation cell size and density are selected to minimize the lattice mismatch, based on the MD equilibrium density in NpT ensemble (constant number of particles, pressure, and temperature). To examine the interfacial thermal transport, we prescribe a temperature difference across the Si/In/Si system and calculate the steady-state temperature distribution and heat flux. The temperature distribution is from the ensemble-average kinetic energy within divided position bins and the heat flux qMD is the net kinetic, potential and pressure energy flow, i.e.,20


where V is the volume, mi is the mass of particle i, ui is the velocity vector, ui,z is the z component of ui, φi is the potential  energy, zij is the z component of the interatomic separation, and Fij is the interaction force between i and j particles.



where 〈 〉 is an ensemble average and τ is taken as 15 ps which is based on the Δω resolution for Dp taken as 2.094 x 1011 rad/s(=π/τ).

FIG 1. (Color online) (a) Interfacial Si and In atomic structures without interactions (initial equilibrium condition) and (b) with in- teractions and restructuring obtained from the MD simulations at 300 K. The bulk resistance (ARp,Si and ARp,In) and the total interfacial resistance( ARp,Si/In) [decomposed into the boundary (ARp,Si/In) and two interfacial-region resistance ( ARp,Si and ARp,In) are shown through the thermal circuit diagram.

We obtain the atomic structure, temperature distribution and heat flux from the NEMD simulations over 2 ns. Si remains solid, while In, with the Morse potential model, melts at 430 K. Even below this temperature, due to force fields from the adjacent-layer atoms, the Si and In atoms adjacent to the boundary are restructured to satisfy the energy minimization, as presented in Fig. 1b, compared to the initial equilibrium nonrestructured lattice structure in Fig. 1a.

From the thermal analysis of Si/In system, we calculate the phonon conductivity of Si and In in phase i (s for solid, and l for liquid phase) [kp,Si or In,i=−qMD/(dT/dz)] and the thermal resistance ARi=ΔTi /qMD) from qMD and the temperature  distribution, as shown in Fig. 2. kp,In,s decreases from 1.7 to 0.72 W/m K, as the temperature increases from 200 to 320 K. These are in good agreement with the Slack relation22,23 using Dp and density from NEMD but as expected  much smaller than the total thermal conductivity including  electronic contribution, 81.6 W/m K at 300 K.17 Upon melting, the thermal conductivity of In drops and slightly decreases with temperature (0.24 W/m K at 450 K, 0.21 W/m K at 700 K), in good agreement with the Kittel model24 for disordered structures. With the phonon mean-free path λp equal to the average interatomic spacing la (from atomic number density n) we use the kinetic theory result21 kp=ncvup,A λp/3, where up,A is the average phonon propagation speed from elastic constants, and cv is the specific-heat capacity calculated from the MD results for Dp and ƒºp. The NEMD kp,Si obtained using Si atoms only, and when considering the finite-size effect,25 is 139.1 W/m K at 300 K (as compared to experiment 149 W/m K in Ref. 17).

FIG 2. (Color online) Temperature distribution with heat flow- ing along direction z ( cold thermostat at 100 K and hot at 300 K).


FIG. 3 ( Color online) MD phonon density of density of state of Si and In, in the bulk and the interfacial regions (within 0.5c,Si-In from the nominal interface). Due to the interpenetration of force fields and atomic restructuring Dp is distorted in the interfacial region. The blue shaded region represents shared phonon states between Si and In in the interfacial region. The interfacial Dps (solid lines) have more shared phonon states participating in the boundary transport, compare to the bulk of Dps (dashed lines)

Despite exclusion of electronic thermal transport, most of the temperature drop occurs near the boundary, as shown in Fig. 2, thus the bottleneck in thermal transport is the interfacial thermal resistance. Adjacent to the interfaces temperature  distributions are nonlinear because the atomic restructuring  and boundary interaction lead to extra phonon scattering and lower local thermal conductivities. The restructured regions (with nonlinear temperature distributions) extend only on the order of the cut-off radii of the force fields on each side. The extra scattering in these regions is represented with the interfacial-region resistance, which is traditionally included in the thermal boundary resistance. So, the total interfacial resistance ARp,Si/In is sum of the boundary resistance ARp,b,Si/In and two interfacial-region resistances ARp,r,Si and ARp,r,In (the corresponding temperature jumps for the three resistances are shown in Fig. 2), ARp,Si/In=ARp,b,Si/In+ARp,r,Si+ARp,r,In.

In these interfacial regions, Dp is modified due to the atomic restructuring and interactions between Si and In atoms,  as presented in Fig. 3. The interfacial-region Dp,Si and Dp,In are calculated using Eq. (3) within the region bounded by the nominal interface and 0.5rc,Si-In, while the bulk values are calculated with 200 atoms in the central region (in z direction), as marked in Fig. 2. Dp,Si shifts to the lower energies  and Dp,In moves to higher energies, and the changes in Dp,In is more pronounced than Dp,Si.

Fig 4. (Color online) (a) Comparison of the variations in the NEMD thermal interfacial-region resistances (in In and Si) as a function of temperature with the approximate relation δr,i/kp,r,i (b) Comparison of the variations in the NEMD Si/In thermal boundary resistance as a function of temperature with the DMM theory. The liquid NEMD simulation results have larger scatter due to the local temperature fluctuations. The scatter increases with increase in temperature (larger kinetic energy) and with samaller heat flux (lower thermal conductivity)


A sample temperature distribution in the restructured regions, shown in Fig. 2, suggests lower thermal conductivity kp,r or reduced λp, due to enhanced scattering. Then the kp,r decreases as the boundary is approached because the λp decreases to la near the diffusely scattering boundary. Moreover, , asymmetric force field near the boundary reduces up,A (depending on the interaction potential), it additionally decreases the phonon conductivity. Using n and cv from the MD density and Dp in the interfacial region with λp=la=(2½ /n)¹⁄³ and MD interfacial region phonon conductivity kp,r=(ncvup,Ala)/3 (kp,r,In≈0.1 and kp,r,Si0.3 W/m K), the interfacial-region up,A is estimated to be 0.4–0.5 for the bulk up,A. We define the restructuring length as δr,i≡kp,r,iARp,r,i, and the dimensionless δr,i(=δr,i/ro,Si-In) from NEMD results are δr,In,s≅1.6, δr,In,I≅0.6, and δr,Si≅1.1 as listed in Fig. 4a. [These solid-phase restructuring lengths are in agreement with the locations where the average (over three principal reciprocal vectors) static structure factor26 is 0.5.] These indicate that the interfacial restructuring for solid/solid extends further into the softer material, and the solid In restructuring extends over the cut-off radius rc,Si/In of the effective force (potential) of interaction with Si. Since the interfacial-region resistance is the decrease in thermal conductivity in the interfacial region, the softer In with extended restructured region has the larger interfacial-region resistance [~5 nK/(W/m2)] compared to the harder Si [~1.3 nK/(W/m2)]. The liquid-phase In interfacial-region resistance is smaller than the solid phase because liquid already has a disordered structure and lower thermal conductivity, as shown in Fig. 4a.

For ARp,b, we assume that all phonons with energy hω are traveling to the boundary region with average speed up,A and are then completely scattered as assumed in the DMM theory because here the diffuse interfacial scattering is dominant at temperatures over 200 K and in the presence of the noncrystalline phases at the boundary. Using this phonon speed and Dp for both interfacial regions with the DMM T1→2(ω) from Eq. (1), the phonon boundary resistance is1


where ωm is the cut-off frequency of the softer material. From Eq. (1), only phonons which can exist on the both sides contribute to the phonon transport with nonzero T1→2(ω). Figure 3 shows in their interfacial regions In and Si have larger range of shared phonon states, about four times compared to the bulk (and these participate in the boundary transport). The retarded up,A in the interfacial region also reduces the phonon energy flux at the boundary. In Fig. 4b, the DMM phonon boundary resistances ARp,b,DMM[Dp,r(T)] calculated with the interfacial properties [11− 13 nK/(W/m2)] are in good agreement with the NEMD results (ARp,b,MD=ΔTp,b /qMD), whereas using the Debye Dp(ARp,b,DMM[DDebye]) and the bulk phonon properties ARp,b,DMM[Dp,bulk(T)]) overestimates ARp,b using the bulk MD properties of Si and In gives 21.55 while the Debye model give 82.43 nK/(W/m2)]. At higher temperatures, the controlling, softer material In has more high-energy phonons participating in the transport, and ARp,b decreases. For the liquid In, the reverse occurs because Dp will have more low-energy (longer wavelength) modes due to increased anharmonicity.

We analyzed the interfacial thermal transport by treating the interfacial atomic restructuring and its phonon transport with NEMD, for Si/In interface. We decomposed the interfacial phonon resistance into boundary and interfacial-region components. The atomic restructuring and the boundary interactions were found to reduce the boundary resistance by providing additional phonon transport channels (modes) while adding extra scattering resistance in the interfacial regions. The DMM using the interfacial-region phonon properties predicts the boundary resistance in good agreement with NEMD, even with low TD materials, and this is also the case when one material is in the melt phase. For the interfacial-region resistances, we provide a simple model based on minimum λp. In the NEMD thermal contact analysis, without the interfacial-region resistances (due to the noncrystalline region which has lower local thermal conductivity than crystalline phase), the bulk behavior away from the interface is not recovered. Adding the boundary resistance from the DMM and the interfacial-region resistances, this is the total interfacial resistance found in experiments. While for the Si/In example considered the total interfacial resistance is less than that predicted by the DMM based on bulk properties, this cannot be generalized for less contrasting material pairs. More accurate predictions of the interfacial region atomic structures and phonon properties, such as the ab initio MD (with the Car-Perrinello method combining the classical dynamics and electronic structures27,28), will improve on these treatments.

We are thankful for financial support by the National Reconnaissance Office, Applied Technology Department.



1E. T. Swartz and R. O. Pohl, Rev. Mod. Phys. 61, 605 (1989).

2R. M. Costescu, M. A. Wall, and D. G. Cahill, Phys. Rev. B 67, 054302 (2003).

3P. E. Phelan, ASME J. Heat Transfer 120, 37 (1998).

4P. Reddy, K. Castelino, and A. Majumdar, Appl. Phys. Lett. 87, 211908 (2005).

5P. E. Hopkins and P. M. Norris, Nanoscale Microscale Thermophys. Eng. 11, 247 (2007).

6T. Beechem, S. Graham, P. E. Hopkins, and P. M. Norris, Appl. Phys. Lett. 90, 054104 (2007).

7E. S. Landry and A. J. H. McGaughey, Phys. Rev. B 80, 165304 (2009).

8R. J. Stevens, L. V. Zhigilei, and P. M. Norris, Int. J. Heat Mass Transfer 50, 3977 (2007).

9X.-G. Wang and J. R. Smith, Phys. Rev. Lett. 95, 156102 (2005).

10J. Fritsch and U. Schröder, Phys. Rep. 309, 209 (1999).

11L. W. da Silva and M. Kaviany, Int. J. Heat Mass Transfer 47, 2417 (2004).

12A. Kikuchi, Phys. Status Solidi A 175, 623 (1999).

13A. V. Sergeev, Phys. Rev. B 58, R10199 (1998).

14H. K. Lyeo and D. G. Cahill, Phys. Rev. B 73, 144301 (2006).

15G. D. Mahan, Phys. Rev. B 79, 075408 (2009).

16F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985).

17W. M. Haynes, CRC Handbook of Chemistry and Physics, 91st ed. (CRC Press, Boca Raton, 2009).

18E. C. Do, Y.-H. Shin, and B.-J. Lee, CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 32, 82 (2008).

19C. S. Wang, J. S. Chen, J. Shiomi, and S. Maruyama, Int. J. Therm. Sci. 46, 1203 (2007).

20G. S. Hwang and M. Kaviany, J. Appl. Phys. 106, 024317 (2009).

21M. Kaviany, Heat Transfer Physics (Cambridge, New York, 2008).

22G. A. Slack, in Solid State Physics, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic, New York, 1979), Vol. 34, pp. 1–73.

23B. L. Huang and M. Kaviany, J. Appl. Phys. 100, 123507 (2006).

24C. Kittel, Phys. Rev. 75, 972 (1949).

25P. K. Schelling, S. R. Phillpot, and P. Keblinski, Phys. Rev. B 65, 144306 (2002).

26T. G. Desai, J. Chem. Phys. 127, 154707 (2007).

27R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).

28T. Luo and J. R. Lloyd, ASME J. Heat Transfer 130, 122403 (2008).

Have a Question or Project to Discuss?