Regular Dynamics Associated with Heat Transfer at the Interface of

1754

J. Phys. Chem. C 2007, 111, 1754-1763

Regular Dynamics Associated with Heat Transfer at the Interface of Model Diamond {111} Nanosurfaces Oleg A. Mazyar,† Tianying Yan,‡ Srirangam V. Addepalli,§ and William L. Hase*,† Department of Chemistry and Biochemistry, Texas Tech UniVersity, Lubbock, Texas 79409-1061, Institute of New Energy Material Chemistry and Institute of Scientific Computing, Nankai UniVersity, China, and High Performance Computing Center, Texas Tech UniVersity, Lubbock, Texas 79409-1061 ReceiVed: August 28, 2006; In Final Form: NoVember 15, 2006

Nonequilibrium molecular dynamics simulations were performed to study the dynamics of heat transfer across the interface of two H-terminated diamond {111} nanosurfaces. It was found that when the surfaces are brought into contact, a coherent low-frequency oscillatory motion, normal to the interface, is initiated for the atomic layers of each surface. This motion lasts more than 50 ps, significantly longer than the relaxation of the interfacial force which occurs on a subpicosecond time scale. Amplitudes of these coherent oscillations gradually increase from the outer atomic layers to the interfacial layer. The temperature of the hot surface does not have a significant effect on the nature of oscillations. Thicker nanosurfaces have oscillations of lower frequency. This coherent oscillation of the surfaces’ atomic layers creates beats in the potential energy of the surface and, thus, in the total energy content of the surface. The amplitude of the energy beat increases as the distance between the interfacial H-atom layers decreases. Thus, the atomic-level energy transfer dynamics from the hot surface to the cold surface is more complex than the exponential dynamics observed for overall heat transfer between the surfaces.

I. Introduction Nanoscale devices, such as miniature electronic integrated circuits and mechanical systems,1 are expected to have a revolutionary impact in all areas of science and engineering.2 To be built and to perform, precise control over the positioning of the nanomachine’s components is needed.3 This requires a fundamental understanding of their large-amplitude internal vibrations, comparable to the size of the system itself. These motions are induced as a result of the device’s operation and affect its performance.4-6 Experimental technologies7-9 to manipulate nanomechanical systems, such as the atomic force microscope (AFM), and to monitor the mechanical motion of the nanosystem, that is, high-resolution scanning (SEM) and transmission (TEM) electron microscopy techniques, provide details of the components’ positions. Theoretical research,4-6 focused on molecular bearings and multiwalled carbon nanotube (MWCN) oscillators,2,4-7,10,11 has contributed to a fundamental, atomic-level description of the large-amplitude motions of nanodevices. Oscillations in molecular bearings were found to adversely affect their performance by creating beats in the shaft angular momentum and in the position of its center of mass.4 Placing a stretching tension on the bearing was found to suppress one of the internal vibrational modes of the bearing and, thus, its beats.4 A rocking motion of the inner tube (shaft) and a wavy deformation in the outer tube (sleeve), undergoing vibrations in the radial breathing mode, were found to be the primary sources of energy dissipation from double-walled carbon nanotube (DWCN) translational degrees of freedom into its disorderly phonon modes.5 This heating of the DWCN na†

Department of Chemistry and Biochemistry, Texas Tech University. Nankai University. § High Performance Computing Center, Texas Tech University. ‡

nooscillator hinders its efficiency and performance. Largeamplitude oscillations, involving motions of all the atoms, are likely to be found in any nanosystem with periodic structure. Thus, it is important to understand the fundamental properties of these large-amplitude regular vibrational motions and their effect on a nanosystem’s performance. As it was noted by Rivera et al.,6 most of the simulations of nanosystems are performed at constant energy, corresponding to simulations in a vacuum, where there is no heat dissipation to the environment. However, any nanosystem will be attached to a surface which can act as a heat sink. Depending on the thermal conductivity of the contact area between the nanosystem and the surface, and the thermal conductivity of both the surface and the nanosystem, heat flow to the surface may or may not be fast enough to maintain the nanosystem at a constant temperature equilibrium condition. It was found that the constant energy and temperature behaviors of DWCN oscillators are qualitatively similar but differ in their quantitative details.6 In this work, we present the results of nonequilibrium molecular dynamics (MD) simulations of heat transfer at the interface of two model diamond {111} nanosurfaces. This research addresses the two above-mentioned issues: (1) the nature and behavior of heat transfer from the nanosystem to the surface and (2) the effect of large-amplitude oscillations of the nanosystem on the dynamics of heat transfer occurring on a time scale comparable to the period of these oscillations. In a previous study, the kinetics of heat transfer across the interface of these two surfaces was studied.12 The coherent, regular-type dynamics considered here was identified but not investigated. This type of motion has been observed in several recent studies of biomolecules and nanostructures.5,13-15 It is of interest to characterize this motion at the interface of nanosurfaces.

10.1021/jp065571d CCC: $37.00 © 2007 American Chemical Society Published on Web 01/10/2007

Regular Dynamics Associated with Heat Transfer

J. Phys. Chem. C, Vol. 111, No. 4, 2007 1755 TABLE 1: Potential Energy Parameters Parameters for the H-Terminated {111} Surfacea C-C stretch Morse function DR βR R0

79.46 kcal/mol 1.858 Å-1 1.544 52 Å C-H stretch Morse function

Dr βr r0

109.94 kcal/mol 1.852 Å-1 1.095 45 Å

diagonal quadratic potential parameters Figure 1. Depiction of the upper (hot) and larger, lower (cold) H-terminated {111} diamond surfaces. The dimensions of each surface are for their 0 K, equilibrium geometries.

II. Computational Procedure A. Model of Contacting Nanosurfaces and Potential Energy Function. Details of the model of the two contacting H-terminated {111} diamond nanosurfaces, the potential energy function for the nanosurfaces, and the nonequilibrium MD simulation procedure were described previously.12 As shown in Figure 1, the upper nanosurface, with a height of 19.65 Å and an interfacial area of 20.19 × 21.86 Å2 at 0 K, sits on a lower surface of the same height, but having a larger area of 60.58 × 65.59 Å2. Each surface consists of 20 layers, that is, 19 carbon layers and 1 hydrogen layer at the interface. The outermost C-atom layer of each surface is held rigid, so that the separation of these two layers determines the normal load applied to the surfaces. The potential energy function for the system is written as

V ) Vsurf,u + Vsurf,l + Vinter

(1)

where Vsurf,u and Vsurf,l are the potentials of the upper and lower surfaces and Vinter is their interaction potential. In this work, Vsurf is the harmonic valence force field potential developed by Tubino, Piseri, and Zerbi16 and Hass et al.17 to fit the diamond phonon spectrum, with the modification that the C-C and C-H stretches are represented by the Morse functions18

V(R) ) DR[1 - exp{-βR(R - R0)}]2

(2)

V(r) ) Dr[1 -exp{-βr(r - r0)}]2

(3)

and

respectively, where DR is the tertiary C-C bond dissociation energy of 79.46 kcal/mol, βR ) (fR/2DR)1/2 ) 1.858 Å was determined from the harmonic force constant (fR) of 3.812 mdyn/ Å, and R0 is the equilibrium C-C bond length of 1.544 52 Å. Similarly, the parameters of the C-H stretch potential given by eq 3 are Dr ) 104.94 kcal/mol, βr ) 1.852 Å, and r0 ) 1.095 45 Å. All parameters of the surface potential are represented in Table 1. Vinter in eq 1 describes the nonbonded intermolecular potential between the top and bottom surfaces. Only interactions between the terminal H-atoms and the C-atoms of the first and second interfacial layers of the two surfaces are included in Vinter, since the distances between other interfacial atoms are too large to contribute to Vinter. These nonbonded H-H, C-H, and C-C potentials are represented by the EXP-6 function of Williams and Starr,19 that is

θ0 φ0 fθ fφ

109.471° 109.471° 0.868 mdyn‚Å/rad2 0.725 mdyn‚Å/rad2 nondiagonal quadratic potential parameters

fRRh fRθ fθθb fθθc

0.163 mdyn/Å 0.39 mdyn/rad 0.177 mdyn‚Å/rad2 -0.0149 mdyn‚Å/rad2

Parameters for the Surface-Surface Intermolecular Potentiald H‚‚‚H H‚‚‚C C‚‚‚C

a

b

c

2790.87 15651.29 87774.86

3.74 3.67 3.60

-32.50 -136.95 -576.96

a The coordinates R, r, θ, and φ are C-C stretch, C-H stretch, C-C-C bend, and C-C-H bend, respectively. b Nondiagonal force constant for two different C-C-C bends, which share one stretch coordinate but not a central atom. c Nondiagonal force constant for two different C-C-C bends, which share one stretch coordinate and a central atom. d The intermolecular potential energy parameters a, b, and c are in units of kcal/mol, Å-1, and kcal‚Å6/mol, respectively.

V(r) ) a exp(-br) +

c r6

(4)

which were derived to represent nonbonded interactions for hydrogen and carbon atoms in experimentally determined crystal structures of 18 hydrocarbon molecules. Parameters for the H-H, C-H, and C-C intermolecular potentials are listed in Table 1. For short distances, the attractive component of eq 4 overwhelms the repulsive component and the potential becomes unphysically attractive. This possibility was checked for the current simulation, and each of the three intermolecular potentials remained repulsive at the shortest nonbonded atomatom distances considered. The closest separation between the H-atom planes of the two surfaces is 0.75 Å, and for the interfacial geometry considered for the simulations (see next section), the resulting shortest H-H separation between two H-atoms of the surfaces is 1.637 Å. The H-H potential energy for this separation is positive and repulsive. The total interfacial potential for an interfacial separation of 0.75 Å is ∼2200 kcal/ mol at 0 K (see Figure 3 of ref 12). B. Molecular Dynamics Simulations. The MD simulations were performed by solving Hamilton’s classical equations of motion with the VENUS computer program.20 The numerical integrations were initiated with a Runge-Kutta-Gill algorithm and then completed with a sixth-order Adams-Moulton routine with an integration time step (∆t) of 0.1 fs. Each nanosurface was equilibrated at the desired temperature before they were brought into contact. The upper surface was

1756 J. Phys. Chem. C, Vol. 111, No. 4, 2007

Mazyar et al.

Figure 2. First and last 5 ps of the 50 ps motion of representative layers of the isolated upper surface, after equilibration to 800 K.

Figure 3. First and last 5 ps of the motion of representative layers of the upper and lower surfaces, after they are brought into contact with Seq ) 0.75 Å. The average interfacial force is 91 nN, and the initial temperatures of the upper and lower surfaces are 800 and 300 K, respectively. For the upper surface, the black line represents the motion of the 4th atomic layer, the red line, the eighth layer, the green line, the 12th layer, the blue line, the 16th layer, and the cyan line, the 20th interfacial H-atom layer. For the lower surface, the blue line represents the motion of the 8th atomic layer, the green line, the 12th layer, the red line, the 16th layer, and the black line, the 20th interfacial H-atom layer. All layers are enumerated starting from the outermost rigid layer of the corresponding surface.

equilibrated at temperatures of Thot ) 800, 1000, or 1200 K and the lower surface at 300 K. During the first step of the MD equilibration, all nonrigid atoms of a surface were coupled to two Berendsen baths21 of a given temperature. One thermal bath was connected to the surface’s carbon atoms, and the second one was connected to the hydrogen atoms to provide fast thermal

equilibration of the surface at the desired temperature, as energy exchange between the C-H stretch modes and the other surface modes is a slow process. This MD simulation was continued until the potential and kinetic energies of the surface become equal. After this equilibration step, the upper surface was further equilibrated for 10 ps to ensure the temperatures of the


J. Phys. Chem. C, Vol. 111, No. 4, 2007 1757

individual surface atomic layers are the same and equal to the desired temperature.22 The lower surface was further equilibrated for 30 ps with its lowest two C-atom layers (excluding the outermost rigid layer) coupled to a 300 K Berendsen thermal bath, since this bath is included in the simulations of the heat transfer from the hot, upper to the cold, lower nanosurface. With their outermost atomic layers rigid, the surfaces were brought into contact and the nonbonded interactions between their interfacial atomic layers turned on. The position of the upper surface was selected so that, at 0 K, each H-atom of the upper surface sat over the center of the equilateral triangle formed by three H-atoms of the lower surface. The outermost layers of the upper and lower surfaces are held rigid. Fixing their separation defines the distance (S) between the planes of the interfacial H-atom layers. When the interactions between the upper and lower surfaces are absent and each surface is in its 0 K equilibrium geometry, this distance is identified as Seq. Values of Seq of 0.75, 1.25, 1.75, 2.25, and 2.75 Å were considered here. Previous work has shown that at Thot ) 800 K these values of Seq give average values of S of 1.31, 1.52, 1.80, 2.19, and 2.67 Å and average interfacial forces of 91.16, 43.22, 10.89, -3.73, and -5.90 nN.12 For Thot ) 1000 and 1200 K and Seq ) 1.25 Å, the average value of S is 1.51 Å, and only slightly different than the S value of 1.52 Å at 800 K. Similarly, average interfacial forces are slightly different than the 800 K value and are 43.43 and 43.72 nN, respectively. The two lowest C-atom layers (excluding the rigid layer) of the lower surface were connected to a 300 K Berendsen thermal bath providing a constant temperature heat sink for the lower surface as heat is transferred from the hot, upper surface. To study the energy transfer dynamics from the hot to cold nanosurfaces, a 52.5 ps MD simulation of energy transfer from the hot to cold model diamond {111} nanosurfaces was performed. C. Simulation Analyses. The simulations were analyzed for regular, oscillatory dynamics in the motion of each nanosurface and in the heat transfer at their interface. Fixed in its 0 K geometry, the outermost (first) layer of C-atoms of the upper surface was placed in the x,y-plane of the coordinate system. When the upper surface is in its 0 K equilibrium geometry, and in the absence of interactions between the upper and lower surfaces, the z-coordinate for the x,y-plane of a C-atom layer is defined as Z0J . The rigid, outer layer of the upper surface is identified as J ) 1. For T * 0 K, the average z-coordinate of the NJ atoms in layer J at time t is NJ

ZJ(t) )

zi(t)/NJ ∑ i)1

(5)

The value of Z0J - ZJ (t) allows convenient monitoring of the collective motion of atoms in layer J of the upper surface. If this value is negative, layer J is compressed in comparison to its 0 K geometry. Further decrease of this value means additional compression of the layer. For the lower surface, the convenient quantity for monitoring motion of the atomic layers is ZJ(t) Z0J , where the J ) 1 layer is the rigid outer layer of the lower surface. A negative value of ZJ(t) - Z0J means that the Jth layer of the lower surface is compressed with respect to its 0 K equilibrium geometry. Vibrational frequencies for the interlayer motions of a surface were determined by taking the Fourier transform23 of the 〈∑J ZJ(t)∑J ZJ(0)〉 autocorrelation function, that is,

I(f) )

∫-∞∞〈∑ZJ(t)∑ZJ(0)〉e2πift dt J

(6)

J

The fast-Fourier-transform (FFT) algorithm used for this analysis was taken from work by Press et al.23 III. Simulation Results A. Interlayer Vibrational Dynamics and Frequencies for the Isolated Upper Nanosurface. Before studying the atomiclevel dynamics of the surfaces during heat transfer, the dynamics of the isolated hot, 19.65 Å thick, upper surface (Figure 1) was investigated. Figure 2 depicts the first and last 5 ps of motion of this surface during a 50 ps long trajectory, after equilibration to 800 K. The motion of each of the three representative layers is periodic, with a vibrational frequency of 74 cm-1. There are two important features of this motion. First, it is periodic and not chaotic. At 1200 K (not shown), the motion is also periodic, with the same frequency and slightly larger amplitude. It is possible that the motion becomes chaotic at higher temperatures, but this was not investigated. The second important feature is the observation of only one dominant vibrational motion for the interlayer dynamics. It is a symmetric motion, with each of the layers moving up or down coherently. It is of interest to determine whether this motion corresponds to a normal mode of the upper surface or is a stable superposition of the surface’s normal modes, which relaxes very slowly.24 To answer this question, the surface’s normal modes of vibration were determined by diagonalizing the surface’s 5238 × 5238 mass-weighted Cartesian force constant matrix.25 This analysis gives a symmetric normal mode of vibration with all of the carbon atoms moving coherently in the same direction, and the frequency for this mode is 74.8 cm-1, nearly identical to the frequency obtained for this motion from the simulation. Thus, apparently the symmetric motion present in the simulation is a normal mode of vibration. The demands on both memory and cpu time make the above matrix diagonalization approach impractical for a much larger surface; that is, the “cost” of diagonalizing a N × N matrix scales by approximately N4. To study the symmetric, coherent vibration for a much larger diamond surface, a model for the above 19.65 Å thick upper surface was constructed consisting of a linear chain of “united atoms”,26 which represent the H-atom layer, the carbon bilayers, and the terminal, rigid carbon layer. For this model of the diamond {111} nanosurface, the 19 C-atom layers and 1 H-atom layer are represented by a linear chain of 11 united atoms, with the end united atoms representing the outermost rigid C-atom and the interfacial H-atom layers and the remaining 9 united atoms representing the interior C-atom bilayers. The mass of each united atom, in this model, is the total mass of the atoms in an interior bilayer or an end layer. Harmonic potentials were determined between each pair of interior adjacent bilayers, the end C-atom layer and its adjacent bilayer, and the interfacial H-atom and its adjacent bilayer by making small z-separation displacements of each of the bilayer/bilayer and bilayer/layer separations and fitting their resulting potential curves with V ) f∆z2/2, where ∆z is the displacement from the 0 K, equilibrium separation. A normal mode analysis, using these force constants, was performed for this linear chain, and a frequency of 85.7 cm-1 was obtained for the symmetric, coherent motion in which all the united atoms move in-phase. This frequency is close to the 74 cm-1 value found from the simulation. A hot surface with a 25.83 Å thickness is also considered in this work, and it is of interest to determine how the frequency


Mazyar et al.

Figure 4. Dynamics of the first (left side) and last (right side) 5 ps of a trajectory with Seq ) 0.75 Å and initial temperatures for the upper and lower surfaces of 800 and 300 K, respectively. Plotted are the total, potential, and kinetic energies of the hot, upper surface and the surface thickness.

for this coherent motion varies with surface thickness. This thicker surface has 25 C-atom layers and the H-atom interfacial layer, and the above linear chain model gives a frequency of 65.2 cm-1 for the coherent motion, which is 85.7/65.2 ) 1.3 times smaller than that for the 19.65 Å thick surface. As discussed below in section IV.C, this is the difference found from the simulations for the two surfaces. Application of this model to chains with 49 and 97 C-atom layers gives frequencies

of 33.3 and 16.9 cm-1. Thus, the frequency of this coherent motion is predicted to decrease in a near-linear manner with an increase in surface thickness. B. Dynamics of the Nanosurfaces with an Interfacial Force. As discussed above, the hot, upper surface undergoes a coherent oscillatory motion when it is isolated from the cold, lower surface. It is of interest to study the interfacial motions of the upper and lower surfaces when they are brought into



Figure 5. Dynamics of the hot, upper surface for different interfacial separations (Seq), with Thot ) 800 K. The upper graph is the thickness of the hot, upper surface versus time and the lower the total energy content of this surface versus time.

Figure 6. Frequency spectrum (eq 6) of the hot, upper surface for different applied loads. The loads of 0, 43, and 91 nN correspond to an isolated surface, Seq ) 1.25 Å, and Seq ) 0.75 Å, respectively.

contact, in the presence of an interfacial force, and heat transfer occurs. Fixing the distance between the outermost layers of each surface determines the initial distance (S) between the H-atom interfacial layers of the surfaces and, thus, the initial interfacial force. Presented here are the dynamics of the surfaces when

they are brought into contact with a separation of their outermost layers so that Seq ) 0.75 Å. The upper surface has a temperature of Thot ) 800 K, and the lower surface temperature is 300 K. As discussed above, in section II.B, with Thot ) 800 K and Seq ) 0.75 Å, the average separation between the interfacial H-atom layers becomes 1.31 Å. As shown in Figure 3, once the simulation begins, the interfacial force rapidly relaxes with a decrease in its value,12 resulting in compression of both surfaces and an increase in the interfacial separation (S). The relaxation of this compression leads to the coherent oscillatory motion of the individual atomic layers illustrated by the ZJ,eq - ZJ(t) and ZJ(t) - ZJ,eq plots for the upper and lower surfaces, respectively, in Figure 3. At t ) 0, both surfaces are thermally expanded so that their respective values of ZJ,eq - ZJ(t) and ZJ(t) - ZJ,eq are positive. The initial repulsive interfacial force compresses both surfaces, with a maximum compression at ∼0.2 ps. The maximum compression of the upper surface reaches ∼0.7 Å, while that for the lower surface is only ∼0.4 Å. Then, both surfaces start simultaneously expanding. As shown in Figure 3, the process of compression and expansion of both surfaces is periodic. The expansion and compression peaks of the upper and lower surfaces coincide; that is, the layers of the upper surface move in counter-phase


Mazyar et al.

to the motion of the atomic layers of the lower surface. At the same time, the atomic layers within each surface move in-phase. The energy associated with relaxation of the interfacial force is rapidly transferred to the layers of each surface, so that the amplitudes of the layers in a surface become nearly equivalent. Figure 3 shows that the coherent oscillations of the atomic layers, within each of the surfaces, are present up to the end of the 52.5 ps long MD trajectory. At the end of the trajectory, the amplitudes of the upper surface’s layers are somewhat smaller than those for the first 5 ps, as a result of heat transfer; for example, for the 16th and 20th layers, the difference, in the maximum and minimum of the amplitudes, is about 0.05 Å smaller at the end of the trajectory. A major difference between the motions of the atomic layers during the first 5 ps and the last 5 ps of the 52.5 ps long trajectory is that by the end of the simulation all atomic layers of both the upper and lower surfaces are oscillating in-phase; that is, when the upper surface is expanded, the lower one is compressed and vice versa. C. Oscillatory Dynamics of Energy Transfer from the Hot Surface to the Cold Surface. From the previous study,12 it was found that the overall energy transfer from the hot surface to the cold surface is a first-order process, represented by the rate law

-

d[E(t) - Ef] ) k[E(t) - Ef] dt

(7)

where E(t) is the energy content of the upper (hot) surface at time t and Ef is its energy content after equilibration between the upper and lower surfaces is attained. Integration of this equation gives

E(t) ) (Ei - Ef) exp(-kt) + Ef

(8)

An explicit assumption of this model, which provides an overall excellent representation of the simulation results, is that the rate constant of energy transfer is independent of the excess energy content, that is, [E(t) - Ef], of the hot, upper surface. A detailed investigation of the energy transfer from the upper surface shows that the dynamics of the trajectory depicted in Figure 3 is considerably more complex than that described by the above first-order, exponential decay model. The short-time dynamics of energy transfer from the upper surface, preheated to 800 K with the lower surface at 300 K, is shown on the left side of Figure 4. For this trajectory, the 0 K equilibrium interfacial distance (Seq) is 0.75 Å. When the simulation is initiated and the interfacial potential between the two surfaces is turned on, an additional energy of 3220 kcal/mol is added to the system. This energy is initially localized at the interface but is quickly distributed to both surfaces. The initial energy of the upper surface, equilibrated at 800 K, is 7864 kcal/mol. However, after the interfacial potential is turned on, the energy of the upper surface reaches 9489 kcal/mol at 0.2 ps. This additional energy is primarily potential (see Figure 4) arising from surface compression as a result of the interfacial force. In the initial 0.2 ps of the simulation, the thickness of the upper surface is compressed from 19.8 to 19.0 Å. The instantaneous “turning on” of the interfacial force and ensuing compression of the surfaces gives rise to coherent oscillations in the separations between the surfaces’ atomic layers, as discussed above and illustrated in Figure 3. However, as shown on the left side of Figure 4, there is also a coherent oscillation in the thickness of the upper surface, which gives rise to a beat in the surface’s potential energy. Each peak in the surface’s potential energy corresponds to a minimum in its

thickness. This coherent oscillation in the surface’s potential energy and thickness is maintained out to 52.5 ps when the trajectory is terminated, and the stability of this motion suggests it lasts much longer. There is not a periodic oscillation in the kinetic energy, as found for the potential energy. The average potential energy is ∼200 kcal/mol higher than the kinetic energy, which results from the compression of the upper surface. With Seq ) 0.75 Å and the interfacial potential “turned on”, the 0 K potential minimum of the upper surface has an energy of 193 kcal/mol. The beats in the potential energy combined with the smoothly decreasing kinetic energy of the upper surface lead to the beats in the total energy of the upper surface as it cools down. Thus, the hot, upper surface periodically acquires energy as it cools. Beats in the potential energy and, thus, the total energy appear as long as there are coherent oscillations of the upper surface’s layers. As shown in Figure 4, the amplitude of the beat in the total energy exceeds 60 kcal/mol at t ) 50 ps. IV. Discussion The above simulations illustrate coherent, periodic dynamics as heat is transferred from hot to cold diamond nanosurfaces following their contact in the presence of an applied force. The simulations were performed with a separation of the outermost layers of the surfaces so that the 0 K separation of the H-atom interfacial layers is Seq ) 0.75 Å, the temperature of the hot surface is Thot ) 800 K, and the 0 K thickness of this surface is 19.65 Å. Below, the regular dynamics associated with the heat flow is investigated in more detail by varying the interfacial separation (Seq), the temperature of the upper surface (Thot), and the thickness of this surface. A. Effect of Interfacial Separation. For each simulation, the separation of the outer layers of the two surfaces is held fixed, that is, specified by the 0 K equilibrium interfacial distance (Seq). Figure 5 illustrates the effect of the interfacial separation on the thickness of the upper surface versus time. For Seq < 2.0 Å, the initial interfacial force is repulsive,12 so both surfaces start compressing when the interfacial potential is turned on. The initial compression of the upper surface increases from 0.2 to 0.8 Å as Seq decreases from 1.75 to 0.75 Å. However, after ∼1 ps, the amplitudes of the oscillations become almost similar for all Seq. For the large interfacial separation of 2.75 Å, the initial interfacial force is attractive12 and there is a small initial 0.02 Å expansion of the upper surface. As the interfacial separation (Seq) is increased, the frequency of the oscillation in Figure 5 becomes slightly lower. The short-time (i.e., 0-5 ps) dynamics of energy transfer from the hot surface is illustrated in Figure 4 for Seq ) 0.75 Å. This dynamics is shown in the lower graph of Figure 5 for the larger Seq values of 1.25, 1.75, 2.25, and 2.75 Å. It is seen that decreasing the initial interfacial separation and, thus, increasing the interfacial interaction enhances heat transfer from the hot surface to the cold surface.1 The amplitude of the beats in the total energy of the hot surface decreases with an increase in the interfacial separation (Seq) At the same time, the amplitudes of the oscillations of surface’s thickness are almost similar after ∼1 ps for all Seq, that is, under different compression conditions of the surfaces. This effect results from an increasingly higher interfacial energy added to the system when the simulations are initiated at smaller equilibrium interfacial separation (Seq) (see section III.C). The frequency spectrum for the layers of the hot, upper surface (eq 6) was calculated for the surfaces not in contact and an Seq value of 0.75 and 1.25 Å, to quantify the interlayer



Figure 7. Effect of temperature (Thot ) 800 and 1200 K) on the dynamics of the hot, upper surface; Seq ) 1.25 Å. Upper graphs, motion of representative surface layers; lower graphs, total energy of the surface versus time.

frequencies versus interfacial separation. As discussed in section II.B., the Seq values of 0.75 and 1.25 Å correspond to applied loads of 91 and 43 nN for Thot ) 800 K. The calculated spectra are shown in Figure 6. With no load, the surface has a single peak at ∼74 cm-1. With a 43 nN load, the principal peak is shifted to ∼87 cm-1 and there are now side peaks, ranging from ∼75 to 115 cm-1. Increasing the load shifts the peak further to ∼89 cm-1 with the intensity of the side peaks, ranging from ∼75 to ∼120 cm-1, becoming stronger. B. Effect of Surface Temperature. Simulations were compared with the hot, upper surface at 800 and 1200 K to determine how the temperature of this surface affects the coherent oscillatory motion of the surface’s layers. The interfacial separation (Seq) was set at 1.25 Å. Values of ZJ,eq - ZJ(t) are plotted in the upper graphs of Figure 7 for the two temperatures. The higher the temperature of the surface, the greater its thermal expansion,12 giving rise to a somewhat smaller initial separation between the interfacial H-atom layers of the two surfaces. Figure 7 shows this does not have a significant effect on the surface’s oscillatory motion and the interlayer dynamics is very similar at 800 and 1200 K. Previous work12 has shown that the rate constant for heat transfer from the hot to cold surfaces (eq 7) is independent of the hot surface’s temperature. As discussed above (i.e., section III.C and Figures 4 and 5), there are beats in the total energy of the upper surface associated with the coherent oscillatory motion

of the surface’s layer. The lower graphs of Figure 7 show that the frequency of this beat is nearly the same at Thot values of 800 and 1200 K. In addition, the amplitudes of the beats at these two temperatures are similar. Increasing the temperature of the hot surface from 800 to 1200 K does not have an important effect on its regular dynamics. C. Effect of Surface Thickness. For the above simulations, a thickness of 19.65 Å was used for the hot, upper surface. To determine whether increasing the surface’s thickness affects the coherent oscillatory motion of the atomic layers, a simulation was performed with a 25.83 Å thick upper surface. The interfacial area of this surface remained at 20.19 × 21.86 Å2, and no changes were made to the size of the lower surface. This larger, upper surface contains 2274 atoms compared to the 1746 atoms for the standard surface model and has 25 C-atom layers in addition to the interfacial H-atom layer. The simulation was performed for the upper surface preheated to 800 K, and the Seq value was 1.25 Å. Representative values of ZJ,eq - ZJ(t) for these simulations (see sections II.C and III.A) are plotted in Figure 8. The frequency of the coherent oscillatory motion is lower in comparison to the frequency of this oscillation for the surface with a thickness of 19.65 Å, occurring under the same Thot and Seq conditions (i.e., see Figure 7). The frequency is approximately 87 and 66 cm-1, respectively, for the surfaces with thicknesses of 19.65 and 25.83 Å. The ratio of these frequencies,


Mazyar et al.

Figure 8. Dynamics of the hot, upper surface with a thickness of 25.83 Å, Seq ) 1.25 Å, and Thot ) 800 K. Upper graphs, motion of the surface’s layers for the initial and final 5 ps; lower graphs, total energy content of the surface for the initial and final 5 ps.

that is, 87/66 ) 1.3, is the same as that found from the linear chain model, discussed above in section II.A. The amplitudes of the layer oscillations are similar for the two surfaces, and this is shown for the first 5 ps of the simulations in Figures 7 and 8. As shown in Figure 8, coherent oscillations of the atomic layers remain intense during the last 5 ps of the 52.5 ps simulations with the amplitude of the interfacial H-atom layer reaching 0.05 Å, which is somewhat smaller than the approximately 0.10 Å for the first 5 ps. This decrease in amplitude versus time is akin to what is found for the smaller 19.65 Å thick, upper surface (see section III.B) and arises from the heat transfer. Similar to what is found for the 19.65 Å thick surface, there are also beats in the energy content of the thicker, upper surface lasting to at least 52 ps (see Figure 8). The amplitudes of these beats for the thicker surface are nearly the same as those for the 19.65 Å thick surface. The surface’s thickness might have a significant effect on the long-time behavior of energy transfer. As shown in Figure 8, the interfacial H-atom layer of the thicker surface is shifted by ∼0.2 Å from its equilibrium position at 0 K in the absence of the lower surface. This shift was found to be only ∼0.18 Å for the surface with a thickness of 19.65 Å. This means that the average distance between the interfacial H-atom layers of the hot, upper and cold, lower surfaces is larger for the thicker upper surface. Trajectories show that the average distance between the H-atom layers of the upper and lower surfaces is

1.54 Å for the 25.83 Å thick surface and 1.52 Å for the surface with the 19.65 Å thickness. This difference in the distances between the surfaces results in an interfacial force of 40 nN acting between the surfaces in the case of the thicker upper surface. At the same Seq value, the interfacial force was found to be 43 nN for the thinner surface of 19.65 Å. These effects are expected to be even more significant at smaller Seq, when the difference in thicknesses of the hot, upper surfaces is larger. Since it was found that the rate constant for the energy transfer from the hot to cold surfaces strongly depends on the interfacial force acting between the two surfaces in contact,12 the effect of the thickness of the surfaces on the rate constant for the energy transfer will be studied in future work. V. Summary In the work presented here, nonequilibrium molecular dynamics simulations were performed to study the dynamics of heat transfer across the interface of two H-terminated diamond {111} nanosurfaces. The model used in this study is one in which a small, hot nanosurface is placed on a large, cold surface coupled to a thermal bath. The important finding of the study is that when the surfaces are brought into contact, a coherent low-frequency oscillatory motion, normal to the interface, is initiated for each surface. This motion lasts more than 50 ps, significantly longer than the relaxation of the interfacial force which occurs on a subpicosecond time scale.12 With time, the

Regular Dynamics Associated with Heat Transfer amplitudes of these coherent oscillations gradually increase from the outer atomic layers to the interfacial layer. The temperature of the hot surface does not have a significant effect on the nature of oscillations. Nanosurfaces of larger size have oscillations of lower frequency. This coherent oscillation of the surface layers creates a beat, with the same coherent frequency, in the potential energy of the surface and, thus, in the total energy content of the surface. The amplitude of the energy beat increases as the distance between the interfacial H-atom layers decreases. Thus, the energy transfer dynamics from the hot surface to the cold surface is more complex than the exponential of eq 8, which gives an overall description of the heat transfer. It is of interest to determine whether this periodicity in the energy content of the hot surface can be harnessed by a nanodevice. In addition, it is important to determine the generality of the coherent dynamics observed here for diamond nanosurfaces by investigating other materials. Some may have significantly more chaotic dynamics than that for diamond nanosurfaces. The rigid layers used in this simulation study, for the boundaries of the upper and lower surfaces, are models for heterogeneous interfaces with other materials. In a more detailed study, one would want to explicitly include these materials and their intermolecular interactions with the diamond surfaces. Such a study would supplement the current work by providing information regarding heat transfer at such heterogeneous boundaries and the nature of the diamond vibrations at the interfaces. It is of interest that regular dynamics, similar to those observed here, have been detected and imaged for a homogeneous domain within a heterogeneous material.13 Finally, the times required for the coherent, symmetric vibration to undergo one period of motion and for sound to traverse the upper surface in both the forward and reverse directions are quite similar. The frequency for the coherent oscillation for the isolated 19.65 Å thick surface is 74 cm-1, which corresponds to a 0.45 ps vibrational period. The vibrational period for the thicker 25.83 Å surface is estimated as 0.59 ps (see section III.A). What is interesting is that these times are quite similar to those required for sound to traverse these surface slabs. The speed of sound in diamond may be estimated27 as c ) (E/F)1/2, where E is Young’s modulus and F the density. Young’s modulus for the {111} diamond surface is 1223 GPa28 and the density of diamond at 300 K is 3515.25 kg/m,3,29 which gives 18 652 m/s for the speed of sound. In comparison, reported values for the speed of sound vary from 12 00030 to 18 190 m/s.31 Using these latter two values, the times required for sound to traverse the 19.65 and 25.83 Å surfaces in the forward and reverse directions are in the ranges 0.220.33 and 0.28-0.44 ps, respectively. The implication from this analysis is that the coherent, symmetric vibration observed in this study may be instrumental in propagating sound in diamond, an issue to investigate in more detail in future work.

J. Phys. Chem. C, Vol. 111, No. 4, 2007 1763 Acknowledgment. This research was supported by the Office of Naval Research, Department of Energy, the Robert A. Welch Foundation, and the High Performance Computing Center (HPCC) at Texas Tech University. References and Notes (1) Forro, L. Science 2000, 289, 560. (2) Drexler, K. E. Nanosystems: Molecular Machinery, Manufacturing, and Computation; Wiley-Interscience: New York, 1992. (3) Drexler, K. E. Trends Biotechnol. 1999, 17, 5. (4) Tuzun, R. E.; Noid, D. W.; Sumpter, B. G. Nanotechnology 1995, 6, 64. (5) Zhao, Y.; Ma, C.-C.; Chen, G. H.; Jiang, Q. Phys. ReV. Lett. 2003, 91, 175504. (6) Rivera, J. L.; McCabe, C.; Cummings, P. T. Nanotechnology 2005, 16, 186. (7) Cumings, J.; Zettl, A. Science 2000, 289, 602. (8) Yu, M.-F.; Lourie, O.; Dyer, M. J.; Moloni, K.; Thomas F. Kelly, T. F.; Ruoff, R. S. Science 2000, 287, 637. (9) Yu, M.-F.; Yakobson, B. I.; Ruoff, R. S. J. Phys. Chem. B 2000, 104, 8764. (10) Merkle, R. C. Nanotechnology 1993, 4, 86. (11) Zheng, Q.; Jiang, Q. Phys. ReV. Lett. 2002, 88, 045503. (12) Mazyar, O. A.; Hase, W. L. J. Phys. Chem. A 2006, 110, 526. (13) Bargheer, M.; Zhavoronkov, N.; Gritsai, Y.; Woo, J. C.; Kim, D. S.; Woerner, M.; Elsaesser, T. Science 2004, 306, 1771. (14) Moritsugu, K.; Miyashita, O.; Kidera, A. Phys. ReV. Lett. 2000, 85, 3970. (15) Ma, C.-C.; Zhao, Y.; Chen, G.-H.; Jiang, Q. Presented at the Conference on Theory and Applications of Computational Chemistry, Gyeongju, Korea, Feb 2004. (16) Tubino, R.; Piseri, L.; Zerbi, G. J. Chem. Phys. 1972, 56, 1022. (17) Hass, K. C.; Tamor, M. A.; Anthony, T. R.; Banholzer, W. F. Phys. ReV. B 1992, 45, 7171. (18) de Sainte Claire, P.; Barbarat, P.; Hase, W. L. J. Chem. Phys. 1994, 101, 2476. (19) Williams, D. E.; Starr, T. L. Comput. Chem. 1977, 1, 173. (20) Hase, W. L.; Duchovic, R. J.; Hu, X.; Komornicki, A.; Lim, K. F.; Lu, D.-h.; Peslherbe, G. H.; Swamy, K. N.; Vande, Linde, S. R.; Varandas, A.; Wang, H.; Wolf, R. J. Quantum Chemistry Program Exchange (QCPE) Bulletin 1996, 16, 671. (21) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (22) The temperature of the atoms in a layer was calculated by dividing the total kinetic energy of the layer’s atoms by 3NkB, where N is the number of atoms in the layer and kB is Boltzmann’s constant. (23) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in C, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1999. (24) Baer, T.; Hase, W. L. Unimolecular Reaction Dynamics. Theory and Experiments; Oxford: New York, 1996; p 67. (25) Califano, S. Vibrational States; Wiley: New York, 1976. (26) Yan, T.; Hase, W. L. J. Phys. Chem. A 2001, 105, 2617. (27) Pain, H. J. The Physics of Vibration and WaVes; Wiley: London, 1968. (28) Wang, S.-F.; Hsu, Y.-F.; Pu, J.-C.; Sung, J. C.; Hwa, L. G. Mater. Chem. Phys. 2004, 85, 432. (29) Singh, J. Physics of Semiconductors and Their Heterostructures; McGraw-Hill: New York, 1993. (30) http://hyperphysics.phy-astr.gsu.edu/hbase/tables/soundv.html (HyperPhysics, Carl R. (Rod) Nave, Department of Physics and Astronomy, Georgia State University, Atlanta, Georgia 30303-3088). (31) Saviot, L.; Murray, D. B. Phys. ReV. Lett. 2004, 93, 05506.

Regular Dynamics Associated with Heat Transfer at the Interface of

Recommend Documents