Frictional Energy Dissipation due to Phonon Resonance in Two-Layer Graphene System

The frictional energy dissipation mechanism of a supported two-layer graphene film under the excitation of the model washboard frequency is investigated by molecular dynamics simulations. The results show that two local maxima in the energy dissipation rate occur at special frequencies as the excitation frequency increases from 0.1 to 0.6 THz. By extracting the vibrational density of states of the graphene, it is found that large numbers of phonons with frequencies equal to the excitation frequency are produced. A two-degree of freedom mass-spring model is proposed to explain the molecular dynamics results. Since the washboard frequency for atomically surfaces in wearless dry friction can be analogous to the excitation frequency in the molecular dynamics simulations, our study indicates that the phonon resonance would occur once the washboard frequency is close to the natural frequency of the frictional system, leading to remarkable local maxima in energy dissipation.


Introduction
Friction is a common phenomenon that occurs in aspects of our lives. Nevertheless, it is difficult to give a universal understanding of the energy dissipation mechanism due to its complexity. For example, the friction generally shows a gradually decreasing trend with the increase of sliding velocity on the macro-scale [1], while on the atomic scale the friction usually presents a trend of gradual increase with the increase in the sliding velocity [2,3]. The variation of friction is closely related to the energy dissipation mechanism in the friction process [4]. For a long time, the researchers generally use the potential barrier height [5,6] and shear strength [7,8] of friction interface to compare and quantify the magnitude of friction force when two interfaces slide relative to each other. The maximum friction generally implies a unique energy dissipation channel [9]. However, the fundamental physics involved in the friction process, such as the phonon [10][11][12][13] and/or electron [14,15] related dissipation, remains unclear.
The friction process can excite large numbers of nonequilibrium phonons (elastic waves) [16], which dissipates the mechanical kinetical energy through the phonon scattering [11,17,18] with other phonons, boundaries, and/or impurities. As for the two-wall carbon nanotube oscillators, Tangney et al. [19] showed that the friction between the inter-tube and the out-tube would have a maximum value when the group velocity of excited phonons was equal to the sliding velocity. Panizon et al. [20] deduced a formula to calculate the friction force based on linear response theory. Their results show that the resonance would occurs and can cause a local maximum in the friction force when the group velocity of the excited phonons is equal to the phase velocity. Our recent theoretical and experimental studies [21] also suggest that when the atomic force microscope (AFM) tip slides on the atomically flat surface, the phonon mode of the substrate excited by the slider also resonates with the entire friction system, resulting in multiple local maxima in the friction force with the increasing sliding velocity. The similar resonance principle is also adopted by quartz crystal microbalance (QCM) to measure the energy dissipation rate of adsorbed films [22] and two-dimensional materials [23]. Although it is common sense that a physical resonance can lead to a local maximum in energy dissipation, it would be interesting to relate resonance under the washboard frequency to the energy dissipation mechanism of wearless dry friction, because the reaction rate theory [3], not the resonance, is often used to explain the energy dissipation mechanism traditionally. Therefore, we speculated that the excited phonons will be strengthened and the energy dissipation will increase if the excited phonons in the wearless dry friction process are consistent with the natural frequency of the whole friction system. Conversely, when the excited phonon frequency is far away from the natural frequency of the friction system, the energy dissipation may decrease or be suppressed.
The phonon is the physical particle representing mechanical vibration and is responsible for the transmission of everyday sound and heat [24]. The classical frictional force and power dissipation are essentially due to the excitation of phonons in the dry friction [20]. In general friction processes, such as on rough surfaces, friction excites many phonon modes with different frequencies because rough surfaces are random and aperiodic. But on atomically flat surfaces, only the phonons with washboard frequency and/ or its harmonics may be excited in the wearless dry friction [21,25]. In order to verify the importance of the phonon resonance to dissipate the frictional energy, this work investigated the frictional energy dissipation mechanism caused by inter-layer shear motion of graphene films using molecular dynamics (MD) method. An external periodic excitation is applied to make one layer of graphene vibrate tangentially on the other layer of graphene. The similar method is also used by Sokoloff to explore the possible nearly frictionless sliding for mesoscopic solids [26]. The frequency of the periodic excitation is analogy for the washboard frequency v/a in the wearless dry friction on atomically flat surfaces with sliding velocity v and substrate period a. Thus, only a single phonon mode with excitation frequency dissipates energy in each MD simulation. The resultant energy dissipation due to periodic excitation is also extremely important for the performance of graphene-based MEMS devices [27,28]. Our study shows that two local maxima in the energy dissipation rate of the system does appear at special excitation frequencies. We further explain the simulation results in details by extracting the vibrational density of states and the natural frequency of the system with various system parameters. Figure 1a shows the atomic model for all our MD simulations, including a period-excited upper graphene and a supported lower graphene. The in-plane size is about 90×104 Å 2 and there are 6912 atoms in the simulation system. To illustrate the details of the model more clearly, a schematic diagram of the atomic model is also shown in Fig. 1b. The upper graphene is connected to a support by three independent springs along x, y, and z direction with stiffness k x , k y , and k z , respectively. The k x and k y are the shearing stiffness of the spring and k z is the normal stiffness. By applying periodic excitation to the support end along the x direction, like the washboard excitation in the wearless dry friction, the upper graphene would vibrate relative to the lower graphene and result in energy dissipation. To simplify the simulations, we also used a set of springs along x, y, and z directions with stiffness of k subx , k suby , and k subz to attach to the lower graphene, while the other ends of these springs are fixed. Both k subx and k suby are set as the effective interfacial shearing stiffness k 0 between two graphene layers at the state of equilibrium. The k 0 depends on the used interlayer interactions and can be calculated by differencing the interfacial shearing force with respect to the lateral displacement [21].

Model and Method
In order to distinguish the excited phonons in the friction process from the intrinsic phonons within the graphene at the finite temperature, the simulation temperature is set at 0.001 K. The optimized Tersoff potential [29] is used to describe the carbon-carbon interactions within each graphene layers, while the Lennard-Jones (LJ) potential [30] for the interlayer carbon-carbon interactions between the two graphene layers. Since no normal load is applied in the upper graphene in the MD simulations, the interlayer distance between the two layers of graphene is 3.35 Å. The normal stiffness k z and k subz are set the same as k z = k subz = 18.0 eV/ Å 2 . The shearing stiffness between the two-layer graphene is calculated to be k 0 = 7.47 eV/Å 2 from the used LJ potential. We first equilibrium the simulation system at 0.001 K for 500 ps, and then a periodic excitation with the frequency f 0 and the amplitude Am are applied to the support to imitate the dry friction on an atomically flat surface [21,26]. The Langevin thermostats on both the upper and lower graphene are used to control the temperature of the simulation system. There are no fixed atoms in the simulation system. The periodic boundary conditions are applied along the in-plane directions. The equations of motions are solved with the velocity-Verlet algorithm with a timestep of 0.5 fs.

Results and Discussion
Firstly, we set k x = k y to k 0 , and the amplitude and excitation frequency to Am = 0.1 Å and f 0 = 0.1 THz, respectively. Figure 2 shows the relation between the accumulated energy in the thermostat and the effective simulation time. The results show that the accumulative energy is proportional to the simulation time, indicating the energy dissipation is stable during the periodic excitation. When we increased the excitation frequency, such as to 0.2, 0.4, and 0.6 THz, the cumulative energy in the thermostat still shows a good linear relation with the simulation time, as also shown in Fig. 2. Thus, the energy dissipation rate in the dissipative system under periodic excitation is defined as the proportional coefficient between the accumulated energy and the effective simulation time.
Next, the excitation frequency is gradually increased from 0.1 to 0.6 THz in the MD simulation, and the corresponding energy dissipation rate is calculated as mentioned above. Figure 3c shows that the energy dissipation rate presents two local maxima that corresponds to the excitation frequencies Fig. 1 The atomic model (a) and the corresponding schematic diagram (b) to investigate the energy dissipation of two-layer graphene system under periodic excitations. The stiffness k subx and k suby of the supported springs are set as the effective shearing stiffness k 0 of two-layer graphene, to simulate the lower graphene being supported on the bulk graphite. The stiffness k x and k y of the dragged springs can be adjusted freely based on the effective shearing stiffness k 0 of two-layer graphene in the molecular dynamics simulations. The normal stiffness k z and k subz are set the same as k z = k subz = 18.0 eV/Å 2  In order to confirm our observations, we reduced the shearing stiffness of the spring between the support and the upper graphene layer by two times and four times, respectively. In these cases, the interfacial shearing stiffness between the two graphene layers is greater than that of the spring between support and the upper graphene layer. Figure 3a, b show that two maxima in the energy dissipation rate still appears with the increase of excitation frequency from 0.1 to 0.6 THz. We also increase the shearing stiffness of the spring by a factor of 2, corresponding to the situation on which the interfacial shearing stiffness is less than that of the spring. Figure 3d shows that the energy dissipation rate of the system only has a maximum value with the increasing excitation frequency. However, the energy dissipation rate with the excitation frequency in Fig. 3d can obviously be fitted using two Lorentz functions with center frequency of 0.25 and 0.40 THz, indicating these two special excitation frequencies can result in significant energy dissipation in the friction system.
In order to clearly present the details in Fig. 3, the relation between the maximum energy dissipation rate and the shearing stiffness k x = αk 0 at two characteristic peaks and the relation between the frequency at the maximum energy dissipation and the shearing stiffness are presented in Fig. 4a, b, respectively. Figure 4a shows that the energy dissipation rate at two characteristic peaks increases significantly as the shearing stiffness k x increases. Interestingly, the energy dissipation rate in the first peak at 0.16 THz dominates with the k x of 0.25k 0 , whereas the energy dissipation rate in the second peak at 0.4 THz dominates with the k x of 2k 0 .
To explain the above observation, we first extracted the velocity amplitude of the center of mass of the upper layer graphene along the excitation direction (or x direction). The results in Fig. 5a show that the atomic amplitudes of the two characteristic peaks increase with increasing shearing stiffness. We also extracted the velocity amplitude of the lower graphene in Fig. 5b and found that the velocity amplitude also increases. In atomic-scale friction models, such as the Prandtl-Tomlinson model [2], the energy dissipation rate depends on the atomic velocity amplitude. Generally, the higher the velocity amplitude, the larger the energy dissipation rate will be. Both Fig. 5a, b show that the velocity amplitude increases markedly with increasing shearing stiffness k x for both graphene layers. This is because, when the shearing stiffness k x is small, the shearing force  exerted on graphene under the same excitation amplitude will be small, and the corresponding velocity amplitude and energy dissipation rate will also be small. On the contrary, when the shearing stiffness is increased, the shearing force exerted on graphene will be large, and the corresponding velocity amplitude and energy dissipation rate will also be large. Figure 5a also show that the velocity amplitude of the first characteristic peak is larger than that of the second peak when the shearing stiffness is 0.25k 0 . However, when k x = 2.0k 0 , the velocity amplitude of the first characteristic peak is smaller than that of the second peak. Figure 5b shows that the velocity amplitude of the first peak is slightly larger than that of the second peak. Which peak dominates, the first peak or the second peak, under different shearing stiffness depends on the relative velocity amplitude at the two characteristic peaks. By considering the data in Fig. 5a, b together, it is evident that the velocity amplitude at the first peak dominates at low shearing stiffness k x , while at high shearing stiffness k x the velocity amplitude at the second peak dominates. Figure 4b shows that the frequencies of the two characteristic peaks increase gradually with increasing shearing stiffness k x . For example, when the shearing stiffness increases from 0.25k 0 to 2.0k 0 , the excitation frequency corresponding to the first peak of energy dissipation rate increases from 0.16 to 0.25 THz, and the excitation frequency corresponding to the second peak of the energy dissipation rate also increases from 0.34 to 0.40 THz. This is because the enhancement of the shearing stiffness will increase the effective stiffness of the system assuming that the two layers of graphene in the MD simulations are respectively regarded as two concentrated masses in series. Increasing the effective stiffness of the system certainly increases the resonant frequency. More quantitative explanation is presented with a mass-spring model as followed.
In order to understand the relationship between the energy dissipation rate and the excitation frequency, we first extracted the time-dependent atomic velocities under the excitation for a duration time of 500 ps. Then, these atomic velocities in the lower graphene layer are transformed into the vibrational density of states (vDOS) from the Fourier transform of the velocity autocorrelation function [31]. Figure 6 shows that when the excitation frequency is 0.1 THz, the phonons with frequencies around 0.1 THz are obviously present in the vDOS of the lower layer graphene, while other phonon frequencies are strongly suppressed due to the low simulation temperature. We also separately calculated the vDOS perpendicular to the excitation direction with a duration time of 500 ps. The inset of Fig. 6 shows that the intrinsic phonons within graphene do appear in the vDOS but the intensity is extremely weak. This indicates that the MD simulations at low temperatures are indeed able to suppress the intrinsic phonons within graphene. When changing the excitation frequency, the excited phonons with the same frequency appear in the vDOS of the lower graphene layer. This indicates that large numbers of phonons with the same phonon frequency [THz] vDOS along the z-direction Fig. 6 The vibrational density of states of the lower graphene layer with the excitation frequency of f 0 = 0.1 THz. The high-frequency phonons are suppressed and not presented here due to the low simulation temperature. The inset shows the extremely weak phonon intensity in the z-direction vDOS frequency as the excitation frequency are generated in the lower graphene under external periodic excitation. It is expected that the local maximum of the energy dissipation rate may be resulted from the phonon resonance when the excited phonons in the substrate have the same frequency with the natural frequency of the friction system. Note that the intensity of the phonon mode in the calculated vDOS is almost constant over time in our MD simulations, even on the resonance condition that the excitation frequency equals the natural frequency of the system. This is because the dissipated energy during each excitation period is removed artificially from the simulation system by the thermostat in order to keep the temperature of the simulation system unchanged.
In order to explain the MD results quantitatively, a two-degree of freedom mass-spring model is established as shown in the inset of Fig. 7. By comparing with the MD model, we set both mass m 1 and m 2 in the model as m 1 = m 2 = m = N × m c , where N is the atom number within the upper graphene layer and m c is the mass of each carbon atom. There is no normal load applied in the MD simulation and the energy dissipation mainly depends on the shearing stiffness between the support and the upper graphene layer. Thus, the spring stiffness k in that connects the two mass is set as the effective interfacial shearing stiffness k 0 between the two graphene layers in the state of equilibrium, i.e., k in = k 0 . The spring stiffness k sub that connects the lower mass and the fixed substrate is also set as k 0 , i.e., k sub = k 0 , and the spring stiffness k that connects the upper mass and the vibrating support can be adjusted with the MD model. When a small periodic excitation is applied on the vibration end, the resultant equations of motion for the two-degree of freedom model can be written as where A is a small amplitude, ω = 2πf 0 is the excitation angular frequency, y 1 and y 2 are the corresponding displacements of the lower and the upper degree of freedom, γ 1 and γ 2 are their corresponding damping. By solving the above equation, it is found that the below condition should be satisfied when the amplitude of both degrees of freedom reach the maximum, The detailed process to obtain the above equation can be found in the supporting information. As stated above, the interfacial shearing stiffness between the two graphene layers is equal to the shearing stiffness of the spring that connects the lower graphene layer and the fixed end in the MD simulations. This means that both k in and k sub should be constants and equal to k 0 in the two-degree of freedom mass-spring model, i.e., k in = k sub = k 0 to make the model be consistent with the MD simulations. Thus, we set k sub ∕k = k in ∕k=k 0 ∕k = 1∕ , then the predicted excitation angular frequency that responds to the maximum energy dissipation rate in the friction system should be the two solutions to the Eq. (3), i.e., ωand ω +, , which can be easily obtained. The two special frequencies ωand ω + are the natural angular frequency of the two-degree of freedom, which are shown in Fig. 7 as a function of the shearing stiffness of the spring k = αk 0 . It is shown that the two natural frequencies gradually increase when increasing the shearing stiffness of the dragged spring.
(1) Fig. 7 The natural frequencies of the two-degree of freedom massspring system as a function of the spring stiffness k = αk 0 . The circles on the lines represent four sets of simulation conditions in the MD. The inset is the two-degree of freedom mass-spring model to explain the MD simulations The special frequencies predicted by the two-degree of freedom mass-spring model are further compared with the results of MD simulation. Table 1 shows that, when the controlled shearing stiffness in MD simulations gradually changes from 0.25k 0 to 2k 0 , the excitation frequency corresponding to the maximum energy dissipation rate of the system in Fig. 3 and Fig. 4b is basically consistent with the natural frequency of the system predicted by the two-degree of freedom mass-spring model, which explains our MD simulation results.

Conclusion
In summary, we investigated the frictional energy dissipation in two-layer graphene system under periodic excitation using molecular dynamics simulations. The results show that the energy dissipation rate presents two local maxima when the excitation frequency increases gradually from 0.1 to 0.6 THz. By extracting the vibrational density of states, it is found that the periodic excitation can induce large numbers of phonons in the substrate, which has the same frequency with the excitation frequency. When the frequency of these phonons is equal to the natural frequency of the frictional system, the maximum in the energy dissipation rate appears. The proposed two-degree of freedom massspring model is also confirmed the MD results. Since the excitation frequency in our MD simulation is analogous to the washboard frequency for an atomically flat surface in wearless dry friction, the results indicate that the energy dissipation rate in the frictional system can be modulated by controlling the washboard frequency and the resonant frequency of the system. This may explain the local maximum of the friction force with the sliding velocity [21]. As for the commensurability-dependent friction behavior on an atomically flat surface, the phonon resonance can occur at the commensurate interface because the excited phonon modes on the two sides of the interface are the same due to the same period along the sliding direction. But at the incommensurate interface, the atomic periods on the two sides of the interface are different along the sliding direction, which cannot result in phonon resonance. Thus, the phonon resonance mechanism can also explain the larger friction at commensurate interfaces than that at incommensurate interfaces [32]. Our study establishes the relationship between the excitation frequency and energy dissipation rate in a model two-layer graphene system. The uncovered phonon resonance mechanism can be used to regulate the friction force and energy dissipation in many systems, including the nano-sensors and actuators.