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 *D _{p} *and the phonon speed

*u*of each side, so the transmission probability from material 1 to 2, τ

_{p}_{1}→

_{2}(ω), is

^{1}

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 *T _{D}*. Otherwise, the measured bulk phonon dispersion and

*D*,

_{p}^{3,4}the joint phonon properties for the interfacial region,

^{5}or the virtual crystal representing the disordered interfacial region

^{6 }have 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 structures

^{9}and phonon properties

^{10}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} *AR _{e,k}=AR_{e,e} /(N_{L,o}T)*, where

*AR*is the electronic thermal resistance,

_{e,k}*AR*is the electrical resistance,

_{e,e}*N*is the Lorenz constant, and

_{L,o}*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 10

^{20}cm

^{−3}for

*n*-type Si,

*AR*= 100– 1000 nK/(W/m

_{e,k}^{2}) 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 experiments

^{14}and theoretical calculations

^{15}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 interactions^{16} 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, *a*In-In= 0.961 Å^{−1}, ro,In-In= 3.41 Å, and *r*_{c,In-In}(cut-off radius)= 5.235 Å in the Morse potential model {φMorse(rij)=φ_{o}[(1−exp−a(r_{ij}−r_{o})])^{2}− 1]} while fitting the experimental thermal properties specific-heat capacity *c _{p}*= 25.8 J/mole K [26.74 J/mole K at 300 K (Ref. 17)], melting point

*T*= 430 K [429.8 K (Ref. 17)], etc.). Other In potentials such as the modified embedded atom model (with stronger bond)

_{m}^{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, a

_{Si-In}= 0.8965 Å

^{−1},

*r*

_{o,Si-In}= 3.402 Å, and

*r*

_{c,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 nm^{3} (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 q_{MD} is the net kinetic, potential and pressure energy flow, i.e.,^{20}

where V is the volume, *m _{i}* is the mass of particle

*i*,

*is the velocity vector,*

**u**_{i}*u*is the

_{i,z}*z*component of

*,*

**u**_{i}*φi*is the potential energy,

*z*is the

_{ij}*z*component of the interatomic separation, and

*is the interaction force between*

**F**_{ij}*i*and

*j*particles.

where 〈 〉 is an ensemble average and *τ* is taken as 15 ps which is based on the Δω resolution for *D _{p}* taken as 2.094 x 10

^{11}rad/s(=π/

*τ*).

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) [k_{p,Si} _{or In,i}=−q_{MD}/(d*T*/dz)] and the thermal resistance AR_{i}=ΔT_{i} /q_{MD)} from q_{MD} and the temperature distribution, as shown in Fig. 2.* k _{p,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 relation

^{22,23}using

*D*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.

_{p}^{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 model

^{24}for disordered structures. With the phonon mean-free path λ

*equal to the average interatomic spacing*

_{p}*l*(from atomic number density

_{a}*n)*we use the kinetic theory result

^{21 }

*k*where

_{p}=nc_{v}u_{p},_{A}λ_{p}/3,*u*is the average phonon propagation speed from elastic constants, and

_{p,A}*c*is the specific-heat capacity calculated from the MD results for

_{v}*D*and ƒº

_{p}_{p. }The NEMD

*k*obtained using Si atoms only, and when considering the finite-size effect,

_{p,Si}^{25}is 139.1 W/m K at 300 K (as compared to experiment 149 W/m K in Ref. 17).

** **

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 *AR _{p,Si/In}* is sum of the boundary resistance

*AR*and two interfacial-region resistances

_{p,b,Si/In}*AR*and

_{p,r,Si}*AR*(the corresponding temperature jumps for the three resistances are shown in Fig. 2),

_{p,r,In}*AR*=

_{p,Si/In}*AR*+

_{p,b,Si/In}*AR*+

_{p,r,Si}*AR*

_{p,r,In.}In these interfacial regions, *D _{p}* is modified due to the atomic restructuring and interactions between Si and In atoms, as presented in Fig. 3. The interfacial-region

*D*and

_{p,Si}*D*,In are calculated using Eq. (3) within the region bounded by the nominal interface and 0.5

_{p}*r*, while the bulk values are calculated with 200 atoms in the central region (in

_{c,Si-In}*z*direction), as marked in Fig. 2.

*D*shifts to the lower energies and

_{p,Si}*D*moves to higher energies, and the changes in

_{p,In}*D*is more pronounced than

_{p,In}*D*.

_{p,Si}** **

A sample temperature distribution in the restructured regions, shown in Fig. 2, suggests lower thermal conductivity *k _{p,r}* or reduced λ

*p*, due to enhanced scattering. Then the

*k*decreases as the boundary is approached because the λ

_{p,r}_{p}decreases to

*l*near the diffusely scattering boundary. Moreover, , asymmetric force field near the boundary reduces

_{a}*u*(depending on the interaction potential), it additionally decreases the phonon conductivity. Using

_{p,A }*n*and

*c*from the MD density and

_{v}*D*in the interfacial region with λ

_{p}*p*=

*l*=(2

_{a}^{½}/

*n)¹⁄³*and MD interfacial region phonon conductivity

*k*=

_{p,r}*(nc*)/

_{v}u_{p,A}l_{a}*3 (k*0.1 and

_{p,r,In}≈*k*0.3 W/m K), the interfacial-region

_{p,r,Si}≈*u*is estimated to be 0.4–0.5 for the bulk

_{p,A}*u*. We define the restructuring length as

_{p,A}*δ*, and the dimensionless

_{r,i}≡k_{p,r,i}AR_{p,r,i}*δ*(=

_{r,i}*δ*/

_{r,i}*r*) from NEMD results are

_{o,Si-In}*δ*≅1.6,

_{r,In,s}*δ*≅0.6, and

_{r,In,I}*δ*≅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 factor

_{r,Si}_{26}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

*r*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/m

_{c,Si/In}^{2})] compared to the harder Si [~1.3 nK/(W/m

^{2})]. 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 *AR _{p,b}*, we assume that all phonons with energy

*hω*are traveling to the boundary region with average speed

*u*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

_{p,A}*D*for both interfacial regions with the DMM

_{p}_{T1→2}(ω) from Eq. (1), the phonon boundary resistance is

^{1}

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

*u*in the interfacial region also reduces the phonon energy flux at the boundary. In Fig. 4b, the DMM phonon boundary resistances

_{p,A}*AR*

_{p,b,}

_{DMM}_{[}

_{Dp,r(T}_{)]}calculated with the interfacial properties [11− 13 nK/(W/m

^{2})] are in good agreement with the NEMD results (

*AR*=

_{p,b,MD}*ΔT*/

_{p,b}*q*), whereas using the Debye

_{MD}*D*(

_{p}*AR*

_{p,b,DMM}_{[DDebye]}) and the bulk phonon properties

*AR*[

_{p,b,DMM}*D*

_{p,bulk}_{(T)}]) overestimates

*AR*using the bulk MD properties of Si and In gives 21.55 while the Debye model give 82.43 nK/(W/m

_{p,b}^{2})]. At higher temperatures, the controlling, softer material In has more high-energy phonons participating in the transport, and

*AR*decreases. For the liquid In, the reverse occurs because

_{p,b}*D*will have more low-energy (longer wavelength) modes due to increased anharmonicity.

_{p}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 *T _{D}* 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 structures

^{27,28)}, will improve on these treatments.

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

*kaviany@umich.edu

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

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

^{3}P. E. Phelan, ASME J. Heat Transfer **120**, 37 (1998).

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

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

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

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

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

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

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

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

^{12}A. Kikuchi, Phys. Status Solidi A **175**, 623 (1999).

^{13}A. V. Sergeev, Phys. Rev. B **58**, R10199 (1998).

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

^{15}G. D. Mahan, Phys. Rev. B **79**, 075408 (2009).

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

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

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

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

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

^{21}M. Kaviany, *Heat Transfer Physics* (Cambridge, New York, 2008).

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

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

^{24}C. Kittel, Phys. Rev. **75**, 972 (1949).

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

^{26}T. G. Desai, J. Chem. Phys. **127**, 154707 (2007).

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

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