Spectral Analysis of the Heat Flow Across Crystalline and Amorphous

Apr 27, 2017 - Bladimir Ramos-Alvarado and Satish Kumar. The George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, ...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Spectral Analysis of the Heat Flow Across Crystalline and Amorphous Si−Water Interfaces Bladimir Ramos-Alvarado* and Satish Kumar The George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, United States ABSTRACT: Nonequilibrium classical molecular dynamics simulations were employed to investigate thermal transport across crystalline and amorphous silicon (a-Si) surfaces in contact with water for different interfacial bonding strengths. A spectral analysis of heat transfer across the different interfaces revealed the characteristics of the phonon modes contributing to thermal transport. Low-frequency modes contributed the most in hydrophobic interfaces, while a shift toward contribution from higher frequency modes was found for hydrophilic surfaces. The shift to higher frequency modes was not significant for a-Si and crystalline Si(111) interfaces. In-plane phonon modes significantly contributed to heat transfer in Si(100), less significantly in a-Si, and had a minimum contribution in Si(111) hydrophilic interfaces. While the wettability and solid−liquid bonding strength failed in explaining these observations, the interfacial liquid density depletion helped to understand the differences between Si(100) and a-Si interfaces with respect to the Si(111) interface. The interface liquid structure observed in the Si(100) but not in the a-Si system served as an explanation for the dominant contribution of in-plane modes in Si(100). These observations posed the density depletion and liquid structure at solid−liquid interfaces as useful parameters for explaining the underlying mechanisms of phonon transport at solid−liquid interfaces.

1. INTRODUCTION Thermal transport across solid−solid and solid−liquid interfaces has caught the attention of investigators for many years, mainly due to the need to enhance the heat dissipation capability in microelectronics and nanoscale devices.1−3 The reason for poor thermal energy dissipation in nanoscale devices is the large number of interfaces created during their fabrication. Thermal transport across solid−liquid interfaces has recently gained more attention in the field of biotechnology for applications such as cancer treatment and drug delivery, friction reduction in microscale gaps, and cooling using nanofluids. Kapitza4 reported the anomaly of heat transfer across interfaces. This phenomenon had been proposed and studied before, but ref 4 is one of the first scientific reports on this topic.5 Nowadays, it is well-known that a temperature discontinuity is observed when heat is transferred between different materials across an interface due to a mismatch in the vibrational properties. Such a mismatch gives origin to a thermal boundary resistance (TBR) or otherwise characterized as a thermal boundary conductance (TBC) defined as q = GΔTjump, where G is the TBC and q and ΔTjump are the heat flux and temperature discontinuity observed at the interface, respectively. Numerical simulations using classical molecular dynamics (MD), particularly nonequilibrium MD (NEMD), have been widely used to understand thermal transport across solid− liquid interfaces. Early investigations focused on characterizing the TBR as a function of the nonbonded solid−liquid interaction strength.6−9 The main takeaway from these © 2017 American Chemical Society

investigations is that the TBR decreases as the wettability of the solid surfaces goes from hydrophobic to hydrophilic; however, assuming that wettability only depends on the interfacial bonding strength gives a partial understanding of the wettability phenomenon.10,11 More recent investigations have focused on the correlation of the contact angle obtained from equilibrium MD (EMD) simulations and the TBC obtained from NEMD simulations of heat transfer. Merabia et al.12 verified that hydrophilic gold has a lower TBC than hydrophobic gold. Shenogina et al.13 found that the TBC scales linearly with the work of adhesion as G ∼ 1 + cos(θ), where θ is the macroscopic contact angle derived from nonsize affected EMD simulations of wettability. Acharya et al.14 confirmed the previous scaling law, showing that G ∼ cos(θ) when the wettability is modified using surface chemical functionalization. More recently, Ramos-Alvarado et al.15 verified the TBCwettability scaling laws and demonstrated that these are not universal. This was accomplished by finding a scaling law that reconciled the anisotropy (crystalline directional dependence) of the wettability of different silicon planes but not the anisotropy of the TBC for different silicon planes. Liquid layering and depletion are observed at solid−liquid interfaces depending on the solid−liquid affinity and the planar concentration of solid atoms at the interface. Xue et al.16 noted that the interfacial liquid layering did not influence the thermal conductivity of the liquid but indicated the possibility of having Received: February 21, 2017 Revised: April 16, 2017 Published: April 27, 2017 11380

DOI: 10.1021/acs.jpcc.7b01689 J. Phys. Chem. C 2017, 121, 11380−11389

Article

The Journal of Physical Chemistry C an effect on the TBC. Tori et al.,17 after conducting a directional analysis of interfacial thermal transport, found that the in-plane modes had a strong contribution to the total heat transfer at solid−liquid interfaces and formulated the possibility of a TBC−liquid layering dependence. Similarly, Murad and Puri18,19 observed the reduction of the TBR at interfaces where the interfacial liquid layering was more pronounced, leading to more fluid adsorption. Ramos-Alvarado et al.15 added to this discussion and reported for the first time a quantification of the relationship between the TBC and the liquid depletion at the interface, where the density depletion was reported as a more general parameter than the contact angle to explain thermal transport at solid−liquid interfaces. In order to enhance the thermal management efficiency of devices where heat transfer occurs at interfaces, an accurate analysis of the phonon modes contributing to thermal transport is crucial, especially for the rather ignored solid−liquid interfaces. Carlborg et al.20 numerically investigated thermal transport in single-wall carbon nanotubes (SWCNTs) surrounded by solid and liquid argon matrices. They found that low-frequency phonons in SWCNT couple with those of the adsorbed atoms onto the SWCNT. This was observed by analyzing the phonon density of states (DOS). Goicochea et al.21 found that low-frequency phonons (0−5 THz) dominated heat transfer across quartz−water interfaces. They also found an attenuation and down-shifting of the frequency peak in the given frequency range by analyzing the DOS. Giri and Hopkins22 observed that the highest density peak of the DOS shifted toward lower frequencies, when compared with the DOS in the bulk, under poorly wetting conditions, indicating a resemblance with a free surface. Saaskilathi et al.23 recently developed a model to perform a spectral analysis of heat transfer without relying on inferences made from the DOS as done in previous investigations. Saaskilathi et al.24 used the model reported in ref 23 and observed a maximum contribution to the interfacial heat transfer from low-frequency modes. The peak contribution for heat transfer in the bulk was shifted toward higher frequencies, a similar observation as that derived from the DOS at a free surface indicated in ref 22. Giri et al.25 recently implemented the model by Saaskilathi et al.23 to solid−solid, solid−liquid, and solid−gas interfaces. It was found that by increasing the solid−liquid interaction strength higher frequency modes contribute to the interfacial heat transfer. For their particular system, ∼80% of the total heat flow was carried by low-frequency phonons. Spectral analyses of interfacial heat transfer have mostly been focused on generic and computationally affordable LennardJones (LJ) systems. In this investigation, a spectral analysis of heat transfer was conducted on a realistic system composed of water interacting with crystalline and amorphous silicon. A reliable water model was used, where van der Waals and electrostatic interactions are considered. The silicon interactions were modeled with a multibody interaction potential. Different wettability conditions were investigated, and the model reported in ref 23 was used to map the spectral contributions to the interfacial heat transfer. The results indicated that some of the analyses performed on generic systems can be extrapolated to these realistic interfaces; however, the qualitative contributions to heat transfer are different. The concept of wettability and the solid−liquid interaction strength parameter were used in an effort to explain the different TBC calculations observed at the different solid− liquid interfaces. These provided a partial description by merely

indicating that more hydrophilic surfaces (stronger solid−liquid bonding strength) are more conductive than hydrophobic surfaces (weaker solid−liquid bonding strength). Alternatively, the liquid depletion at the interface provided a better explanation for why some interfaces having the same wettability have different thermal conductance. Likewise, the wettability of surfaces could not be used to justify the drastic differences between the phonon modes contributing to interfacial heat transfer. The hydrophilic Si(111) interface had a minimum contribution of in-plane phonon modes, while the a-Si and Si(100) counterparts presented a more important contribution of these modes. The liquid density depletion at the interface helped to get a better understanding of these underlying mechanisms. Lastly, the combination of density depletion and interface water structure was used to explain the conditions necessary for in-plane phonon modes to dominate interfacial heat transfer at hydrophilic Si(100) interfaces.

2. METHODOLOGY 2.1. Generation of the Amorphous Silicon Structure. Amorphous silicon (a-Si) can be generated by a melting− quenching process in MD simulations26−29 or by randomization methods such as the Wooten−Winer−Weaire bond switching algorithm.30,31 Due to its simplicity of implementation, the melting−quenching process was applied herein to generate a-Si. First, a supercell of crystalline Si made of 8000 atoms was melted at a constant temperature and pressure of 4000 K and 1 bar (NPT), respectively, for 50 ps. Two quenching steps followed: a quick one to lower the temperature from 4000 to 1500 K in 50 ps and a slow one to bring the system from 1500 to 300 K in 10 ns. The pressure was kept at 1 bar during both quenching processes. After the system reached a constant density, the simulation was switched to a canonical integration (constant volume and temperature or NVT) where the temperature was kept at 300 K and the energy variation of the system was tracked for 20 ns. Lastly, to verify the equilibrium of the structure, the integration was switched to a purely Newtonian integration (constant energy), and the temperature of the system was tracked for 1 ns. The Tersoff interaction potential32 was used for the silicon, and the simulations were run using LAMMPS33 with a time step of 1 fs. A slab of the supercell was taken and duplicated; afterward, both slabs were used to enclose a water block. The process outlined above closely follows that reported in ref 29; however, the implementation of the methodology reported in ref 29 led to an unstable system when the a-Si was in contact with water and running under a microcanonical ensemble (NVE), a fundamental step in the NEMD simulation of heat transfer. In order to solve this problem, a second slower quenching step was introduced, in addition to tracking the energy and temperature of the system during the last NVT and NVE equilibration steps. 2.2. NEMD Simulation of Heat Transfer. Heat transfer was investigated at the interface between crystalline and amorphous silicon with water. The SPC/E34 water model was implemented due to its relatively accurate prediction of thermal conductivity when compared with other water models,35 for being computationally affordable, and because of its widespread application, which allows for an easier comparison with previous investigations. The rigidity of the SPC/E model was enforced using the SHAKE36 algorithm, and the long-range electrostatic interactions were handled using the PPPM37 method with an accuracy of 1 × 10−6. The Si−Si interactions in 11381

DOI: 10.1021/acs.jpcc.7b01689 J. Phys. Chem. C 2017, 121, 11380−11389

Article

The Journal of Physical Chemistry C Si and a-Si were handled using the multibody Tersoff32 potential. The solid−liquid interactions were modeled as van der Waals forces through a 12-6 Lennard-Jones potential where only the interacting pair Si−O was considered. The parameters were σSi−O = 2.63 Å, and εSi−O was varied in order to emulate different wettability conditions according to ref 15. Figure 1 depicts renderings of the systems used in this investigation, where two silicon−water interfaces are created by

per bin in which the water slab was divided. A total of 50 slabs were used to discretize the 6 nm water slab. For water T = KE/ 3kB since the rigid SPCE water model has only 6 degrees of freedom. The temperature profiles were fitted with linear functions, and ΔTjump was calculated from the extrapolation of the temperature profiles at the interfaces. Simulation runs using different heat fluxes (q) were carried out to verify the linear response of these nanoscale systems as q = GΔTjump, at the same time that G was calculated from the slope of the linear fit to the data. This process was implemented for the different solid−liquid interfaces under a wide range of wetting conditions.

3. SPECTRAL ANALYSIS OF INTERFACIAL HEAT TRANSFER The first step in the description of the spectral analysis of interfacial heat transfer was the calculation of the phonon density of states (DOS) DOS(ω) ∼ 1⟨v(t )·v(0)⟩

(1)

where the DOS is proportional to the Fourier transform (1 ) of the ensemble average of the velocity autocorrelation function (⟨v(t)·v(0)⟩).42 Equation 1 can be simplified by acknowledging that 1 cancels out with its inverse operator involved in the calculation of the autocorrelation function; therefore, the only operation left is obtaining the magnitude of the 1 of the velocity vectors. The DOS was calculated for solid atomic monolayers in the bulk and at the interfaces in contact with water. The velocity components of the atoms were sampled during the NVE equilibration step every 25 fs, and 250 time steps were taken for each Fourier transform and averaged for the total equilibration run. Second, the spectral heat flux, q(ω), reported in refs 23 and 24, was evaluated as

Figure 1. Renderings of the systems under investigation: (a) a-Si and water, (b) Si(100) facing water, and (c) Si(111) facing water. Heat transfer occurs in the longitudinal direction of the figures (z-direction).

two silicon slabs confining a water block. Figure 1(a) illustrates the a-Si−water system. For this particular case, the a-Si slabs were obtained from the super cell described in the previous section. Figures 1(b) and (c) represent Si slabs confining water where the planes (100) and (111) are facing the liquid molecules, respectively. In Figure 1 the z-direction unit vector is normal to the solid−liquid interfaces and corresponds to the heat flow direction. A solid slab length of 10 nm and areas of 2.7 × 2.2 nm2 for a-Si and Si(100) and 3.1 × 2.7 nm2 for Si(111) were the dimensions that resulted in nonsize affected calculation of the TBC. Periodic boundary conditions were imposed in the x- and y-directions (in-plane), while the outermost atoms were kept fixed in both confining solid slabs. Heat input\output regions 1.5 nm thick were defined right after the location of the fixed atoms. The size of the water confinement was 6 nm, and the number of water molecules was changed according to the wettability condition in order to eliminate any compressibility effects.38 The LAMMPS33 code was used for the MD simulations and VMD39 for visualization. The simulations followed the steps: (1) energy minimization, (2) equilibration at a constant temperature of 300 K for 1.0 ns using a Nose−Hover40,41 thermostat with a time constant of 0.1 ps, (3) equilibration in the microcanonical ensemble (NVE) for 1.0 ns in order to verify the stability of the system’s temperature without having any thermostat, (4) thermal energy addition/removal in the heat input regions for 5 ns, and (5) data sampling of the kinetic energy and coordinates of the atoms every 10 ps during a production run of 5 ns. The temperature profiles of the solid slabs were obtained by time-averaging the kinetic energy (KE) per atomic plane for Si and per bin for a-Si, where T = KE/ 1.5kB. The temperature profile of the water slab was obtained by means of time-averaging and particle-count averaging the KE

q(ω) =

2 Re ∑ ∑ 1⟨Fij(t ) ·vi(0)⟩ A j∈l i∈s

(2)

where q(ω) is calculated as the real part of the Fourier transform of the force−velocity cross-correlation function ⟨Fij(t)·vi(0)⟩. In eq 2 A is the cross-section area of heat flow; Fij is the total force acting on the solid atom i due to all liquid particles; and vi is the velocity vector of the solid atom i. The data required to calculate q(ω) was recorded during the nonequilibrium production run, i.e., the 5 ns where the KE energy was sampled according to the previous section. A region in the solid slabs near the interfaces was defined within the cutoff of the solid−liquid interaction potential (13 Å). The velocity was recorded for the solid atoms in this region and the coordinates of all the atoms every 25 fs. The coordinate file contained the dynamics of the system, such that it was possible to reset the interaction potentials in order to isolate Fij while getting this information in a per atom basis as required in eq 2. The data were averaged for the total production run.

4. RESULTS AND DISCUSSION 4.1. Amorphous Silicon Properties. Prior to conducting any analyses, it was of the utmost importance to verify the properties of a-Si since this material was generated by following an empirical methodology. The first property that was verified was the density of the supercell at 300 K after the slow NPT quenching. The density of a-Si is saturated at ∼2.3 g/cm3, 11382

DOI: 10.1021/acs.jpcc.7b01689 J. Phys. Chem. C 2017, 121, 11380−11389

Article

The Journal of Physical Chemistry C

response of the system, meaning that the q = GΔTjump relationship was not followed. From the data shown in Figure 3, it could be inferred that both the Si(100) and Si(111) interfaces are equally conductive; however, a sampling of the TBC at larger values of εSi−O demonstrates that the Si(100) interface is more conductive than the Si(111) interface. An explanation for this anisotropic behavior was given in terms of the interfacial liquid density depletion in ref 15. Figure 3 indicates that the a-Si interface is more conductive than the two crystalline Si planes even though the same εSi−O range is considered and similar wettability conditions exist. In Section 4.4, a spectral decomposition of the phonon modes contributing to interfacial heat transfer is conducted in an effort to explain the TBC behavior of the surfaces under investigation. In Section 4.5, the analysis of the interfacial liquid structure is used to provide further understanding of the TBC in different surfaces. It has been demonstrated that εs1 ∼ 1 + cos(θ) and εs1 ∼ 180° − θ,38 where εsl is an arbitrary nonbonded interaction strength between a solid (s) and a liquid (l); likewise it has been numerically and experimentally reported that G ∼ 1 + cos(θ) and G ∼ 180° − θ.13,38,44 Figure 3 illustrates that a linear relationship G ∼ εSiO is obtained for every interface as expected, a nongeneral scaling law, just like the TBC− wettability and εsl−wettability relationships.15 The contact angle obtained from EMD simulations does not only depend on the set of parameters that determine the interfacial bonding strength but also on the characteristics of the solid structure. Additionally, the combination of these variables determines the liquid depletion at the interface, which in turn plays a major role in explaining the interfacial heat transfer.15 In an early investigation,8 two regimes of the TBR as a function of εsl were found, but no characterization of the resulting contact angle was conducted for the given values of εsl. Giri and Hopkins22 and Giri et al.25 reported a linear regime of G for large εsl and nonlinear for small εsl. The distinction between large and small εsl was given by the relationship between εsl and the liquid−liquid interaction parameter. This characterization can be done for simple generic systems as the ones investigated in refs 9, 22, and 25; therefore, a contact angle characterization could be a more appropriate approach for realistic systems as indicated above. 4.3. Phonon Density of States. The first part of the spectral analysis of heat transfer across the Si−water and a-Si− water interfaces was the calculation of the DOS at interfacial and bulk atomic monolayers using eq 1. Figure 4(a) depicts these calculations for the Si(111) and Si(100) interfaces. It can be said that there is no drastic change in the DOS in the range of 6−14 THz for both interfaces in comparison with the bulk. It is observed that the density peak of the lower-frequency modes increases and down-shifts in frequency, when compared with the bulk. The most drastic change of the interfacial DOS is observed at ∼15 THz. The large characteristic peak of DOS of the bulk Si45 is damped in both Si(111) and Si(100) interfaces, the latter one having the most significant change. According to Giri and Hopkins,22 this could be interpreted as an increasing transmission of these phonon modes across the interface, which was verified in Figure 5. Giri and Hopkins22 found that a shift toward lower frequencies in DOS peaks occurred for the phonon modes under poorly wetting conditions (very small εsl), indicating a resemblance with a 2-D noninteracting surface. Additionally, the magnitude of the DOS peak of the lower-frequency modes

consistent with that of Si. Second, the structure factor, a parameter used to characterize the atomic structure of materials, especially crystals, was used to analyze the atomic arrangement of atoms in the a-Si structure. Figure 2 depicts a

Figure 2. Comparison between the numerical and experimental structure factor for a-Si at 300 K. The experimental data were taken from ref 27.

comparison between the structure factor calculated for the a-Si super cell and experimental data reported in ref 27. The good match between the numerical and experimental data observed in Figure 2 suggests a proper generation of the a-Si structure. Lastly, and jumping ahead of the upcoming sections, the thermal conductivity of the a-Si slabs was calculated using the heat flux and the temperature gradient for the vast number of simulations conducted. The thermal conductivity of a-Si was 1.40 ± 0.22 W/m K, a value close to 1.55 ± 0.22 W/m K43 and 1.7 W/m K.28 4.2. Solid−Liquid Bonding Strength. The effect of the solid−liquid bonding strength on the TBC (G) was investigated by means of varying the energy parameter in the 12-6 Lennard-Jones potential εSi−O. The contact angles of the surfaces investigated are 147.4°−73.9° for Si(100), 153.1°− 83.4° for Si(111), and 141.3°−87.4° for a-Si and correspond to the range from small to large εSi−O. Figure 3 depicts the

Figure 3. Thermal boundary conductance dependence on the strength of the solid−liquid nonbonded interaction.

relationship between G and εSi−O, where linear functions were the most appropriate fits for the data obtained for Si in different planes and a-Si. A narrower range of εSi−O than that reported in ref 15 was used in this investigation. The reason was that, as indicated in Figure 3, the a-Si−water interface is highly conductive; therefore, significantly large heat fluxes had to be applied in order to distinguish a noise-free temperature discontinuity at the interface. Large heat fluxes created unrecognizable temperature discontinuities and a nonlinear 11383

DOI: 10.1021/acs.jpcc.7b01689 J. Phys. Chem. C 2017, 121, 11380−11389

Article

The Journal of Physical Chemistry C

Figure 4. Phonon density of states (DOS) in arbitrary units for bulk and interfacial atoms in (a) crystalline and (b) amorphous silicon.

decreased as the bonding strength increased, suggesting an opening of a channel for energy transport from the solid to the liquid due to the migration or redistribution of energy between energy carriers at the interface. The aforementioned shift and damping of the low-frequency peak were not observed for the Si−water interfaces, probably because the values of εsl were calibrated such that a set of macroscopic contact angles could be observed from droplet wettability simulations or due to the different level of complexity of the model used in these simulations as opposed to LJ particles. The utilization of generic LJ systems allows for using any value of εsl, such that a free surface condition can be observed. Additionally, the drastic shift of the DOS, from hydrophobic to hydrophilic surfaces reported in ref 22, was not observed in this investigation. The changes between hydrophobic and hydrophilic conditions were minimum, thus only one curve per plane is illustrated in Figure 4(a). Figure 4(b) depicts the DOS for interfacial and bulk a-Si for two wettability conditions. Similar to the Si case, a transition as clear as that reported in ref 22 for different wetting conditions was not observed for the amorphous interface. The DOS for the a-Si bulk atoms is quite different and smoother than that for crystalline Si, mainly due to the lack of well-defined atomic planes. The shape of the DOS for a-Si is similar to that reported in ref 26. Unlike the interfacial DOS for crystalline Si, the a-Si DOS depicts a damping of >4THz phonon modes and a less drastic DOS reduction of high-frequency modes. Similarly, Carlborg et al.20 found that the lower-frequency phonons of SWCNT are damped out when SWCNTs are surrounded by liquid and solid argon matrices. Such a frequency range corresponded to that of the argon atoms adsorbed onto the SWCNTs, demonstrating a resonant coupling between phases. The results reported so far and the comparison with previous investigations conducted on generic and realistic systems indicate that conclusions drawn from observations of the

Figure 5. Spectral analysis of the normalized heat flux across solid− liquid interfaces and bulk solids: (a) Si(100) facing water, (b) Si(111) facing water, and (c) a-Si facing water.

DOS for generic LJ systems cannot be readily applied to more realistic structures such as Si and a-Si in contact with water. 4.4. Spectral Mapping of the Interfacial Heat Flux. The spectral heat flux, q(ω), was calculated at the solid bulk and interfaces following eq 2. The calculations were normalized with respect to the maximum value of q(ω) in order to obtain a qualitative description of the spectral contributions regardless of the wetting condition. Figure 5(a) depicts q(ω) for the Si(100) interface for hydrophilic and hydrophobic conditions and for a 3 nm slab region in the bulk away from the interface. The shape of q(ω) for the bulk Si follows a similar shape as that of the DOS, with a sharp q(ω) peak at ∼15 THz, indicating a strong contribution of TO modes to the heat transfer.45 Alternatively, large q(ω) peaks at ∼5 THz are observed for the interfacial regions, hydrophobic and hydrophilic, with a subtle frequency shift of the peak position. Unlike previous findings22,24,25 for generic LJ solid−liquid systems, significant q(ω) peaks are observed at high frequencies for the Si(100) interface under hydrophilic conditions and a smaller peak for the hydrophobic interface. The location of these high frequency contributions to q(ω) matches the frequency range where damping of phonon modes at the interface with respect to the 11384

DOI: 10.1021/acs.jpcc.7b01689 J. Phys. Chem. C 2017, 121, 11380−11389

Article

The Journal of Physical Chemistry C bulk was observed in Figure 4(a). According to ref 22, the transmission of these modes can be implied from the damping of the DOS peaks at the interface with respect to the bulk. Figure 5(b) depicts the calculation of q(ω) for bulk and interfacial Si(111) under different wetting conditions. As in the previous case, q(ω) is mostly carried by low frequency phonon modes, and the q(ω) peak is observed in the same region as in Si(100), approximately at ∼5 THz. A subtle down-shift of the q(ω) peak is observed between wetting conditions. Highfrequency modes have a significant contribution to q(ω) for hydrophilic Si(111), and the position of this peak matches the region of the damping and shifting of the high-frequency modes indicated in the DOS in Figure 4(a). As was indicated in ref 22, the damping of the DOS at high frequencies at the interface with respect to the bulk indicates transmission of these modes. Figure 5(c) illustrates q(ω) for the bulk a-Si and hydrophilic (large εSi−O) and hydrophobic (small εSi−O) a-Si interfaces. Similar to the crystalline Si interfaces, the peak contribution to q(ω) is close to 5 THz. Contrariwise to the crystalline interfaces, there is no peak contribution to q(ω) from highfrequency modes, and q(ω) exhibits a minimum upshift in peak from hydrophobic to hydrophilic surfaces. Additionally, at frequencies 10−18 THz there is a slightly higher contribution of these modes for the hydrophilic interface. These observations are in accord with the information deducted from the DOS in Figure 4(b). According to the DOS, there should not exists any q(ω) peak at high frequencies, due to the minimum damping of the DOS at the interface with respect to the bulk. However, only a cumulative calculation of q(ω) as a function of ω can give a clearer picture of the most contributing modes, which is presented in Figure 6. Saaskilahti et al.24 found a major contribution to q(ω) from low-frequency modes and an upshift of the q(ω) peak of the interface with respect to the bulk regions, a similar change to that of the DOS at an interface as indicated in refs 22 and 25 and eventually in the calculation of q(ω). Saaskilahti et al.24

showed that, in general, high-frequency phonon modes do not significantly contribute to the interfacial and bulk heat transfer (decoupling of high-frequency modes), but these increase their contribution as the wetting conditions become more hydrophilic. It must be remarked that the investigations reported in refs 22 and 24 were conducted in generic LJ systems. Some of the findings reported for the generic systems can be observed in the more realistic interfaces investigated herein, such as the generally observed peak contributions to q(ω) by the lowerfrequency modes; however, the high-frequency peaks of q(ω) observed in Si (e.g., in Figure 5a) have never been reported in the generic systems. The relative contributions of these modes with respect to the total are evaluated next. Giri et al.25 found that by increasing the solid−liquid interaction strength or by increasing the wettability of a surface higher-frequency modes contribute to the interfacial heat transfer, and for their particular LJ system, ∼80% of the total heat flow was carried by low-frequency phonons. Figures 6(a) and (b) indicate that for regions in the bulk (crystalline and amorphous) modes with frequencies >15 THz must be considered to reach ∼80% of the total heat flux. However, the interfaces depict a completely different picture. Figure 6(a) indicates that in order to reach ∼80% of the total heat flux at a Si(100) interface, phonon modes of frequencies up to ∼13.5 THz and 10 THz for hydrophilic and hydrophobic surfaces, respectively, have to be considered, and for a Si(111) interface, phonon modes of frequencies up to ∼10.75 THz and 9.25 THz for hydrophilic and hydrophobic surfaces, respectively, have to be considered. As for the amorphous interfaces in Figure 6(b), modes of up to ∼11.5 THz and ∼10 THz must be considered in order to reach ∼80% of the heat flux for hydrophilic and hydrophobic conditions, respectively. Merabia et al.46 used a modified acoustic mismatch model that accounted for a finite bonding strength between two materials. In this model, the combination of the transmission and dynamics of phonons in strongly bonded interfaces provided a simple but elegant explanation for the greater contribution of high frequency modes to heat transfer for hydrophilic surfaces. The following can be summarized from Figure 6: (1) as indicated in refs 24 and 25, increasing the solid−liquid affinity increases the frequency range of phonon modes contributing to the total heat transfer; (2) the gap between the frequency modes contributing to ∼80% of the total heat flux in hydrophobic and hydrophilic surfaces is strongly dependent on the type of surface facing the liquid phase; therefore, the concept of wettability, characterized by the contact angle, does not help in explaining the different spectral contributions to heat transfer; (3) as the surfaces become more hydrophobic, the contribution to the total heat transfer by modes with frequencies