A comprehensive review on the molecular dynamics simulation of the novel thermal properties of graphene

Abstract

This review summarizes state-of-the-art progress in the molecular dynamics (MD) simulation of the novel thermal properties of graphene. The novel thermal properties of graphene, which include anisotropic thermal conductivity, decoupled phonon thermal transport, thermal rectification and tunable interfacial thermal conductance, have attracted enormous interest in the development of next-generation nano-devices. Molecular dynamics simulation is one of the main approaches in numerical simulation of the novel thermal properties of graphene. In this paper, the widely used potentials of MD for modeling the novel thermal properties of graphene are described first. Then MD simulations of anisotropic thermal conductivity, decoupled phonon thermal transport, thermal rectification and tunable interfacial thermal conductance are discussed. Finally, the paper concludes with highlights on both the current status and future directions of the MD simulation of the novel thermal properties of graphene.

---

1. Introduction

In recent years, graphene, a single layer of carbon atoms densely packed in a sp² bonded honeycomb lattice, has attracted remarkable attention all over the world. The strong and anisotropic sp² bonding and low-mass carbon atoms in the microscopic structure give graphene exceptional physiochemical characteristics compared with traditional carbon- and silicon-based materials. These extraordinary properties, e.g., good deformation beyond the linear regime,¹ room temperature superconductivity,² ballistic electronic propagation,³ realization of the Klein paradox,⁴ and metal free magnetism,⁵ etc., have made graphene a potential candidate for the next generation nanoelectronics.

Compared to other unique features of graphene, its novel thermal properties have gained the most interest because of their direct and widespread applications in both academy and industry.^6–10 The novel thermal properties of graphene mainly refer to its anisotropic thermal conductivity, decoupled phonon thermal transport, thermal rectification and tunable interfacial thermal conductance. These four novel thermal properties have found extensive usage in practical applications and distinguish graphene from other traditional materials explicitly. Thus, it is necessary to discuss the current status of the study of the novel thermal properties of graphene to help better utilize their superiorities in future applications.

The novel thermal properties of graphene have been extensively studied by both experiment and numerical simulation. Although experiment can provide direct measurements of related novel thermal properties, the high experimental costs and limited capability of micro/nanoscale temperature probing and thermal detection have posed great challenges to further comprehensive investigations. Numerical simulations, where established theories are employed to describe the fundamental interactions between atoms and complex phenomena are solved through conservation equations on computers, is a promising tool to complement experiments to advance our understanding of the novel thermal properties possessed by graphene. Therefore, recent years have seen a dramatic increase in numerical simulations applied to the investigation of the novel thermal properties of graphene.

So far, the main numerical approaches applied to the investigation of the novel thermal properties of graphene are first-principle quantum mechanics (FPQM),^11,12 the Boltzmann transport equation (BTE),^13,14 and molecular dynamics (MD).^15–19 Though FPQM describes the graphene system most accurately, its modeling ability is highly restricted by its spatiotemporal scale. The BTE method has less statistical noise and can tune various parameters simultaneously; however, it is subject to great numerical instabilities. Classical MD simulations, in which the system is driven by the FPQM-derived interatomic energy potential, provides a reasonable balance between system size, modeling accuracy, and computational speed. Therefore, so far, a great volume of numerical work on the novel thermal properties of graphene has employed MD and considerable new phenomena have been revealed.

This review summarizes comprehensively the numerical efforts devoted to MD simulation of the novel thermal properties of graphene and its derived materials, such as graphene nanoribbons (GNR). Without further clarification, the studies in the following discussion were all performed by MD. The paper is organized as follows. In the next section, the basic theory of MD will be described and discussed. Then the following four sections will summarize relevant MD work on each novel thermal property. In these sections, the discussed novel thermal property is first defined and then a state-of-the-art discussion of the progress is presented. Finally, the paper is concluded and future research directions are highlighted.

2. The principle of molecular dynamics for graphene

Molecular dynamics solves the motions of Lagrangian particles (atoms, molecules, granules, etc.) through the classical Newton’s second law. It is by far the most efficient and widely-used numerical approach for studying thermodynamic properties in atomic systems. In MD, atoms with initial positions and velocities are exposed to collisions governed by the empirical interatomic potentials (EIP). Principally, the force acting on each atom can be obtained from the summation of the spatial derivative of interatomic energy potentials. Thus, it can be concluded that the evolution of a system is ultimately determined by the type of interatomic energy potential. With the calculated force, the position and velocity of each atom are updated within each time step. The system behavior is concisely described by the evolution of all the atoms and desired properties can be easily extracted. Because of its inherent simplicity in mathematical formulation and numerical algorithm, MD has been widely viewed as an “indispensable bridge” to connect the microscopic atomic structures to macroscopic material characteristics.

Specific to the MD simulations in graphene, several types of EIPs are widely used to describe the C–C interactions. The Tersoff potential²⁰ is the original and popular potential used to model graphene, which can be expressed as

where E is the three-body system energy, and f_C(r_ij) is the cutoff function to limit the range of the potential since the short-ranged functions permit a tremendous reduction in computational effort. The function f_R(r_ij) represents the repulsive pair potential, which includes the orthogonalization energy when atomic wave functions overlap, and f_A(r_ij) is an attractive pair potential associated with bonding. The b_ij term is a monotonically decreasing function of the coordination of atoms i and j. The Tersoff EIP provides an effective approach for the modelling of the graphene system. However, since its original parameters were optimized to describe the bulk silicon and carbon systems, its accuracy in describing 2D membrane systems remains controversial.

The Brenner EIP,²¹ which is another option to describe the solid-carbon structures, shares the same control functions as those in original the Tersoff EIP. The main differences between the Tersoff and Brenner EIPs are that the Brenner EIP includes two additional exponential terms with corresponding adjustable parameters in the attractive pairwise term f_R(r_ij). It includes a screened Coulomb term in the repulsive pairwise term and uses a fifth-order polynomial spline between the bond orders for diamond and graphite. It also contains a dihedral bending term for bond energies which plays a role in carbon nanotubes (CNT) and graphene.

The second generation of Brenner EIP,^21–23 which is the reactive empirical bond-order (REBO) potential based on the Tersoff EIP, is another widely used potential for graphene. The REBO EIP is chosen because its functions and parameters are known to give reasonable predictions for the thermal properties of graphene, whereas the adaptive intermolecular reactive empirical bond order (AIREBO) was reported to underestimate the dispersion of out-of-plane acoustic (ZA) phonons in graphene. The AIREBO potential is expressed as

When the LJ and TORSION terms are turned off, it becomes the second generation of the REBO potential. The original parameter sets of Tersoff and Brenner EIP cannot reproduce the phonon dispersions of graphene accurately, as the velocities of the three acoustic branches near the center of the Brillouin zone are misrepresented.²⁴ A set of optimized potential parameters for Tersoff and Brenner EIPs are calculated by Lindsay et al.²⁵ to better describe the C–C interactions in graphene and carbon-nanotube structures.

3. Anisotropic thermal conductivity

Graphene is inherently thermally anisotropic.^26–28 This is mainly due to the different phonon couplings along a given crystal axis, which can exhibit diverse thermal conductivities along different crystal directions. Phonon thermal transport in graphene is highly dependent on the chiral angles and structural formations. For an isotropic material, its thermal conductivity (κ) can be defined as follows

q′′ = −κ·ΔT, (3)

where q′′ is the heat flux and ΔT is temperature gradient. However, the conductivity of graphene should be a second order tensor with nine components as

Phonon scattering at GNR’s rough edges is partially diffusive and can be accurately described by a momentum-dependent specularity parameter that represents the probability of certain phonon modes being scattered from the edges. This angular variation becomes more significant with smaller GNR widths, which is due to increased edge–surface ratio and stronger phonon boundary scatterings. Due to the D_6h (antiprismatic point group) symmetry in the honeycomb structure in graphene, group velocities of the phonon modes are direction dependent. The calculated κ depends weakly on the direction of the thermal flux periodically with a period of π/3.²⁹ Using a non-equilibrium molecular dynamics (NEMD) approach, the thermal conductivities of the armchair GNR (a-GNR) and zigzag GNR (z-GNR) can be characterized. It has been calculated that the thermal conductivity of graphene in the zigzag direction is much higher than that in the armchair direction.³⁰ Recently, it has been argued that the rough edge scattering theory isn’t sufficient to explain the huge κ differences between zigzag and armchair edges since the dominant phonon wavelength can be orders of magnitude longer than the edge roughness profiles.³¹ By decomposing the heat flux along the width direction of GNR, strong suppression of thermal transport is observed in the edge areas, which is due to the localization of phonon populations. The results on anisotropic thermal conductivity of armchair and zigzag GNRs with temperature and system width are shown in Fig. 1. It is found that the high participation rate of the phonon mode can enhance the phonon scattering, especially armchair chiral, resulting in the reduction of GNR’s thermal conductivity at the boundaries. By edge passivation and isotope engineering, the in-plane thermal conductivity of graphene can be tuned and the anisotropic effect can be either suppressed or magnified.

Fig. 1 Anisotropic thermal conductivity of armchair and zigzag GNRs with temperature and system width; the length of the graphene system is 15 nm for all cases (with permission from AIP Publishing LLC for use in this paper).³¹

The intrinsic ripples in graphene exhibit inharmonic couplings between bending and stretching modes, which can suppress the collapse of long-range order in 2D crystals and avoid crumpled structures. When the ripples grow into large-magnitude wrinkles, the mean bond length in graphene is elongated and the phonon density of states increases significantly. Due to the chirality dependent bending stiffness and wave velocity difference, transverse wave propagation in GNR was found to be highly sensitive to the excitation frequency of the vibration source and chirality directions.³² The minimum permissible wavelength was found to be 3.69 Å in the zigzag direction while in the armchair direction the counterpart is 5.68 Å. Therefore, by adjusting the wrinkle level and direction into certain textures, the κ of graphene along different chiral directions can be tuned. Demonstration of one possible textural configuration is shown in Fig. 2(a). By applying a vertical shear strain on the left boundary of GNR, clear and well-defined out-of-plane wrinkles are created. Temperature profiles from the NEMD simulation are shown on the left y-axis. Atom displacements in the out-of-plane direction are shown on the right y-axis. It is discovered that by carefully designing the out-of-plane wrinkle textures, κ of GNR can be reduced by 80% compared with unstrained structures, as is shown in Fig. 2(b).³³ It can be observed that thermal conductivity along the textual direction is severely suppressed only subtly changes parallel to the wrinkling direction. Both the configuration and strain of wrinkles determine the κ of GNR. Therefore, the anisotropic thermal conductivity can be either suppressed or enhanced with proper wrinkle designs.

Fig. 2 (a) Wrinkling cross-section configuration, temperature gradient along the texture direction and wrinkling model description; (b) relative thermal conductivity along the texture, and wrinkling directions under different wrinkling levels (with permission from Royal Society of Chemistry for use in this paper).³³

The extent of anisotropic thermal conductivity can also be affected by structural defects such as non-hexagonal rings and vacancies. At low defect concentrations, the defect areas will behave as scattering centers which induce high local phonon scatterings. In addition, at high defect concentrations, delocalized scattering will further suppress the thermal conductivity due to the interactions between scattering centers. The effects of these defects on a-GNR might have been somewhat less due to the presence of edge effects, which already posed significant phonon scattering. By introducing Stone–Thrower–Wales (STW) defects in the heat flux direction of GNR, the thermal conductivity of z-GNR is reduced by 54% and a-GNR by 43%.³⁴ Enhanced phonon scatterings at the defect locations and reduced phonon mean free paths were found to lead to a one to two orders of magnitude reduction of κ in the presence of such defects. Structural defects such as pentagonal, heptagonal, and octagonal carbon atom rings preserve the threefold coordination of carbon atoms.

The thermal conductivity of GNR depends sensitively on the arrangements (isolated, lines, or haeckelites) of these defects. In systems containing arrays of parallel line defects, the κ decrease in the direction of defect lines is similar to that caused by narrow confined spaces. By introducing linearly formed defects in graphene’s lateral directions, the calculated κ parallel to the defects is significantly larger than that in the perpendicular directions. The symmetry breaking caused by non hexagonal rings in the graphene will enhance the anisotropic thermal conductivity in GNR. The reduction of κ due to the presence of non hexagonal rings in the graphene lattice is comparable in magnitude to that caused by divacancies at the same concentration. The results indicate that the main reason for the reduction of the thermal conductivity in defective systems is a decrease in the phonon mean free path.³⁵ By introducing elliptical pores in GNRs, a preferred heat conduction direction can be created. The scale of anisotropy increases monotonically with the elliptical pores’ aspect ratio.³⁶

Another effective approach to manipulate the anisotropic thermal conductivity of graphene is chemical functionalization. It is well known that thermal transport in GNR is largely dominated by low frequency (0–30 THz) acoustic phonons.³⁷ Upon graphene hydrogenation, both the peak positions and intensities of the graphene phonon densities of state (PDOS) will be altered significantly. A significant PDOS shift in hydrogenated GNR (HGNR) by ∼10 THz together with a reduction of the peak intensity is calculated.³⁸ Another important issue related to the κ decrease is the significant reduction of all the elastic constant of HGNR with respect to GNR due to the sp² to sp³ transition. The C–C interatomic bonds in graphene will be converted from in-plane sp² bonding into out-of-plane sp³ bonding with surface hydrogenation. The HGNR can be viewed as a mixture of a high thermal conductivity sp² C–C structure and a low thermal conductivity sp³ C–C assembly. Therefore, the degradation in κ can be attributed to the introduction of sp³ bonds into the structure, which can be regarded as randomly distributed defects in a sp² plane matrix to scatter phonons. With an increasing coverage ratio of hydrogen atoms, the sp² planar network gradually diminishes and deteriorates. As a consequence, the phonon scattering increases and the thermal conductivity decreases. By randomly attaching hydrogen atoms on the surface of GNR, κ of both z-GNR and a-GNR reduce by 35% at a coverage ratio of 5%. When the coverage increases to 100% (graphane), κ reduction reaches its maximum at 68%. The random hydrogenation has the same effect in both chiral directions.³⁹ It is worth noting that the thermal conductivity of GNR becomes insensitive to the increase of hydrogen atoms when the coverage reaches 30% or above, which can be explained by the percolation theory. Since the site percolation threshold for a honeycomb lattice is 0.697,⁴⁰ the sp² bond network is interrupted when the hydrogen coverage reaches ∼30%. The thermal conductivity of GNR will be dominated by the sp³ C–C domains at this percentage or higher and become insensitive to further increases. When the hydrogen distribution on GNR is patterned, its effects on the κ of z-GNR and a-GNR will be diverse. If hydrogen atoms are distributed in stripe formations perpendicular to the heat flux direction, the κ decrease in both chiral directions will be the same. This is because the sp³ regions scatters the phonon transport along the thermal transport direction and dominated the thermal conductivity. While the hydrogen stripe is parallel to the heat flux direction, the κ of a-GNR reduces faster than that of z-GNR, since both the sp² and sp³ bonding contribute equally to the thermal transport but a-GNR has greater phonon boundary scatterings.³⁹ Similar conclusions can be also be achieved by oxygen functionalization in graphene structures.⁴¹

4. Decoupled phonon thermal transport

The exceptional thermal properties of graphene are partially due to its unique phonon transport mechanism in the 2-D system. There are three acoustic phonon branches, i.e. the in-plane longitudinal (LA) and transverse (TA) branches and the out-of-plane flexural (ZA) branch, which contribute to the thermal conductivity of graphene. Unlike isotropic materials such as bulk silicon and crystal argon, where different phonon branches always share the same energy level, the ZA phonon branch in graphene is decoupled with the lateral LA/TA branches during the thermal transport process. The peculiar phenomenon is caused by two major factors: first, the in-plane/out-of-plane phonon coupling rates are much lower than those between in-plane phonons; second, the ZA phonon branch has much higher thermal transport capability than the LA/TA phonons. Since phonon diffusion is the dominant mechanism of thermal transport in graphene, the challenge of accounting for these phonon features gives the starting impetus to its study.

For a long time it has been tacitly accepted that in-plane acoustic phonons are dominant in the thermal transport of graphene.^42–44 However, recent studies have proven that the fact is quite different. Saito et al.⁴⁵ calculated the ballistic thermal conductance of graphene by investigating the dispersion relation of phonons and electrons. They proved that the ballistic phonon conductance of graphene below about 20 K is mainly determined by the out-of-plane acoustic mode (ZA branch) and the in-plane acoustic modes (LA and TA branches) cannot be ignored above 20 K. After measuring the thermal transport of single layer graphene (SLG) supported by amorphous SiO₂, Seol et al.⁴⁶ performed a revised calculation of the calculated thermal conductivity for suspended graphene. They found that the ZA branch can contribute as much as 77% at 300 K and 86% at 100 K due to the high specific heat and long mean scattering time of ZA phonons. Recent studies by Zhang et al.^47–49 revealed the fact that, in a GNR system, the ZA branch has an unusually higher thermal conductivity than the LA and TA branches. Their study also showed that ZA ↔ ZA energy transfer is much faster than the ZA ↔ LA/TA phonon energy transfer. It has been proved that under the influence of a moving or static localized heat source, the flexural mode phonons dissipate heat much faster than phonons with a longitudinal mode or transverse mode, which gives rise to an energy inversion phenomenon at the system level.

In general, to define the local temperature, a mode-wide thermal equilibrium among different phonon modes should be reached. However, since there are huge differences in the thermal conductivities among different phonon modes in graphene, it is possible that there is a mode-wide energy difference during steady state heat conduction. A unique heating and cooling technique was developed by Zhang et al.⁴⁸ to investigate this particular phenomenon. In their study, the thermal energies are added/subtracted to/from the lateral and flexural phonons, respectively, at opposite directions. Under such a scenario, a heat flux carrying out-of-plane phonon energies is generated and another heat flux carrying in-plane phonon energies was found to coexist in the opposite direction. Configuration of the bi-directional heat conduction is illustrated in Fig. 3(a). An apparent thermal conductivity (κ_app) can be calculated given the overall temperature gradient and net heat flux in the system. After carefully adjusting the thermal energies added to the flexural and lateral phonons (μ), both positive and negative κ_app were calculated in the graphene system, as shown in Fig. 3(b). A close inspection of individual phonon modes proved that this new phenomenon doesn’t violate the Fourier’s law of heat conduction, and graphene is the only material known to have this special thermal property.

Fig. 3 (a) Configuration of the bi-directional heat conduction phenomenon in GNR. (b) Apparent thermal conductivity topology for the 2.0 × 50.1 nm² GNR system: the upper x-axis represents the net heat flux in the GNR system and the lower x-axis stands for the corresponding out-of-plane (FM) versus in-plane (LM/TM) energy ratios μ; negative κ_app is observed when μ is within the range 1.05–1.38 (with permission from Elsevier for use in this paper).⁴⁸

The ability of graphene to retain two heat currents is a direct proof of the decoupling nature among its in-plane and out-of-plane phonons. The calculated thermal relaxation time between ZA and TA/LA was found to be 4.7 times longer than that between TA and LA.⁴⁷ Moreover, it was found that the coupling strength is highly dependent on the system temperature. During transient thermal transport processes, this coupling lag would lead to a phonon energy inversion phenomenon in graphene between the in-plane and out-of-plane phonons. When phonon excitation is triggered in one of the lateral directions, the initial low-energy flexural phonons will quickly absorb the excitation energy from the in-plane directions and slowly release it back. An illustration of this energy inversion phenomenon in graphene is depicted in Fig. 4. The energy inversion can be observed when a static or moving localized hot spot exists in the graphene system. These discoveries provide helpful knowledge for designing graphene-based nano-devices working under non-uniform heating conditions.

Fig. 4 Illustration of the energy inversion phenomenon in graphene; after an instant thermal excitation in the lateral x direction (E_k,x), the initial “cold” out-of-plane phonons (E_k,z) quickly absorb the thermal energies from in-plane phonons and become the highest level phonon mode (energy inversion) in the following energy relaxation process (with permission from Elsevier for use in this paper).⁴⁷

5. Thermal rectification

For the last few decades, thermal rectification (TR) has attracted increased attention in fields like applied physics and thermal sciences. TR is defined as the capability of a material or a structure to transfer heat asymmetrically. Through a better understanding of TR in thermal interface materials (TIMs) like graphene, devices such as thermal diodes, thermal transistors and thermal logic circuits could be developed and utilized in micro/nanoelectronic cooling as well as thermal memory and computation.

Several models have been proposed to explain the mechanism of TR in graphene, which includes the thermal potential barrier at interfaces, the dependency of thermal conductivity on temperatures between adjacent materials, rough surface contacts, nanostructure asymmetry, and anharmonic lattices. Graphene-based thermal diodes have become an important building block of modern phononic devices, which employ logic calculations performed with phonons.⁵⁰ By properly designing the structure of a graphene system, a high TR ratio and fast switching speed can be obtained to improve significantly the performance of nanoelectronic devices.

Substrate coupling is the most commonly used perturbation for the thermal transport of graphene. Due to the significant contributions of a flexural phonon to graphene’s thermal transport, TR in graphene can be achieved by boosting/suppressing the out-of-plane phonons. The population of long-wavelength flexural phonons will significantly decrease in the presence of a substrate. A TR efficiency of 40% was calculated between suspended and inhomogeneous encased graphene structures.⁵¹ It was found that perturbations induced by substrate coupling significantly suppress the long wavelength flexural phonon mode in the encased graphene when compared with that in the suspended graphene. As a result, the in-plane phonons can transmit through the interfaces well, whereas low frequency flexural phonon modes are reflected, leading to a nontrivial interfacial thermal resistance. Similarly, Zhong et al.⁵² calculated the TR ratio in thickness-asymmetric GNRs. The suppression of the flexural phonons in thicker regions reduces the phonon transmission rate from single layer GNR, which causes direction-dependent heat transport.

For symmetric graphene systems, TR can also be induced by lattice distortions which can both alter the characteristic vibrational frequencies and the degree of anharmonicity. If a tensile stress or strain is non-uniformly applied to the graphene system, phonon scattering and participation rates would be different along the GNR structures. In three-dimensional bent or twisted GNRs, the phonon scatterings are furious in the deformed areas due to high local strains.⁴⁹ By applying tensile strains perpendicular to the heat flux direction of a C¹⁴ isotope doped GNR, a 45% TR was obtained. Moreover, the ratio of TR was found to increase monotonically with the strain’s magnitude.⁵³ A significant TR over 70% was calculated in a rectangular a-GNR by applying a transverse force asymmetrically.⁵⁴ The heat transport was found to be favorable from the less stressed region to the more stressed region, and TR is only weakly dependent on edge defects and randomly distributed vacancies. While on the other hand, if the vacancy defect is distributed with patterns in GNR, the TR rate can be altered depending on the density and formation of the vacancies. By introducing a 1% concentration of single-vacancy or substitutional silicon defects, and a moderate partition of the pristine and defected regions, a TR ratio of 80% was found in a 14 nm long 4 nm wide GNR at a temperature of 200 K.⁵⁵ It was found that under various sizes of triangular vacancy defect and various temperatures, the heat flux runs preferentially from the vertex to the bottom of the triangular vacancy. A TR of 16% was reported in the triangular-vacancy-defected a-GNRs.⁵⁶ Other defect types such as grain boundaries, boron/nitrogen doping and hydrogen-terminating can also induce strong TR effects in GNR depending on the distribution patterns and MD temperatures.^57–59

Another popular approach to create TR in graphene is asymmetric atomic configurations. The structural variation can cause different inelastic phonon scatterings at the boundaries, which can change the phonon power spectrum formations at different locations. The localized phonon populations can increase the three-phonon scatterings in the boundary regions and affect the thermal conductivity of graphene by changing the phonon participation rate during the thermal transport process.⁶⁰ The spatial distribution of a specific range (Λ) of normal modes can be expressed as

where N is the total number of atoms with atom number i, and α is the direction of interest. ε_iα,λ is the eigenvector component of the λth phonon mode. A larger value of ϕ_iα,Λ indicates stronger localization of modes Λ on the ith atom. The participation ratio (p_λ) characterizes phonons on a per atom basis and is defined for a specific mode λ as

The participation ratio characterizes each mode individually and serves as a useful discriminant of spatial localization. The participation ratio measures the fraction of atoms participating in a mode and hence varies between O(1) for delocalized states to O(1/N) for localized states and effectively indicates the fraction of atoms participating in a given mode. The phonon participation ratio for bulk graphene and a trapezium-shaped GNR (TGNR) is shown in Fig. 5(a). It can be observed that the p_λ of TGNR is less than that of bulk graphene, indicating stronger phonon localizations in the confined nanoribbons. The spatial distribution of the localized modes (LM) is shown in Fig. 5(b) and (c), which illustrate the edge-preference phenomenon for LM. In Fig. 5(b), phonon localization is enhanced on the wider side when heat flows in the decreasing width direction (left to right), but has a limited effect on thermal transport since the channel bottleneck is from the narrower end. On the other hand, when heat flows in the direction of increasing width (right to left), more LMs can be found on the narrower end, which acts as the bottleneck of thermal transport channel in such narrow GNRs. Therefore, delocalized modes have a narrower channel of propagation, and hence the effective κ is reduced. Based on above explanations, the TGNR selects the direction of decreasing width as the favored direction of thermal rectifiers. When the dimensions of GNR increase to bulk size, the three-phonon scattering dominates and the edge or chiral scattering becomes unimportant. Therefore, no TR is observed in bulk graphene structures. Thermal rectification in triangular asymmetric and rectangular symmetric GNRs are shown in Fig. 6.⁶¹ The TR factor can be as large as 120% at MD temperature of 180 K as shown in the inset of Fig. 6(a). A TR up to 74% was calculated in a grain boundary graphene system.⁵⁹ Thermal rectification in TGNR, U-shaped GNR (UGNR), Y-junctions GNR (YGNR) and three-terminal graphene nanojunctions (TGNJ) were investigated by several MD studies.^62–65 The calculated TR results from the match/mismatch of the PDOS between the two ends, which controls the heat current in the system. The UGNR is composed of a beam and two asymmetric arms. Due to the asymmetric boundary scattering of phonons, the heat flux runs preferentially from the wide arm to the narrow arm, which indicates a strong rectification effect. It was found that the rectification ratio is insensitive to the heat flux and beam, and the structural asymmetries can be properly designed to obtain the maximum rectification ratio. In YGNR with lengths of 16.7 nm, it was found that the heat flux runs preferentially from the branches to the stem. More interestingly, due to the presence of layer–layer interactions, a larger rectification ratio can be achieved in double-layer structures compared to single-layer graphene Y junctions, which suggests that few-layer graphene systems could be applied to develop more reliable and robust directional phononics devices. Similar conclusions can be achieved in TGNRs, which shows that heat flux runs preferentially along the direction of decreasing width. It is shown that the graded geometric asymmetry is of remarkable benefit to increase the rectification ratio in TGNR.

Fig. 5 (a) Participation ratio of bulk graphene and T-shaped GNRs; spatial distribution of localized modes when heat flows in the direction of (b) decreasing width (left to right) and (c) increasing width (right to left) (with permission from ACS Publications for use in this paper).⁶⁰

Fig. 6 Thermal conductivity of (a) triangular asymmetric and (b) rectangular symmetric GNRs; the right inset of (a) shows the thermal rectification factor η as a function of temperature (with permission from ACS Publications for use in this paper).⁶¹

The convergence of rectification ratio with the increase in length is of particular importance for practical thermal phononics applications. An extraordinarily high TR ratio of about 200% can be achieved in three-terminal graphene nanojunctions. It was shown that the heat flux runs preferentially along the direction from narrow to wide terminals. The corner form of the TGNJs plays an important role in the rectification effect. Meanwhile, the rectification efficiency increases rapidly with the width discrepancy between the left and right terminals and is strongly dependent on the formation of the nanojunctions. The reported TR ratio increases monotonically with the MD temperature. Phonon boundary scatterings become more furious at high temperatures and the high frequency optical phonons might break down into numerical acoustic phonons, which have a higher transmission probability due to the limited width of the asymmetric grain boundary. Furthermore, the increased temperature would stimulate more acoustic phonons in the graphene system and directly increase the boundary thermal conductance.

TR has also been reported in various graphene-based heterostructures. The existence of the interface TR can be explained by the asymmetry in the phonon transmission caused by the mismatch of phonon spectra across the interface. A significant TR ratio of 180% was observed in carbon-nanotube/GNR junctions due to the localized phonon modes at the boundary and at the junction region. Since rectification is a result of non-equilibrium transport, which is mainly determined by optical phonons, the calculated TR ratio increases significantly with the interfacial temperature difference. When the temperature difference is large enough, the optical phonons are excited and start contributing to the heat conduction, which then induces greater TR ratio.⁶⁶ A 44% TR is calculated at the graphene and silicene monolayer interface.⁶⁷ Since graphene’s phonon frequency range is higher than that of silicene, it can transmit higher frequency phonons. Therefore, both the low frequency and high frequency phonons originating from silicene can be accepted by graphene, whereas the high frequency phonons are blocked from graphene to silicene. In a similar study, a 23% TR was obtained at the graphene/silicene bilayer heterostructure.⁶⁸ The interface thermal conductance is larger when heat transfers from graphene to silicene than in a reverse transfer direction. Since the atom density per area of the graphene is 2.25 times larger than that of the silicene. During the temperature increase, the larger atom density per area of the graphene also brings about more phonon excitation, contributing to the interface thermal transport.

The temperature dependency of heat capacity and atom density mismatch gives rise to the TR in the heterostructure. A high rectification ratio of 140% at room temperature was achieved in GNR–carbon nanotube (CNT) junction structures.⁶⁶ The localized phonon modes at the boundary of GNR are responsible for the thermal conductivity reduction from CNT to GNR. In vertically aligned GNR–silicon interfaces, it was found that the TR factor increases with raised heat flux values, indicating the importance of temperature bias across interfaces in controlling thermal rectification.⁶⁹ Rajabpour et al.⁷⁰ investigated the thermal transport across graphene–graphane interfaces using NEMD method. The hydrogen atoms on graphane were found to behave as vibrating centers for the carbon atoms, which affects the phonon power spectrum of the graphane system and induces TR in the hybrid structure. Similarly, asymmetric chemical functionalization and isotopic doping in graphene can also induce thermal rectification at various temperatures due to the mismatch of the phonon spectra.^53,71

6. Tunable interfacial thermal conductance

When the size of a device is reduced to a length scale on the order of the energy carrier’s mean free path, device-level thermal transport is no longer determined by the thermal properties of the materials comprising the devices, but rather by the energy transport across the interfaces between adjacent materials. Thermal contact resistance is a measurement of an interface’s resistances to thermal flow and it is the most common quantity used to characterize interfacial thermal transport. Understanding the thermal resistance between two materials is of great significance to study their thermal properties. This is crucial for micro/nanoscale systems where interface phonon behaviors could directly affect the properties relative to bulk materials.

Recently there has been a great interest in the study of the interfacial thermal conductance across graphene based nanocomposites.^68,72–76 A summary of interfacial thermal resistances (R) between graphene and various substrates calculated by MD simulation are listed in Table 1. The limited internal phonon coupling and transfer within graphene in the out-of-plane direction significantly affects graphene–substrate interfacial phonon coupling and scattering, leading to unique interfacial thermal transport phenomena. In many electronic applications, the surface of a supporting material is often carved with grooves to achieve maximum thermal performance and special electrical functions. For pliable 2D materials like graphene, such structures can be easily bent to fit the surface formations of the substrate. Under such scenarios, it is generally considered that the interfacial thermal resistance can increase due to the loose contacts at the interface. However, by performing classical MD simulations, it was proved that by carefully designing the roughness patterns on the substrate surface, the thermal contact resistance between graphene and silicon can decrease by more than 10%.⁷⁶ The observed results were soundly explained with the interatomic bond strength tuning by the surface roughness. For small nanogroove depths, the C–Si bond length in the suspended graphene region is large, and graphene experiences a strong pulling-down force (attractive) from the Si substrate. On the other hand, the graphene in the supported region is in direct contact with Si and experiences an extremely strong repulsive force to balance the pulling-down force in the suspended region. The repulsive force on graphene in the supported region can reach a level of 228 MPa. This significantly increases the local energy coupling and offsets the energy coupling reduction in the suspended graphene region. In a similar MD study by Hong et al.,⁷⁷ a surface nanoengineering design was reported to reduce the R between the supported graphene and copper substrate by 17%. It was proved that interfacial thermal resistance is dependent on the widths (d) and height (δ) of the nano-pillars on the substrate surface as well as the roughness formations. Zhang et al.⁴⁹ discovered furious phonon scatterings at the bent regions in a right-angle graphene structure by performing NEMD simulations. Detailed radial distribution function (RDF) analyses were conducted to explain the results. The compressive strain in the deformed areas enhances the three-phonon scattering and indirectly facilitates the in-plane and out-of-plane phonon couplings. The bending angle (θ) plays a vital role in determining the thermal resistance, i.e., R increases with the decrease in θ. Therefore, it can be speculated that the interfacial thermal resistance between graphene and its substrate is dependent on the shapes of the surface nanobumps due to different in-plane phonon transport and surface contacts. Cylindrically shaped roughness patterns were proved to give better out-of-plane thermal dissipations than the rectangular shaped nanobumps, as shown in Fig. 7(a)–(c).

Table 1 Interfacial thermal resistance R between graphene and substrate materials calculated by molecular dynamics (MD) simulation

Team	Materials	Temperature/K	R × 10^−8/K m^2 W^−1
Zhang et al.^76	Graphene/Si	300	3.1–4.9
Yu et al.^86	Graphene/Si	300	7.5
Shen et al.^87	Graphene/Si	300	∼4.3–5
Vallabhaneni et al.^69	Graphene/Si	300	∼0.068–0.085
Li et al.^88	Graphene/SiC	300	0.091–5.6
Wang et al.^89	Graphene/SiC	300	∼0.01–0.07
Yue et al.^90	Graphene/SiC	300	0.07–0.085
Wang et al.^91	Graphene/SiC	300–600	∼1.03–1.75
Wang et al.^92	Graphene/SiC	100–700	∼0.1–3.5
Hong et al.^77	Graphene/copper	300	2.1–3.6
Wei et al.^93	Graphene/copper/nickel	230–430	∼0.2–0.83
Luo et al.^94	Graphene/polymer	298	∼1.7
Liu et al.^68	Graphene/silicene	200–700	∼4.2–16.8
Hu et al.^95	Graphene/phenolic-resin	300	∼11.2
Park et al.^96	Graphene/CNT	10–600	∼0.018–0.045
Wei et al.^97	Graphene/graphene	400–1200	∼0.021–0.025

---

Fig. 7 Atomic configurations of (a) cylindrical nanobumps and (b) rectangular nanobumps for d = 2 nm, δ = 0.83 nm heterostructures at a steady state; (c) dependence of R on nanogroove depth for cylindrical and rectangular shaped nanobump systems (with permission from Royal Society of Chemistry for use in this paper).⁷⁷

Functionalization of graphene through chemical methods has attracted great interest and has been explored recently due to the possibility of nanostructuring graphene into complex patterns.^78–80 The thermal properties of hydrogenated graphene (H-GNR) and oxidized graphene (O-GNR) were investigated by several studies to seek better thermal management and thermoelectric applications.^81–83 Hopkins et al.⁸⁴ investigated the role of interfacial bonding at metal/graphene interfaces by introducing chemical absorbates on graphene surfaces in order to increase the density of covalent bonds bridging the metal and the SLG. Interfacial thermal conductance is found to increase with oxygen functionalization due to the strengthened chemical bonding. Interfacial thermal conductance (G) across stacked GNR and silicene can be tuned by adjusting the hydrogen coverage (f) on graphene. When hydrogen atoms are randomly distributed on both sides of the GNR, G has a non-monotonic trend with f. With a hydrogen coverage of 50%, thermal conductance between GNR and silicene can be increased by 50 MW m⁻² K⁻¹.⁶⁸ The physical mechanism can be explained by phonon power spectrum analyses on the hybrid structure. Enhanced phonon couplings between C–Si atoms and improved intra-mode phonon couplings within GNR are the main reasons for the thermal conductance increase. Compared with pristine graphene, it has been proved that the interfacial thermal resistance between H-GNR and hexagonal boron nitride (h-BN) can be reduced by 76.3%, as shown in Fig. 8.⁸⁵ Since the in-plane and out-of-plane phonons in graphene are highly decoupled, the thermal energies are restrained in the lateral direction and can only be released slowly to the flexural phonons. By hydrogenating graphene in the heterostructure, it has been observed that the couplings between in-plane and out-of-plane phonons have been greatly enhanced. As the hydrogenation ratio increases, the overlap area between in-plane and out-of-plane phonons in H-GNR greatly increases. The improvement of in-plane and out-of-plane phonon couplings in graphene indirectly facilitates the thermal transport across the interface and reduces the thermal contact resistance. In the meantime, the overlap area between H-GNR and h-BN’s out-of-plane direction phonons also increases with higher hydrogenation ratio, which directly contributes to the thermal conductance across the interface.

Fig. 8 Effects of hydrogen functionalization on interfacial thermal transport between graphene and h-BN; the calculated R result for fully hydrogenated GNR is decreased by 76.3% compared with pristine GNR structures (with permission from AIP Publishing LLC for use in this paper).⁸⁵

7. Conclusions and future work

In this review, the novel thermal properties of graphene studied by classical MD simulations, including anisotropic thermal conductivity, decoupled phonon thermal transport, thermal rectification and tunable interfacial thermal resistance, are summarized and discussed. The high thermal transport capability of the flexural phonon mode and the decoupled nature among in-plane and out-of-plane phonons plays a vital role in these peculiar properties. Flexural acoustical phonons are dominant for graphene’s thermal conductivity at a broad range of system temperatures and dimensions. The phonon couplings between lateral and flexural phonons are highly temperature dependent, which causes unbalanced thermal transport rates within localized hot spots and untreated areas. With the existence of static or moving hot regions, the flexural phonons will harvest the majority of thermal energies and affect the interfacial thermal transport in nano-composites. Factors like system asymmetry, edge effects and temperature dependent thermal conductivity are investigated to explain thermal rectification. A better understanding of thermal rectification is crucial for the design and development of next generation thermal logic circuits and high performance micro/nano thermal radiators. The 2D and pliable nature of graphene makes it possible to remain in close contact with its substrate at rough interfaces. Thermal contact resistance in graphene-based heterostructures can be greatly reduced by patterned surface nano-engineering, which can be achieved by laser processing of the substrate materials in practical applications. Chemical functionalization of graphene is another effective approach to manipulate its interfacial thermal properties. The enhanced intra-mode phonon coupling facilitates the cross-plane thermal transport indirectly and reduces the thermal contact resistance.

Despite this summary of the large volume of studies using MD, many problems of the thermal properties of graphene remain unsolved. In graphene-based nano-devices, thermal energies, as in most scenarios, are non-uniformly distributed across the surface. The effects of localized hot sports on interfacial thermal transport are unclear. Also, the thermal transport across graphene and metal materials has not been well investigated. Due to the computational restrictions of the ab initio method, classical MD simulation doesn’t consider quantum effects and electron thermal transport. Some theoretical approaches, like the two-temperature model and stochastic force are proposed to solve this problem; however, the non-trivial parameter settings and complicated pre/post-processing make wide implementation difficult. To further improve the thermal performance in the cross-plane directions of graphene heterostructures, the design of the optimum substrate surface needs to be explored. There are countless pattern combinations and enormous effects could be obtained. The combined effects of chemical functionalization and rough surfaces on interfacial thermal transport are of great interest and its potential to significantly reduce the thermal contact resistance is promising.

References

1. C. Lee, X. D. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385–388.
2. A. Ballestar, P. Esquinazi, J. Barzola-Quiquia, S. Dusari, F. Bern, R. R. da Silva and Y. Kopelevich, Carbon, 2014, 72, 312–320.
3. Y. P. Bliokh, V. Freilikher and F. Nori, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 245134.
4. B. van Duppen and F. M. Peeters, Europhys. Lett., 2013, 102, 27001.
5. V. V. Ilyasov, B. C. Meshi, V. C. Nguyen, I. V. Ershov and D. C. Nguyen, Solid State Commun., 2014, 199, 1–10.
6. S. Subrina, J. Nanoelectron. Optoelectron., 2013, 8, 317–336.
7. V. Dhand, K. Y. Rhee, H. J. Kim and D. H. Jung, J. Nanomater., 2013, 2013, 158.
8. K. M. F. Shahil and A. A. Balandin, Solid State Commun., 2012, 152, 1331–1340.
9. A. A. Balandin, Nat. Mater., 2011, 10, 569–581.
10. A. A. Balandin, 2011 27th Annual Ieee Semiconductor Thermal Measurement and Management Symposium (Semi-Therm), 2011, p. 89.
11. F. Ma, H. B. Zheng, Y. J. Sun, D. Yang, K. W. Xu and P. K. Chu, Appl. Phys. Lett., 2012, 101, 111904.
12. F. Ma, Z. K. Guo, K. W. Xu and P. K. Chu, Solid State Commun., 2012, 152, 1089–1093.
13. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari and F. Mauri, Nano Lett., 2014, 14, 6109–6114.
14. M. Trushin and J. Schliemann, Phys. Rev. Lett., 2007, 99, 216602.
15. H. J. Zhang, G. Lee and K. Cho, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 115460.
16. W. Li, H. Sevincli, G. Cuniberti and S. Roche, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 041410.
17. Z. G. Fthenakis and D. Tomanek, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 125418.
18. B. Mortazavi and G. Cuniberti, Nanotechnology, 2014, 25, 215704.
19. V. Varshney, A. K. Roy, G. Froudakis and B. L. Farmer, Nanoscale, 2011, 3, 3679–3684.
20. J. Tersoff, Phys. Rev. B, 1988, 37, 6991–7000.
21. D. W. Brenner, Phys. Rev. B: Condens. Matter, 1990, 42, 9458–9471.
22. J. Tersoff, Phys. Rev. Lett., 1988, 61, 2879–2882.
23. B. W. Dodson, Phys. Rev. B, 1987, 35, 2795–2798.
24. V. K. Tewary and B. Yang, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 075442.
25. L. Lindsay and D. A. Broido, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 205441.
26. P. H. M. Böttger, L. Braginsky, V. Shklover, E. Lewin, J. Patscheider, D. G. Cahill and M. Sobiech, J. Appl. Phys., 2014, 116, 013507.
27. Y. Wang, C. Gu and X. Ruan, Appl. Phys. Lett., 2015, 106, 073104.
28. Y. Peng, J.-U. Thiele, G. Ju, T. P. Nolan, Y. Ding and A. Q. Wu, Recording layer for heat assisted magnetic recording, US pat., No. 8,507,114, 2013.
29. J. W. Jiang, J. S. Wang and B. W. Li, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 205418.
30. Z. Aksamija and I. Knezevic, Appl. Phys. Lett., 2011, 98, 141919.
31. Y. Wang, B. Qiu and X. L. Ruan, Appl. Phys. Lett., 2012, 101, 013101.
32. X. Y. Liu, F. C. Wang and H. A. Wu, Appl. Phys. Lett., 2013, 103, 071904.
33. C. Wang, Y. Liu, L. Li and H. Tan, Nanoscale, 2014, 6, 5703–5707.
34. T. Y. Ng, J. J. Yeo and Z. S. Liu, Carbon, 2012, 50, 4887–4893.
35. Z. G. Fthenakis, Z. Zhu and D. Tomanek, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 125421.
36. L. Hu and D. Maroudas, J. Appl. Phys., 2014, 116, 184304.
37. B. Qiu and X. Ruan, Appl. Phys. Lett., 2012, 100, 193101.
38. G. Barbarino, C. Melis and L. Colombo, Carbon, 2014, 80, 167–173.
39. Q. X. Pei, Z. D. Sha and Y. W. Zhang, Carbon, 2011, 49, 4752–4759.
40. P. N. Suding and R. M. Ziff, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 275–283.
41. X. Shen, X. Y. Lin, J. J. Jia, Z. Y. Wang, Z. G. Li and J. K. Kim, Carbon, 2014, 80, 235–245.
42. P. G. Klemens, Int. J. Thermophys., 2001, 22, 265–275.
43. D. L. Nika, S. Ghosh, E. P. Pokatilov and A. A. Balandin, Appl. Phys. Lett., 2009, 94, 203103.
44. D. L. Nika, E. P. Pokatilov, A. S. Askerov and A. A. Balandin, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 155413.
45. K. Saito, J. Nakamura and A. Natori, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 115409.
46. J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. S. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff and L. Shi, Science, 2010, 328, 213–216.
47. J. C. Zhang, X. W. Wang and H. Q. Xie, Phys. Lett. A, 2013, 377, 721–726.
48. J. C. Zhang, X. W. Wang and H. Q. Xie, Phys. Lett. A, 2013, 377, 2970–2978.
49. J. C. Zhang and X. W. Wang, Nanoscale, 2013, 5, 734–743.
50. S. Pal and I. K. Puri, Small, 2015, 11, 2910–2917.
51. W. Xu, G. Zhang and B. W. Li, J. Appl. Phys., 2014, 116, 134303.
52. W. R. Zhong, W. H. Huang, X. R. Deng and B. Q. Ai, Appl. Phys. Lett., 2011, 99, 193104.
53. Q. X. Pei, Y. W. Zhang, Z. D. Sha and V. B. Shenoy, Appl. Phys. Lett., 2012, 100, 101901.
54. K. G. S. H. Gunawardana, K. Mullen, J. N. Hu, Y. P. Chen and X. L. Ruan, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 245417.
55. Y. Wang, S. Chen and X. Ruan, Appl. Phys. Lett., 2012, 100, 163101.
56. P. Yang, X. L. Li, Y. F. Zhao, H. Y. Yang and S. T. Wang, Phys. Lett. A, 2013, 377, 2141–2146.
57. Q. Liang and Y. Wei, Phys. B, 2014, 437, 36–40.
58. P. Yang, X. L. Li, H. Y. Yang, X. N. Wang, Y. Q. Tang and X. M. Yuan, Appl. Phys. A: Mater. Sci. Process., 2013, 112, 759–765.
59. H.-Y. Cao, H. Xiang and X.-G. Gong, Solid State Commun., 2012, 152, 1807–1810.
60. Y. Wang, A. Vallabhaneni, J. N. Hu, B. Qiu, Y. P. Chen and X. L. Ruan, Nano Lett., 2014, 14, 592–596.
61. J. Hu, X. Ruan and Y. P. Chen, Nano Lett., 2009, 9, 2730–2735.
62. J. Cheh and H. Zhao, J. Stat. Mech.: Theory Exp., 2012, 2012, P06011.
63. G. Zhang and H. S. Zhang, Nanoscale, 2011, 3, 4604–4607.
64. N. Yang, G. Zhang and B. W. Li, Appl. Phys. Lett., 2009, 95, 033107.
65. T. Ouyang, Y. P. Chen, Y. E. Xie, X. L. Wei, K. K. Yang, P. Yang and J. X. Zhong, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 245403.
66. X. X. Ni, G. Zhang and B. W. Li, J. Phys.: Condens. Matter, 2011, 23, 215301.
67. B. Liu, J. A. Baimova, C. D. Reddy, S. V. Dmitriev, W. K. Law, X. Q. Feng and K. Zhou, Carbon, 2014, 79, 236–244.
68. B. Liu, J. A. Baimova, C. D. Reddy, A. W. K. Law, S. V. Dmitriev, H. Wu and K. Zhou, ACS Appl. Mater. Interfaces, 2014, 6, 18180–18188.
69. A. K. Vallabhaneni, B. Qiu, J. Hu, Y. P. Chen, A. K. Roy and X. Ruan, J. Appl. Phys., 2013, 113, 064311.
70. A. Rajabpour, S. M. Vaez Allaei and F. Kowsary, Appl. Phys. Lett., 2011, 99, 051917.
71. H. Fan, L. Deng, X. Yuan, J. Guo, X. Li and P. Yang, RSC Adv., 2015, 5, 38001–38005.
72. S. W. Chang, A. K. Nair and M. J. Buehler, J. Phys.: Condens. Matter, 2012, 24, 245301.
73. E. Hogfors-Ronnholm and T. Wiklund, Microb. Pathog., 2010, 49, 369–375.
74. Y. Liu, J. S. Huang, B. Yang, B. G. Sumpter and R. Qiao, Carbon, 2014, 75, 169–177.
75. K. F. Mak, C. H. Lui and T. F. Heinz, Appl. Phys. Lett., 2010, 97, 221904.
76. J. Zhang, Y. Wang and X. Wang, Nanoscale, 2013, 5, 11598–11603.
77. Y. Hong, L. Lei, X. C. Zeng and J. Zhang, Nanoscale, 2015, 7, 6286–6294.
78. T. Kuila, S. Bose, A. K. Mishra, P. Khanra, N. H. Kim and J. H. Lee, Prog. Mater. Sci., 2012, 57, 1061–1105.
79. D. W. Boukhvalov and M. I. Katsnelson, J. Phys.: Condens. Matter, 2009, 21, 344205.
80. A. K. Singh and B. I. Yakobson, Nano Lett., 2009, 9, 1540–1543.
81. H. J. Zhang, A. F. Fonseca and K. Cho, J. Phys. Chem. C, 2014, 118, 1436–1442.
82. J. Y. Kim, J. H. Lee and J. C. Grossman, ACS Nano, 2012, 6, 9050–9057.
83. S. K. Chien, Y. T. Yang and C. K. Chen, Appl. Phys. Lett., 2011, 98, 033107.
84. P. E. Hopkins, M. Baraket, E. V. Barnat, T. E. Beechem, S. P. Kearney, J. C. Duda, J. T. Robinson and S. G. Walton, Nano Lett., 2012, 12, 590–595.
85. J. Zhang, Y. Hong and Y. Yue, J. Appl. Phys., 2015, 117, 134307.
86. X. Yu, L. Zhang, X. Song, T. Xi, Y. Zhao, J. Liu, X. Yang, M. Chen and P. Yang, Int. J. Mater. Struct. Integr., 2012, 6, 65–73.
87. M. Shen, P. K. Schelling and P. Keblinski, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 045444.
88. M. Li, J. Zhang, X. Hu and Y. Yue, Appl. Phys. A: Mater. Sci. Process., 2015, 119, 415–424.
89. Q. Y. Wang, C. L. Wang, Y. Zhang and T. S. Li, Nucl. Instrum. Methods Phys. Res., Sect. B, 2014, 328, 42–47.
90. Y. N. Yue, J. C. Zhang and X. W. Wang, Small, 2011, 7, 3324–3333.
91. H. Wang, J. Gong, Y. Pei and Z. Xu, ACS Appl. Mater. Interfaces, 2013, 5, 2599–2603.
92. Z. Wang, K. Bi, H. Guan and J. Wang, J. Mater., 2014, 2014, 479808.
93. C. Shu-Wei, K. N. Arun and J. B. Markus, J. Phys.: Condens. Matter, 2012, 24, 245301.
94. T. Luo and J. R. Lloyd, Adv. Funct. Mater., 2012, 22, 2495–2502.
95. L. Hu, T. Desai and P. Keblinski, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195423.
96. J. Park and V. Prakash, J. Appl. Phys., 2014, 116, 014303.
97. Z. Y. Wei, Z. H. Ni, K. D. Bi, M. H. Chen and Y. F. Chen, Phys. Lett. A, 2011, 375, 1195–1199.

---

Footnote

† Disclaimer: This submission was written by the author acting in their own independent capacity and not on behalf of UT-Battelle, LLC, or its affiliates or successors.

---