Van der Waals interaction-tuned heat transfer in nanostructures

Abstract

Interfaces usually impede heat transfer in heterogeneous structures. Recent experiments show that van der Waals (vdW) interactions can significantly enhance thermal conductivity parallel to the interface of a bundle of nanoribbons compared to a single layer of freestanding nanoribbon. In this paper, by simulating heat transfer in nanostructures based on a model of nonlinear one-dimensional lattices interacting via van der Waals interactions, we show that the vdW interface interaction can adjust the thermal conductivity parallel to the interface. The efficiency of the adjustment depends on the intensity of interactions and temperature. The nonlinear dependence of the conductivity on the intensity of interactions agrees well with experimental results for carbon nanotube bundles, multi-walled carbon nanotubes, multi-layer graphene, and nanoribbons.

---

1. Introduction

Harnessing heat transfer at the nanoscale is essential for the development of microelectronic, photovoltaic, thermoelectric^1–4 and information processing devices.^5,6 Recent years have witnessed important progress in the utilization of heat transfer on the nanoscale. For example, a thermal diode that rectifies heat current is proposed in nonlinear lattices^7–10 and experimentally realized with gradual mass-loaded carbon and boron nitride nanotubes.¹¹ Models of thermal logic gates⁵ and even thermal memories⁶ which provide the possibility to design a “phononic computer” are proposed as well. All of this progress is based on control of heat transfer in multi-component nanostructures. However, heat transfer in non-metallic and semiconductor nanomaterials is dramatically affected by scattering of phonons at surfaces and interfaces.¹² It is also strongly adjusted by the interactions between different components. The van der Waals (vdW) interaction,^13–16 ubiquitous in nanostructures, may dominate heat transfer^17–23 in a variety of configurations such as multi-walled carbon nanotubes, multi-layer graphene, interacting nanoribbons, and composite structures of these elements (Fig. 1a–c). While most work was focused on understanding heat transfer perpendicular to interfaces in heterogeneous materials, which lead to a reduction in heat transfer, recent experiments of Yang et al.²⁴ demonstrate that the longitudinal thermal conductivity, parallel to the interface of bonding boron nanoribbons (Fig. 1c), can change from a significantly higher value to the value of a single freestanding nanoribbon by adjusting the interface interaction. Yang et al.²⁴ explain the enhancement of the conductivity of the interacting nanoribbons using a phonon transmission model. On the other hand, Ghosh et al.²⁵ experimentally show that the in-plane thermal conductivity of multi-layer graphene (MLG) (Fig. 1b) is lower than that of a single-layer graphene (SLG) and the conductivity further decreases as the number of layers increases. However, so far, to our best knowledge, a unified picture for the role of vdW interaction in heat transfer for a wide spectrum of nanostructures is yet to be revealed. We show in this paper that the differences between the thermal conductivities of single-walled carbon nanotubes (CNT), multi-walled CNTs (MWCNTs) as well as CNT-bundles (Fig. 1a) actually suggest the important role of the vdW interface interactions in heat transfer in these nanostructures.^26–28 The interface vdW interaction, which glues multi-component nanostructures together, can be used to engineer the contact resistance on the nanoscale to control heat transfer, and further to create nonlinear behavior such as rectification and switching.²⁹

Fig. 1 Schematics of nanotube bundles and MWCNTs (a), multi-layer graphene (b), interacting nanoribbons (c), and a phononic model (d).

2. Simulation methodology

We simulate heat transfer along the longitudinal direction of nonlinear parallel lattices, bonded by vdW interface interaction, with a phononic model illustrated in Fig. 1d. First, we study two identical chains of nonlinear lattices^3,30–32 bundled by vdW interaction (Fig. 2). Each chain is modeled by a 1D Frenkel–Kontorova (FK) lattice³³ which is widely used in condensed-matter physics and nonlinear physics,^{5–9,34–39} whose Hamiltonian readswhere x̃_i and p̃_i denote the position and momentum of the particles. The motion of all the particles is restricted along the chain, and the influence of other atoms belonging to the same atomic layer but not on the chain is summarized by the cosine on-site potential,⁴⁰ which is in agreement with the phonon dynamics model.²⁵m̃, k̃, μ, Ṽ, λ, and N represent the mass, the harmonic coupling constant, the equilibrium distance of the interatomic potential, the strength of the on-site potential, the period of the on-site potential, and the system length, respectively. Introducing non-dimensional parameters m = m̃/m₀, k = k̃/(m₀ω²₀) and x_i = x̃_i/λ, we obtain the conventional form of the FK model,where H_C = H̃_C/(m₀λ²ω²₀), p_i = p̃_i/(m₀λω₀), V = Ṽ/(m₀λ²ω²₀) and a = μ/λ. m₀ and ω₀ are the mass and oscillating frequency of typical atoms with m₀ ∼ 10⁻²⁶ to 10⁻²⁷ kg and ω₀ ∼ 10¹³ s⁻¹.^3,36,39λ is set to be the carbon–carbon bond length λ = 0.142 nm. In our simulations, we set m = a = 1, k = 0.5 and V = 2.5. The Hamiltonian of the system iswhere M stands for the number of chains in a system with more than two chains, and the superscripts α and β stand for the α^th chain and β^th interface, respectively. H_vdW represents the vdW interface, and we adopt the Lennard-Jones (LJ) potential to describe the chain–chain nonbonded interactions:¹⁶where r_i is the distance between the i^th particle pair on the two sides of the vdW interface. The LJ potential is normalized by m₀λ²ω²₀ so that the dimensionless energy parameter ε = [small epsilon, Greek, tilde]/m₀λ²ω²₀ denotes the intensity of the vdW interaction. σ denotes the distance parameter and is set to be σ = 2.5. The vertical distance d = 2.359 between the two chains of the vdW interface is the equilibrium distance of two chains interacting with LJ potential at T₀ = 0.2. This corresponds to a real distance d_r = 0.335 nm, the interlayer distance of graphite.¹⁸ Usually, d is slightly smaller than σ, as reported in previous studies.⁴¹ As a result, the particles on the two sides of the vdW interface are subjected to repulsive and attractive forces during the motion.

Fig. 2 A two-chain system (M = 2). The upper and lower chains are a FK model and bonded by vdW interaction. The four particles on the left and right ends are in contact with two heat baths at temperatures T_L and T_R, respectively.

Fixed boundary conditions are applied, as shown in Fig. 2, and the two ends of the chains are coupled with two Langevin heat baths^30,36 at temperatures T_L = T₀(1 + Δ) and T_R = T₀(1 − Δ), where T₀ denotes the average temperature of the system and Δ is the dimensionless temperature difference. We have checked that the thermal conductivity of the FK model is finite and independent of the system length when N > 400, which is in agreement with previous calculations.^36,38,39 Unless otherwise stated, we fix Δ = 0.4 and N = 500, so that the main adjustable parameters are only T₀, ε and M. We integrate the equations of motion using the velocity Verlet algorithm described in ref. 42 and the time step is set to δt = 0.005. We have checked that the Hamiltonian of the system is conserved throughout the computations, which is in contrast to the Hamiltonian dissipation using other algorithms, e.g., the Runge–Kutta algorithm. The local temperature for a particle is defined as T_i = m〈x²_i〉, where 〈〉 stands for time average. The local heat flux, as defined through the continuity equation,³⁵ is modified as

Simulations are performed long enough to allow the system to reach a non-equilibrium steady state in which the local heat flux is constant along the chains. To obtain a steady state, the total integration is more than 10⁹ time steps in all simulations. In the final state, the time-averaged heat flux J_i reaches a constant value independent of the site i, which suggests that the system reaches a steady state. Accordingly, the thermal conductivity for a finite system can be evaluated as

where S₀ is the cross section of a chain.

3. Results and discussion

Fig. 3a shows the ratio of the thermal conductivity (κ) of the chain bundle of two chains (M = 2) to that of a single chain (κ_S at T₀ = 0.2) with respect to the average temperature (T₀) for different normalized energy parameters ε/ka² = 0.0, 1.0, 2.0, 3.0, 4.0, 5.0 and 6.0 with fixed Δ = 0.4. Fig. 3b depicts the tuning factor of thermal conductivity defined as α = ((κ − κ_S)/κ_S) × 100%. The factor decreases with the temperature for all the considered values of ε/ka². In order to compare our simulations with experiments, it is necessary to establish the relation between the real temperature and the parameters used in the simulations. Based on the analysis for a typical atomic lattice in the previous studies,^3,7,36,39 the real temperature T_r is related to the dimensionless temperature T by T_r = τT, where the proportionality factor τ and the energy unit m₀λ²ω²₀ in our simulations are ∼1500 K and ∼126 meV, respectively; therefore the room temperature T_r = 300 K and the energy magnitude 126 meV approximately correspond to a dimensionless temperature T = 0.2 and a normalized energy parameter ε/ka² = 2.0, respectively. As shown in Fig. 3(a), κ is higher than κ_S in the low temperature region for weak vdW interface interactions, and this is also true in the whole considered temperature region for strong vdW interface interactions (e.g. ε/ka² = 6.0). In the experiment,²⁴ the adhesion energy in the interfaces of the nanoribbon bundles is estimated as 320 mJ m⁻². It corresponds to a normalized energy parameter ε/ka² = 4.26 in our simulations. Fig. 3 shows that κ > κ_S when ε/ka² ∈ (4.0–5.0) and T₀ < 0.26, which corresponds to T_r < 390 K. Thus the simulations in Fig. 3 are in agreement with the experimental results.²⁴

Fig. 3 (a) Comparison of κ and κ_Sversus different values of ε/ka² and T₀ for dimensionless temperature difference Δ = 0.4. (b) Dependence of tuning factor on T₀versus different values of ε/ka².

On the other hand, the adhesion energy between two graphene sheets or collapsed CNTs is ∼40 meV atom⁻¹,⁴³ equivalent to about 200 mJ m⁻², which corresponds to ε/ka² = 2.66 ∈ (2.0–3.0). Fig. 3a shows that κ < κ_S in this range of ε/ka² for T₀ near or above 0.2. The simulation is in agreement with the recent experimental results for graphene.²⁵ However, even for weak vdW interactions (e.g. ε/ka² = 2.0), κ can still be higher than κ_S at low temperatures; i.e., T₀ < 0.14. These findings show that the vdW interface can adjust thermal conduction in a wide range of system parameters.

To understand the underlying mechanism of the dependence of the thermal conductivity on T₀, let us start from the temperature dependence of the particles motion.³⁴ For the FK model, the particles are confined in the valleys of the on-site potential at very low temperature since the kinetic energy from the heat bath is low, and the particles are in the localized status. The role of the vdW interface is like a channel in the system when the heat flows parallel to the interface and permits the phonons to transmit through the interface with weak scattering. Therefore, the vdW interface interaction leads to an enhancement in thermal conduction. However, the situation is strikingly different in the high temperature region where the particles are in the thermal fluctuation status. Fig. 3a shows that κ < κ_S in this temperature region if the intensity of the vdW interface interaction is low (e.g. ε/ka² = 2.0), and the difference becomes larger as T₀ further increases. It strongly suggests that the vdW interface no longer contributes to κ; instead, it plays a role like a barrier that strongly scatters phonons since the particles lie in the thermal fluctuation status and leads to a reduction of thermal conductivity; i.e., κ < κ_S.

Fig. 4 shows the dependence of κ on ε/ka² for T₀ = 0.2. κ is almost unchanged until ε/ka² = 0.03 since an overly weak vdW interface interaction has no effect on the thermal conduction. There is a depression area, i.e., κ < κ_S when 0.03 < ε/ka² < 3.0. The adhesion energy inside the MLG or collapsed CNTs is in the range 4–210 mJ m⁻² from the previous measurements and calculations;^18,43–45 this corresponds to 0.053 < ε/ka² < 2.8. Previous theoretical studies show that the thermal conductivity of the CNT-bundles is lower than that of the SWCNTs.²² Experimentally, Pettes et al.²⁷ show that the thermal conductivity of MWCNTs is lower than that of the SWCNTs and further decreases as the number of shells increases; and Ghosh et al.²⁵ show that the thermal conductivity of the MLG is lower than that of the SLG. These results are all in agreement with our simulations as shown in Fig. 4, i.e., κ < κ_S. As ε/ka² further increases into the range 3.0 < ε/ka² < 40, κ increases sharply with ε/ka². The adhesion energy of the nanoribbons in the experiment of Yang et al.,²⁴ showing that strong vdW interface interactions increase the thermal conductivity of nanoribbons, evaluated as 320 mJ m⁻² (ref. 24) and stronger than those in the MLG and CNTs, corresponds to ε/ka² = 4.26, i.e., ε/ka² > 3.0. Our simulation also agrees with the experiment.

Fig. 4 Dependence of normalized thermal conductivity κ of chain bundle on normalized energy parameter ε/ka² for M = 2 and T₀ = 0.2 which corresponds to room temperature T_r = 300 K.

The dependence of heat transfer on ε/ka² can be analyzed from the power spectra.^3,7,8 The power spectra, which describes the distribution of a particle's kinetic energy in the frequency domain, is calculated as: , where f and u are the frequency and velocity of the particle, and t and t_final are the time and the final time of the calculation. Based on the Parseval's theorem one can obtain , where 〈E_k〉 is the average kinetic energy of the particle.

Fig. 5 shows that the power spectra of the mid-particle of the upper chain for ε/ka² = 0 and ε/ka² = 0.002 almost overlap, which implies κ/κ_S = 1, in agreement with the trend in Fig. 4. However, the results are strikingly different in the moderate domain (e.g. ε/ka² = 1.0): the power spectra in the moderate-frequency domain become more dispersed and shift towards the high-frequency region. The thermal conductivity can be estimated from κ = Cvl/3, where C is the specific heat, v is the average phonon velocity and l is the mean free path (MFP) of the phonon. The phonon in the high-frequency region, which means shorter MFP, will lead to a decreasing thermal conductivity, i.e., κ/κ_S < 1.

Fig. 5 Power spectra of the mid-particle of the upper chain for different normalized energy parameter ε/ka² with M = 2 and T₀ = 0.2.

However, the power spectrum splits into two branches when ε/ka² further increases (e.g. ε/ka² = 10.0). The branch of the low-frequency keeps the major features of ε/ka² = 0 while a branch of high-frequency is excited. Phonons with a high frequency indicate that the MFP is short and this will lead to a negative effect on the thermal conductivity. On the other hand, the branch of the high-frequency will lead to an enhanced energy exchange along the chain since the on-site potential can be easily overcome by the phonon with the high frequency. Moreover, when ε/ka² further increases to the extreme case of a welded contact, the phonon scattering at the interface will vanish. Therefore, the excited branch of the high-frequency due to the high intensity of the vdW interaction leads to an increase in the thermal conductivity, i.e., κ/κ_S > 1.

In order to reveal the size effect on heat transfer, we consider a system containing multiple chains bundled by the vdW interactions. The dependence of κ on M for T₀ = 0.2 and ε/ka² = 1.6 is plotted in Fig. 6. It clearly shows that κ decreases with the number M of layers, exhibiting a significant size effect. These results are in qualitative agreement with the recent experiment of the thermal conduction for MLGs²⁵ and MWCNTs.²⁷ Moreover, as shown in Fig. 6, κ tends to an asymptotic value as M increases. In particular, κ after M ≈ 5 almost remains constant. This observation agrees well with the theory for lattice thermal conductivities which is based on the Boltzmann equation for phonons and predicts that the conductivity of MLG of more than five layers approaches that of graphite.⁴⁶ In addition, the asymptotic value obtained from the present phononic model falls into the range given by Ghosh et al.²⁵

Fig. 6 Dependence of thermal conductivity κ on number M of chains with T₀ = 0.2 and ε/ka² = 1.6, which corresponds to adhesion energy 120 mJ m⁻², used here to simulate the adhesion energy inside the MLGs and MWCNTs. The green squares and blue triangles are from the experimental measurements and the first-principles theory, respectively.²⁵ The navy triangles are calculations based on the Boltzmann equation.⁴⁶ The magenta and orange triangles are obtained from the molecular dynamics (MD) simulations for the zigzag and armchair graphene, respectively.⁴⁷ The dotted straight lines denote the range of thermal conductivities of bulk graphite.²⁵

To illustrate the underlying mechanism of size effect, we present the power spectra of the mid-particle of the uppermost chain for different numbers of chains in Fig. 7. The power spectra markedly expands to the moderate-frequency domain when M = 2, which means a short MFP. Furthermore, it has been demonstrated that the low frequency phonon modes dominate heat transfer in graphene.^18,46 Therefore the thermal conductivity will decrease. As M further increases, the power spectra in the moderate-frequency domain increases, although not significantly. This is consistent with the findings shown in Fig. 6 that the decrease of the thermal conductivity becomes slow when M > 2. We have checked that the power spectra remain almost unchanged when M > 5 and this is again in agreement with the calculations that κ almost remains constant after M ≈ 5.

Fig. 7 Power spectra of the mid-particle of the uppermost chain for different numbers of chains with ε/ka² = 1.6 and T₀ = 0.2.

4. Conclusions

In summary, we have illustrated a unified picture of the important role of the van der Waals interaction-adjusted heat transfer in nanostructures using a simple lattice model. We have demonstrated that the vdW interaction can either increase or decrease the thermal conductivity in nanostructures parallel to the interface depending on the intensity of interaction and temperature. The conductivity exhibits a nonlinear and non-monotonic dependence on the intensity of the vdW interaction. The simulations are in agreement with the experimental results of the thermal conductivity of carbon nanotube bundles, multi-walled carbon nanotubes, multi-layer graphene, and nanoribbons, which have different intensities of interaction. Moreover, the size effects predicted by the simulations for vdW interacted structures also agree with the previous predictions. More importantly, the model and the simulations, which reveal the effect of the vdW interaction on heat transfer in a wide spectrum of nanostructures, have applications in tuning heat transfer in nanostructures through a combined size effect and interface effect, which are necessary for the development of nanoscale thermal devices and for the thermal management of micro- and nano-devices.

Acknowledgements

This work is supported by the National Natural Science Foundation of China (Grant no. 10932001), and Major State Basic Research Development Program of China (Grant no. 2011CB013101). The authors thank Dr H. Wang for helpful discussions and suggestions.

References

1. A. I. Hochbaum, et al., Enhanced thermoelectric performance of rough silicon nanowires, Nature, 2008, 451(7175), 163–167.
2. M. D. Losego, et al., Effects of chemical bonding on heat transport across interfaces, Nat. Mater., 2012, 11(6), 502–506.
3. N. Li, et al., Colloquium: phononics: manipulating heat flow with electronic analogs and beyond, Rev. Mod. Phys., 2012, 84(3), 1045–1066.
4. C. Yan, J. H. Cho and J.-H. Ahn, Graphene-based flexible and stretchable thin film transistors, Nanoscale, 2012, 4(16), 4870.
5. L. Wang and B. Li, Thermal logic gates: computation with phonons, Phys. Rev. Lett., 2007, 99(17), 117208.
6. L. Wang and B. Li, Thermal memory: a storage of phononic information, Phys. Rev. Lett., 2008, 101(26), 267203.
7. B. Li, L. Wang and G. Casati, Thermal diode: rectification of heat flux, Phys. Rev. Lett., 2004, 93(18), 184301.
8. B. Li, J. Lan and L. Wang, Interface thermal resistance between dissimilar anharmonic lattices, Phys. Rev. Lett., 2005, 95(10), 104302.
9. B. Hu, L. Yang and Y. Zhang, Asymmetric heat conduction in nonlinear lattices, Phys. Rev. Lett., 2006, 97(12), 124302.
10. M. Terraneo, M. Peyrard and G. Casati, Controlling the energy flow in nonlinear lattices: a model for a thermal rectifier, Phys. Rev. Lett., 2002, 88(9), 094302.
11. C. W. Chang, et al., Solid-state thermal rectifier, Science, 2006, 314(5802), 1121–1124.
12. V. Varshney, et al., A novel nano-configuration for thermoelectrics: helicity induced thermal conductivity reduction in nanowires, Nanoscale, 2012, 4(16), 5009.
13. R. Prasher, Acoustic mismatch model for thermal contact resistance of van der Waals contacts, Appl. Phys. Lett., 2009, 94(4), 041905.
14. L. Zhao, et al., Plasmon-induced modulation of the emission spectra of the fluorescent molecules near gold nanorods, Nanoscale, 2011, 3(9), 3849.
15. G. Gao, et al., Artificially stacked atomic layers: toward new van der Waals solids, Nano Lett., 2012, 12(7), 3518–3525.
16. H. Zhong and J. Lukes, Interfacial thermal resistance between carbon nanotubes: molecular dynamics simulations and analytical thermal modeling, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 74(12), 125403.
17. V. Varshney, et al., Molecular dynamics simulations of thermal transport in porous nanotube network structures, Nanoscale, 2011, 3(9), 3679.
18. J. H. Seol, et al., Two-dimensional phonon transport in supported graphene, Science, 2010, 328(5975), 213–216.
19. G. G. Yadav, et al., Self-templated synthesis and thermal conductivity investigation for ultrathin perovskite oxide nanowires, Nanoscale, 2011, 3(10), 4078.
20. Z. Xu and M. J. Buehler, Nanoengineering heat transfer performance at carbon nanotube interfaces, ACS Nano, 2009, 3(9), 2767–2775.
21. D. Donadio and G. Galli, Thermal conductivity of isolated and interacting carbon nanotubes: comparing results from molecular dynamics and the Boltzmann transport equation, Phys. Rev. Lett., 2007, 99(25), 255502.
22. J. Hone, et al., Thermal conductivity of single-walled carbon nanotubes, Phys. Rev. B: Condens. Matter, 1999, 59, R2514–R2516.
23. P. Sheridan, et al., Device and SPICE modeling of RRAM devices, Nanoscale, 2011, 3(9), 3833.
24. J. Yang, et al., Enhanced and switchable nanoscale thermal conduction due to van der Waals interfaces, Nat. Nanotechnol., 2012, 7(2), 91–95.
25. S. Ghosh, et al., Dimensional crossover of thermal transport in few-layer graphene, Nat. Mater., 2010, 9(7), 555–558.
26. A. A. Balandin, Thermal properties of graphene and nanostructured carbon materials, Nat. Mater., 2011, 10(8), 569–581.
27. M. T. Pettes and L. Shi, Thermal and structural characterizations of individual single-, double-, and multi-walled carbon nanotubes, Adv. Funct. Mater., 2009, 19(24), 3918–3925.
28. M. S. Dresselhaus and P. C. Eklund, Phonons in carbon nanotubes, Adv. Phys., 2000, 49(6), 705–814.
29. C. Dames, Thermal materials: pulling together to control heat flow, Nat. Nanotechnol., 2012, 7(2), 82–83.
30. S. Lepri, Thermal conduction in classical low-dimensional lattices, Phys. Rep., 2003, 377(1), 1–80.
31. A. Dhar, Heat transport in low-dimensional systems, Adv. Phys., 2008, 57(5), 457–537.
32. R. Livi and S. Lepri, Thermal physics – heat in one dimension, Nature, 2003, 421(6921), 327–327.
33. O. M. Braun and Y. S. Kivshar, Nonlinear dynamics of the Frenkel–Kontorova model, Phys. Rep., 1998, 306(1–2), 1–108.
34. J. Wang and Z. G. Zheng, Heat conduction and reversed thermal diode: the interface effect, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81(1), 011114.
35. E. Díaz, R. Gutierrez and G. Cuniberti, Heat transport and thermal rectification in molecular junctions: a minimal model approach, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84(14), 144302.
36. B.-q. Ai and B. Hu, Heat conduction in deformable Frenkel–Kontorova lattices: thermal conductivity and negative differential thermal resistance, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83(1), 011131.
37. W.-R. Zhong, Different thermal conductance of the inter- and intrachain interactions in a double-stranded molecular structure, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81(6), 061131.
38. B. Hu and L. Yang, Heat conduction in the Frenkel–Kontorova model, Chaos, 2005, 15(1), 015119.
39. B. B. Hu, B. W. Li and H. Zhao, Heat conduction in one-dimensional chains, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57(3), 2992–2995.
40. W. Atkinson and N. Cabrera, Motion of a Frenkel–Kontorowa dislocation in a 1-dimensional crystal, Phys. Rev., 1965, 138(3A), A763–A766.
41. Y. Shibuta and J. A. Elliott, Interaction between two graphene sheets with a turbostratic orientational relationship, Chem. Phys. Lett., 2011, 512(4–6), 146–150.
42. Computer Simulation of Liquids, ed. M. P. Allen and D. J. Tildesley, Oxford University Press, Clarendon, 1989, p. 408.
43. L. X. Benedict, et al., Microscopic determination of the interlayer binding energy in graphite, Chem. Phys. Lett., 1998, 286(5–6), 490–496.
44. R. Huang, Graphene: show of adhesive strength, Nat. Nanotechnol., 2011, 6(9), 537–538.
45. C. Lee, et al., Measurement of the elastic properties and intrinsic strength of monolayer graphene, Science, 2008, 321(5887), 385–388.
46. L. Lindsay, D. Broido and N. Mingo, Flexural phonons and thermal transport in multilayer graphene and graphite, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83(23), 235428.
47. W.-R. Zhong, et al., Chirality and thickness-dependent thermal conductivity of few-layer graphene: a molecular dynamics study, Appl. Phys. Lett., 2011, 98(11), 113107.

---