Thermal conductivity of graphene with defects induced by electron beam irradiation

Abstract

We investigate the thermal conductivity of suspended graphene as a function of the density of defects, N_D, introduced in a controllable way. High-quality graphene layers are synthesized using chemical vapor deposition, transferred onto a transmission electron microscopy grid, and suspended over ∼7.5 μm size square holes. Defects are induced by irradiation of graphene with the low-energy electron beam (20 keV) and quantified by the Raman D-to-G peak intensity ratio. As the defect density changes from 2.0 × 10¹⁰ cm⁻² to 1.8 × 10¹¹ cm⁻² the thermal conductivity decreases from ∼(1.8 ± 0.2) × 10³ W mK⁻¹ to ∼(4.0 ± 0.2) × 10² W mK⁻¹ near room temperature. At higher defect densities, the thermal conductivity reveals an intriguing saturation-type behavior at a relatively high value of ∼400 W mK⁻¹. The thermal conductivity dependence on the defect density is analyzed using the Boltzmann transport equation and molecular dynamics simulations. The results are important for understanding phonon – point defect scattering in two-dimensional systems and for practical applications of graphene in thermal management.

---

Introduction

Graphene¹ has exceptionally high intrinsic thermal conductivity, K.^2,3 The measurements of thermal conductivity of large suspended graphene samples using the optothermal Raman technique revealed K values exceeding those of bulk graphite, which is K = 2000 W mK⁻¹ at room temperature (RT).³ Independent measurements with the optothermal Raman technique^4,5 and the scanning thermal microscopy⁶ confirmed the excellent heat conduction properties of graphene. Theoretical considerations suggest that graphene can have higher thermal conductivity than that of the graphite basal planes despite similar phonon dispersions and crystal lattice anharmonicities. This fact is attributed to an unusually long mean free path (MFP) of the long-wavelength phonons in two-dimensional (2-D) lattices.^3,7,8 Recent calculations by different methods suggested that the graphene sample size should be in the 100 μm (ref. 9 and 10) or even 1 mm (ref. 11) range in order to fully recover the intrinsic thermal conductivity limited only by the lattice anharmonicity, i.e. without phonon scattering by defects, polycrystalline grains, and edges of the samples. The intrinsic K values obtained in these works ranged from 4000–6000 W mK⁻¹ near RT.^9–11 In other terms, the high intrinsic K of graphene can be explained by the fact that the phonon Umklapp scattering is less efficient in restoring thermal equilibrium in 2-D systems than in bulk three-dimensional (3-D) systems.^9,12,13

The thermal conductivity of graphene can be degraded by defects such as polymer residue from nanofabrication,¹⁴ edge roughness,⁸ polycrystalline grain boundaries,¹⁵ and disorder from contact with a substrate or a capping layer.^16–18 For this reason, the thermal conductivity of graphene synthesized by the chemical vapor deposition (CVD) is always lower than that of the mechanically exfoliated graphene from highly ordered pyrolytic graphite (HOPG).^2,4,19–21 A possible loss of polycrystalline grain orientation in the average quality CVD graphene can lead to additional degradation of the thermal conductivity.²² However, to date, there have been no quantitative experimental studies of the thermal conductivity dependence on the concentration of defects, N_D, in graphene. The only reported experimental study of the phonon – point-defect scattering in graphene utilized isotopically modified graphene.²³ The phonon scattering on isotope impurities is limited to the mass-difference term only. It does not include the local strain effects owing to missing atoms, bond breaking or presence of chemical impurities. It was established in ref. 23 that the dependence of the thermal conductivity on the isotope impurity (¹³C) concentration is in line with the prediction of the well-established virtual crystal model²⁴ used to calculate thermal conductivity in alloy semiconductors such as Si_xGe_1−x ²⁴ or Al_xGa_1−xAs.²⁵ This model predicts the highest K for the material with either x = 0 or (1 − x) = 0 and a fast decrease to a minimum as x deviates from 0. The situation is expected to be different in materials with defects induced by irradiation.

The knowledge of the K dependence on the concentration of defects induced by irradiation can shed light on the strength of the phonon – point defect scattering in 2-D materials. The change in the dimensionality results in different dependencies of the scattering rates on the phonon wavelengths in the processes of phonon relaxation by defects and grain boundaries.^8,26,27 In bulk 3-D crystals, the phonon scattering rate on point defects, 1/τ_P, varies as ∼1/f⁴ (where f is the phonon frequency).²⁷ Owing to the changed phonon density of states (PDOS), the phonon scattering rate in 2-D graphene has a different frequency dependence, 1/τ_P ∼ 1/f³, which can, in principle, affect the phonon MFP and the thermal conductivity. In addition to the fundamental scientific interest, a quantitative study of the dependence of K on N_D is important for practical applications of graphene in thermal management. The graphene and few-layer graphene (FLG) heat spreaders^28–30 will likely be produced by CVD while FLG thermal fillers in thermal interface materials (TIMs)^31–33 will be synthesized via the liquid phase exfoliation (LPE) technique. Both methods typically provide graphene with a large density of defects than that exfoliated from HOPG.

Experimental details

We report the results of an investigation of the thermal conductivity of suspended CVD grown single layer graphene (SLG) as a function of the density of defects, N_D. The unusual non-monotonic dependence of K on N_D is analyzed within the Boltzmann transport equation (BTE) approach and the molecular dynamics (MD) simulations. The samples for this study were grown by CVD on copper foils³⁴ and transferred onto gold transmission electron microscopy (TEM) grids with 7.5 μm × 7.5 μm square holes. Only the holes fully covered with graphene were chosen for the study to simplify the data extraction in the optothermal Raman technique.^2,35 In this technique, a laser serves as source of heat and local temperature rise is determined from the shift in the Raman G peak position.² The technique and its modifications^4,36 have been verified with other methods^37,38 and extended to a wide range of 2-D materials.^39,40 In the present measurements, the gold TEM grid (diameter – 3.05 mm, thickness ∼ 25 μm) served as both the heat sink and the support for the suspended graphene, in a way similar to the experiments reported in ref. 2 and 4. The high thermal conductivity of gold (K = 350 W mK⁻¹) and strong attachment of graphene to gold grid ensured the accuracy of measurements. A scanning electron microscopy (SEM) image of the TEM grid is shown in Fig. 1. The holes have black color while the gold parts appear yellow. One can notice that some holes are partially covered with graphene as seen from its greenish color. To ensure the accuracy and reproducibility of the results, the present study was conducted on three squares, which were completely covered with graphene.

Fig. 1 Scanning electron microscopy image of graphene transferred on gold TEM grid showing 7.5 μm array of square holes. Some holes are fully or partially covered with the graphene flake. The grid is depicted in gold color, the holes are shown in black and the almost transparent greenish areas are suspended graphene flakes.

The optothermal Raman technique is a non-contact steady-state technique, which directly measures the thermal conductivity.^2,3 The micro-Raman spectrometer acts both as a heater and thermometer. The measurement is done in two steps: the calibration procedure and the power-dependent Raman measurement. During the calibration, the Raman spectrum of graphene sample is recorded under low-power laser excitation in a wide temperature range.² In order to do this, the sample is placed inside a cold–hot cell (Linkam 600), where the temperature is controlled externally with steps of 10 °C and accuracy of ∼0.1 °C. The samples are kept at least five minutes at each step to stabilize the temperature, and then the Raman G peak positions are recorded. The calibration Raman measurements are performed at low laser excitation power of ∼0.5 mW to avoid any local heating caused by the laser. The procedure provides the position of the G peak as a function of the sample temperature. In the second step of the optothermal measurements, the excitation laser power is intentionally increased to cause local heating in the suspended graphene. The spectral position of the Raman G peak reveals the local temperature rise in response to the laser heating with the help of the calibration curve.^2,3

Fig. 2 shows representative calibration (a) and power measurement (b) results. One can see from Fig. 2(a) that the dependence of the G peak spectral position on the sample temperature can be approximated as linear in the examined temperature interval. The extracted temperature coefficient χ_G = −0.013 cm⁻¹ °C⁻¹ is in line with previous reports for graphene.⁴¹ The G-peak shift with increasing laser power is presented in Fig. 2(b). One should note the excellent linear dependence of the G-peak shift on the laser power. The portion of light absorbed by suspended graphene, which causes the local heating, was measured directly by placing a power meter (Ophir) under the sample. To ensure accuracy, the absorbed power was measured for a graphene covered hole and on a reference empty hole. The difference in power readings corresponds to the power absorbed by graphene at a given laser wavelength λ. The measurement was repeated ten times at different laser power levels to determine the absorption coefficient of 5.68% ± 0.72% at the excitation laser wavelength of λ = 488 nm. The light absorption coefficient at λ < 500 nm, used in our experiments, is larger than the well-known long-wavelength limit, and it can increase further owing to surface contamination, defects, and bending.^3,42–44

Fig. 2 Raman spectroscopy data for extraction of thermal conductivity of suspended CVD graphene flakes. (a) Calibration dependence of the Raman G peak position as a function of temperature. The measurement was conducted before graphene exposure to the electron beam. The inset shows a representative Raman spectrum of CVD graphene. (b) Raman G peak position dependence on the power on the excitation laser. The SEM image of this sample is depicted in the inset. The results demonstrate an excellent linear trend.

The slope of the ω_G(ΔP) curve in Fig. 2(b) contains information about the value of thermal conductivity K, which can be extracted by solving the heat diffusion equation, knowing the sample geometry and temperature rise ΔT = χ_G⁻¹Δω_G (where Δω_G is the shift in the spectral position of G peak ω_G). The large sample size ensures that the phonon transport is diffusive or partially diffusive. The “grey” phonon MFP in graphene is around ∼800 nm near RT.³ The sample size of ∼7.5 μm ensures that phonons scatter several times before reaching the edges. The details of the K extraction procedure are provided in the ESI.† The thermal conductivity of suspended CVD graphene before introduction of defects was found to be ∼1800 W mK⁻¹ near RT. This value is in agreement with the previous independent reports for suspended CVD graphene.^4,5 A possible presence of few grain boundaries and defects, introduced during synthesis or transfer, reduce the thermal conductivity of CVD graphene as compared to that of graphene obtained by mechanical exfoliation from HOPG.^2,3,20,21

The additional defects in the suspended graphene were introduced in a controllable way using low-energy electron beam irradiation.^45,46 The samples were exposed to 20 keV electron beams (SEM XL-30) with the beam current varying from ∼3 nA to 10 nA. The irradiated area was kept constant at 6.6 × 10⁷ nm² during the whole process. The irradiation dose was controlled by changing the beam current and irradiation time. The beam current was measured before each irradiation step using a Faraday cup. The details of the irradiation procedures are provided in the Methods section. The Raman spectra of the suspended graphene samples were recorded after each irradiation step. Fig. 3 shows the evolution of the D and D′ peaks in the Raman spectrum of single layer graphene after each irradiation step. One can see from Fig. 3 that the D to G peak intensity ratio, I_D/I_G, increases from 0.13 for as-grown CVD graphene all the way to 1.00, after four steps of the electron beam irradiation. The presence of the D peak in the spectrum before irradiation indicates a background defect concentration characteristic for CVD graphene and explains K values somewhat below the bulk graphite limit.^3,47 The evolution of the Raman spectrum under irradiation was used for quantifying the density of defects, N_D, following a conventional formula:^48,49

Fig. 3 Evolution of Raman spectrum under electron beam irradiation. As the sample is exposed to the electron beam, the Raman D peak intensity increases resulting in a D-to-G peak intensity ratio change from ∼0.13 to ∼1.00. The Raman G peak shifts to higher frequencies and the D′ peak appears at ∼1620 cm⁻¹. The Raman D to G peak intensity ratio is used to quantify the amount of induced defects.

It is known that eqn (1) is valid for a relatively low defect-density regime. This criterion was met in the reported experiments. The defect density increases linearly with the Raman D to G peak intensity ratio. To show the correlation between the density of defects and the electron beam irradiation dose, we have plotted the Raman D to G peak intensity ratio, I_D/I_G, as a function of the total irradiation dose (see Fig. 4). The linear dependence is clearly seen as expected for the low defect-density regime.⁴⁵ The optothermal Raman measurements were performed after each irradiation step. The temperature coefficient of the Raman G-peak, χ_G, was not significantly affected by the defect density. In measuring the ω_G(ΔP) dependence, we had to keep the power level small enough in order to avoid local healing of defects via heating.

Fig. 4 Correlation of the Raman D-to-G peak intensity ratio with the electron beam irradiation dose. The low energy 20 keV electron beam was used to irradiate graphene. The beam current varied from ∼3 to ∼9 nA. The Raman D-to-G peak intensity ratio depends linearly on the irradiation dose.

In Fig. 5 we present the extracted thermal conductivity, K, as a function of the defect density, N_D, by squares, circles and triangles corresponding to three suspended flakes of graphene. The details of the thermal data extraction have been reported by some of us elsewhere⁷ and are briefly summarized in the ESI.† For the small defect densities, N_D < 1.2 × 10¹¹ cm⁻², the thermal conductivity decreased with increasing N_D. It can be approximated with the linear dependence K = 1990 − 116 × N_D [W mK⁻¹]. In the N_D = 0 limit, the thermal conductivity K = 1990 W mK⁻¹ was still smaller than that of the ideal basal plane of HOPG due to the background defects and possible grain boundaries present in CVD graphene before irradiation. The presence of defects before irradiation was evidenced from D peak in the Raman spectrum. At the defect density of N_D ∼ 1.5 × 10¹¹ cm⁻², one can see an intriguing change in the K(N_D) slope. It can be interpreted as a strong reduction in the rate of the decrease of K with increasing concentration N_D or the on-set of saturation. The thermal conductivity in this region is still rather high K ∼ 400 W mK⁻¹. This is clearly above the amorphous carbon limit.³

Fig. 5 Dependence of the thermal conductivity on the density of defects. The experimental data are shown by squares, circles and triangles. The solid curves are calculated using the BTE with different values of the specularity parameter p. Note that the interplay of three phonon relaxation mechanisms – Umklapp, point-defect, and rough edge scattering – gives a thermal conductivity dependence on the defect density close to the experimentally observed trend.

Discussion

For theoretical interpretation of the measured behavior of the thermal conductivity we employed both a BTE analysis and MD simulations. The details of our BTE approach are provided in the Methods section. For analysis of our experimental data, we take into account three main mechanisms of phonon scattering: phonon–phonon Umklapp (U) scattering, phonon – rough edge scattering (also referred to as boundary (B) scattering), and phonon – point-defect (PD) scattering. Within the relaxation time approximation (RTA), the total relaxation rate is given as:

1/τ_tot(s,q) = 1/τ_U(s,q) + 1/τ_B(s,q) + 1/τ_PD(s,q), (2)

where the index s = LA, TA, or ZA enumerates longitudinal acoustic (LA), transverse acoustic (TA), and out-of-plane acoustic (ZA) phonon polarization branches, and q is the phonon wave number. The dependence of the thermal conductivity on the defect density, calculated from BTE within the RTA for different values of the specularity parameter p, is presented in Fig. 5 by solid curves. The specularity parameter depends on the roughness of the edge and defines the probability of specular scattering of the phonons. For p = 1, the scattering of phonons is purely specular, which means that the edge scattering does not introduce extra thermal resistance. For p = 0, the scattering is fully diffuse, which corresponds to the strongest thermal resistance from the graphene edges.^8,9 The experimentally observed trend in thermal conductivity can be recovered with the reasonable values of the specularity parameter changing from p = 0.5 to p = 0.9.⁸

The strength of the phonon scattering on defects is determined by the mass-difference parameter ζ = (ΔM/M)², where M is the mass of carbon atom and ΔM = M − M_D is the difference between masses of a carbon atom and a defect. The value of ζ strongly depends on the nature of defects. In our BTE analysis, we used ζ as a fitting parameter to the experimental data. Within our model assumptions, the agreement with the experimental results is reached for ζ = 590. The perturbation theory calculations⁵⁰ for pure vacancy defects in graphite estimate the value of the parameter to be ζ ∼ 9. This is substantially smaller than our fitting to the experimental data. The latter is explained by the fact that our model assumes only one type of phonon – defect scattering: mass-difference scattering on single vacancies. In reality, our samples contain a variety of defects, including those that were present before irradiation and those induced by irradiation, which are different from simple vacancies. Thus, large ζ imitates the effect of phonon scattering on all other types of defects. The expected defect clustering will also result in higher ζ than that calculated from the perturbation theory under point-defect assumption. The important conclusion from the BTE modeling is that the observed weakening K(N_D) dependence can be reproduced via interplay of the three main phonon scattering mechanisms – Umklapp scattering due to lattice anharmonicity, mass-difference scattering, and rough edge scattering.

Let us now consider a possible nature of defects in our samples and their effect on the thermal conductivity as revealed from MD simulations. The details of our MD computational procedures are given in the Methods section. The electron energies of 20 keV used in the electron beam irradiation process are less than the knockout threshold energy of 80 keV.^45,51–53 Such irradiation is only sufficient to overcome the energy barrier required for breaking of the carbon–carbon bond and initiating reaction with any residual impurities such as H₂O and O₂ on the surface of graphene. This reaction results in functionalization of graphene with –OH and –C=O groups. Prior studies have shown that the –C=O configuration is energetically more favorable than –OH, and the transition of –OH and other functional groups into the energetically stable –C=O configuration can occur especially when they are annealed.⁵⁴ The energy barrier for the diffusion of –OH and epoxy groups is around 0.5–0.7 eV,⁵⁵ which corresponds to a diffusion rate ∼ 10² s⁻¹ as calculated from transition-state theory, assuming a typical phonon frequency range in graphene. For this reason, the functional groups can be mobile at the temperature of the thermal experiments (∼350 K). Upon continuous electron beam irradiation, two epoxy or hydroxyl group can come together and release an O₂ molecule.⁵⁵ When the coverage of functional groups is high, detectable amounts of CO/CO₂ can be released creating vacancies in the graphene lattice.⁵⁶ The presence of –OH and –C=O functional groups can be the reason for stronger phonon – defect scattering than that predicted by BTE models with vacancies only (and the resulting large ζ required for fitting to the experimental data). Our MD simulations show that a combination of single and double vacancy defects can also account for the experimentally observed thermal conductivity dependence on the defect concentration. The absolute value at the zero-defect limit is lower than the experimental due to the domain-size limitation in the simulation.

As one can see from Fig. 6, the thermal conductivity decreases drastically for N_D increasing from 2 × 10¹⁰ cm⁻² to 10 × 10¹⁰ cm⁻² and subsequently reaches a near-constant value at the higher concentrations of defects. This value is substantially above the amorphous carbon limit – in line with the experiment. According to this model scenario, upon irradiation, –C=O and other functionalized defects are formed that strongly reduce the thermal conductivity. Continuous irradiation results in the creation of single and double vacancies. The increase in their concentration does not lead to pronounced K reduction, which approaches an approximately constant value for the N_D range that was investigated. It can be explained in the following way. As more defects are introduced in graphene through irradiation the additional defect sites serve as scattering centers for phonons with wavelengths shorter than the distance between two vacancies. The delocalized long-wavelength phonons, that carry a significant fraction of heat, are less affected by extra defects that are closely spaced compared to those introduced at the previous irradiation step. At some irradiation dose, the increase in the phonon scattering rate of the delocalized modes due to extra defects is substantially smaller than that of the short-ranged localized modes. Hence, after a certain critical N_D the thermal conductivity effectively saturates. The weakening of the K(N_D) dependence observed experimentally and revealed in the present MD simulation is in line with reported computational results performed for graphene and graphene ribbons under various assumptions about the nature of defects.^57–60

Fig. 6 Molecular dynamics simulation results for thermal conductivity of graphene with single and double vacancy defects. The simulated defect structures are depicted in the inset. The results show that the contributions of single and double vacancies are similar in reducing the thermal conductivity of graphene. The results are in line with the experimental trend. The absolute value at the zero-defect limit is lower than the experimental due to the domain-size limitation in the simulation.

We further analyzed experimental Raman data to confirm the presence of vacancies in the irradiated graphene following the methodology developed in ref. 61 In this approach, the type of defects is determined from the ratio of intensities of D and D′ peaks, I(D)/I(D′). It has been shown that I(D)/I(D′) in graphene attains its maximum (≃13) for the defects associated with sp³ hybridization, decreases for the vacancy-like defects (≃7), and reaches a minimum for the boundary-like defects (≃3.5).⁶¹ Following this method,⁶¹ the presence of vacancy type defects has been confirmed in our irradiated graphene sample (I(D)/I(D′) ≃7). The details of this analysis are provided in ESI.† It is known that the intensity of the D band depends not only on the concentration of defects,⁶² but also on the type of defects, and only defects that are capable of scattering electrons between the two valleys K and K′ of the Brillouin zone can contribute to the D band.^63–65 For this reason, not all types of defects in graphene can be detected by Raman spectroscopy. However, our Raman data confirm the presence of vacancies supporting the theoretical assumptions. A recent study⁶⁶ suggested that the actual thermal conductivity of graphene can be even higher than that obtained from the optothermal technique using the standard procedures. Since our study is focused on the relative change in the thermal conductivity due to defect introduction, this possibility does not substantially affect the above discussion.

Conclusions

We investigated the thermal conductivity of suspended CVD graphene as a function of the defect density. The defects were introduced by the low-energy electron beams and quantified by the Raman D-to-G peak intensity ratios. It was found that as the defect density changes from 2.0 × 10¹⁰ cm⁻² to 1.8 × 10¹¹ cm⁻² the thermal conductivity reduces from ∼(1.8 ± 0.2) × 10³ W mK⁻¹ to ∼(4.0 ± 0.2) × 10² W mK⁻¹ near RT. At higher defect density the thermal conductivity revealed an intriguing weakening of the K(N_D) dependence. This behavior was explained theoretically within the Boltzmann transport equation and molecular dynamics approaches. The obtained results contribute to understanding the acoustic phonon – point defect scattering in 2-D materials. Our data indicating rather large values of thermal conductivity for graphene with defects adds validity to the proposed practical applications of graphene in thermal management.

Methods

Graphene synthesis and transfer

The single layer graphene samples were synthesized using ambient pressure chemical vapor deposition (AP-CVD) on a Cu foil.^34,67 A polycrystalline Cu foil (99.8%, Alfa Aesar) was cleaned in acetic acid, acetone and IPA to remove any surface oxides. The cleaned Cu foil was loaded into the CVD chamber and the furnace temperature was ramped to 1030 °C while flowing Ar and H₂ and the foil was annealed for 2 h. For the growth of graphene, methane (90 ppm) along with Ar and H₂ was introduced into the chamber for 20 min. After the growth, the furnace was turned off and cooled to room temperature in Ar and H₂ atmosphere. Then, the SLG grains were transferred on to a gold TEM grid using a direct transfer method to avoid any contamination from the polymer support layer. A TEM grid (G2000, 7.5 μm square holes, TedPella) was placed directly on the Cu foil–graphene stack along with a drop of isopropyl alcohol (IPA). Upon heating, as IPA evaporates, the surface tension draws graphene and the metallic grid together into intimate contact. The Cu foil was then etched in ferric chloride, washed in DI water, and the resulting graphene on TEM grid was dried for use in subsequent Raman measurements.

Electron beam irradiation

The samples were irradiated under 20 keV electron beam using Philips XL-30 FEG field-emission system. The suspended graphene sample was exposed to continuous electron beam from electron gun with current varying from ∼3 nA to ∼10 nA controlled by the beam spot size. Before each irradiation step, the Faraday cup was used to read the beam current at the desired spot size. A constant magnification was maintained during all irradiation steps in order to keep the irradiated area constant (6.6 × 10⁷ nm²). As a result, the dose density was controlled by the irradiation time. The irradiation process was done inside a vacuum chamber with the pressure below 10⁻⁴ Torr.

Boltzmann transport equation approach

In order to analyze the experimental data we used the BTE approach with the relaxation time approximation. In the framework of this BTE–RTA approach the thermal conductivity can be written as:^8,68where h = 0.335 nm is the graphene layer thickness, and the summation is performed over all acoustic phonon branches s = LA, TA or ZA, ω_s is the phonon frequency of the s-th phonon branch, q is the phonon wave number, τ_tot(s,q) is the total phonon relaxation time, T is the absolute temperature, ħ and k_B are Plank's and Boltzmann's constant, respectively. The scattering rates for the three main phonon relaxation processes, phonon–phonon Umklapp (U) scattering, phonon – rough-edge scattering (B) and phonon – point-defect (PD) scattering, are given by:

Here ν_s = dω_s/dq is the phonon group velocity, p is the specularity parameter introduced above, S is the surface per atom, ω_s,max is the maximum cut-off frequency for a given branch, γ_s is an average Gruneisen parameter of the branch s, M is the mass of an unit cell, Γ = ζ(N_D/N_G) is the measure of the strength of the point defect scattering and N_G = 3.8 × 10¹⁵ cm⁻² is the concentration of carbon atoms.

Molecular dynamics

Simulations are performed on a pristine graphene sheet of size 319.5 nm × 54.1 nm containing 660 000 carbon (C) atoms. Defects (single and double vacancies) are introduced in the structure by randomly selecting and removing carbon atoms. The C–C interactions are described using the optimized Tersoff potential for thermal transport in graphene.⁶⁹ Periodic boundary conditions are employed in all directions. The simulations are carried out with the LAMMPS package.⁷⁰ The graphene structure is energy minimized and subsequently simulated under the isothermal–isobaric (NPT) ensemble using the Nose–Hoover thermostat at 300 K and barostat at 0 MPa for 4 ns, followed by equilibration in canonical (NVT) ensemble for 4 ns using the Nose–Hoover thermostat at 300 K. The coupling time for thermostats are 0.1 ps and that for barostat is 1 ps. The thermodynamic constraints are removed and the structure is simulated under the microcanonical (NVE) ensemble for 3 ns to ensure equilibration. Subsequently, thermal conductivity is computed using the reverse non-equilibrium MD technique.^71,72 The time step of integration used in all the simulations is 1 fs.

Author contributions

A. A. B. initiated and coordinated the project, led experimental data analysis and wrote the manuscript; H. M. performed material characterization and Raman optothermal measurements; P. R. prepared suspended CVD graphene samples; A. M. supervised material synthesis and contributed to data analysis; D. L. N performed BTE simulation of heat conduction and contributed to data analysis; S. S. and G. B. performed the atomistic simulations and the corresponding thermal conductivity analysis; R. L. contributed to the computational data analysis. All authors contributed to manuscript preparation.

Competing interest

The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to (A. A. B.): balandin@ece.ucr.edu

Acknowledgements

The work at UC Riverside was supported, in part, by the National Science Foundation (NSF) grant CMMI 1404967 on engineering defects in designer materials; NSF grant ECCS 1307671 on thermal properties of graphene, and by the Defense Advanced Research Project Agency (DARPA) and Semiconductor Research Corporation (SRC) via STARnet Center for Function Accelerated nanoMaterial Engineering (FAME). The work at Iowa State University was supported by the NSF grant CMMI-1404938 on engineering defects in designer materials. AAB acknowledges useful discussions with Professor A. V. Krasheninnikov on the nature of defects in graphene.

References

1. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
2. A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao and C. N. Lau, Nano Lett., 2008, 8, 902–907.
3. A. A. Balandin, Nat. Mater., 2011, 10, 569–581.
4. W. Cai, A. L. Moore, Y. Zhu, X. Li, S. Chen, L. Shi and R. S. Ruoff, Nano Lett., 2010, 10, 1645–1651.
5. S. Chen, A. L. Moore, W. Cai, J. W. Suk, J. An, C. Mishra, C. Amos, C. W. Magnuson, J. Kang, L. Shi and R. S. Ruoff, ACS Nano, 2010, 5(1), 321–328.
6. K. Yoon, G. Hwang, J. Chung, H. G. Kim, O. Kwon, K. D. Kihm and J. S. Lee, Carbon, 2014, 6, 77–83.
7. S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau and A. A. Balandin, Nat. Mater., 2010, 9, 555–558.
8. D. L. Nika, E. P. Pokatilov, A. S. Askerov and A. A. Balandin, Phys. Rev. B: Condens. Matter, 2009, 79(15), 155413.
9. D. L. Nika, A. S. Askerov and A. A. Balandin, Nano Lett., 2012, 12, 3238–3244.
10. S. Mei, L. N. Maurer, Z. Aksamija and I. Knezevic, J. Appl. Phys., 2014, 116, 164307.
11. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari and F. Mauri, Nano Lett., 2014, 14, 6109–6114.
12. S. Lepri, R. Livi and A. Politi, Phys. Rep., 2003, 377, 1–80.
13. D. L. Nika and A. A. Balandin, J. Phys.: Condens. Matter, 2012, 24, 233203.
14. M. T. Pettes, I. Jo, Z. Yao and L. Shi, Nano Lett., 2011, 11, 1195–1200.
15. A. Y. Serov, Z. Y. Ong and E. Pop, Appl. Phys. Lett., 2013, 102, 033104.
16. K. Saito and A. Dhar, Phys. Rev. Lett., 2010, 104(4), 040601.
17. J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido and N. Mingo, Science, 2010, 328, 213–216.
18. Z. Y. Ong and E. Pop, Phys. Rev., 2011, 84(7), 075471.
19. L. A. Jauregui, Y. Yue, A. N. Sidorov, J. Hu, Q. Yu, G. Lopez, R. Jalilian, D. K. Benjamin, D. A. Delk and W. Wu, ECS Trans., 2010, 28, 73–83.
20. S. Ghosh, I. Calizo, D. Teweldebrhan, E. P. Pokatilov, D. L. Nika, A. A. Balandin, W. Bao, F. Miao and C. N. Lau, Appl. Phys. Lett., 2008, 92, 151911.
21. H. Li, H. Ying, X. Chen, D. L. Nika, A. I. Cocemasov, W. Cai, A. A. Balandin and S. Chen, Nanoscale, 2014, 6, 13402.
22. Z. Aksamija and I. Knezevic, Phys. Rev. B: Condens. Matter, 2014, 90, 035419.
23. S. Chen, Q. Wu, C. Mishra, J. Kang, H. Zhang, K. Cho, W. Cai, A. A. Balandin and R. S. Ruoff, Nat. Mater., 2012, 11(3), 203–207.
24. B. Abeles, Phys. Rev., 1963, 131(5), 1906.
25. W. Liu and A. A. Balandin, J. Appl. Phys., 2005, 97(7), 73710–73710.
26. P. G. J. Klemens, Wide Bandgap Mater., 2000, 7, 332.
27. P. G. Klemens and D. F. Pedraza, Carbon, 1994, 32, 735.
28. Z. Yan, G. Liu, J. M. Khan and A. A. Balandin, Nat. Commun., 2012, 3, 827.
29. S. Subrina, D. Kotchetkov and A. A. Balandin, IEEE Electron Device Lett., 2009, 30, 1281.
30. Z. Gao, Y. Zhang, Y. Fu, M. M. F. Yuen and J. Liu, Carbon, 2013, 61, 342–348.
31. K. M. F. Shahil and A. A. Balandin, Nano Lett., 2012, 12, 861–867.
32. V. Goyal and A. A. Balandin, Appl. Phys. Lett., 2012, 100, 073113.
33. K. M. F. Shahil and A. A. Balandin, Solid State Commun., 2012, 152, 1331–1340.
34. Q. Yu, L. A. Jauregui, W. Wu, R. Colby, J. Tian, Z. Su, H. Cao, Z. Liu, D. Pandey, D. Wei and T. F. Chung, Nat. Mater., 2011, 10(6), 443–449.
35. S. Ghosh, D. L. Nika, E. P. Pokatilov and A. A. Balandin, New J. Phys., 2009, 11, 095012.
36. H. Malekpour, K. H. Chang, J. C. Chen, C. Y. Lu, D. L. Nika, K. S. Novoselov and A. A. Balandin, Nano Lett., 2014, 14(9), 5155–5161.
37. V. E. Dorgan, A. Behnam, H. J. Conley, K. I. Bolotin and E. Pop, Nano Lett., 2013, 13, 4581–4586.
38. R. Murali, Y. Yang, K. Brenner, T. Beck and J. D. Meindl, Appl. Phys. Lett., 2009, 94, 243114.
39. R. Yan, J. R. Simpson, S. Bertolazzi, J. Brivio, M. Watson, X. Wu, A. Kis, T. Luo, A. R. H. Walker and H. G. Xing, ACS Nano, 2014, 8(1), 986–993.
40. Z. Yan, C. Jiang, T. R. Pope, C. F. Tsang, J. L. Stickney, P. Goli, J. Renteria, T. T. Salguero and A. A. Balandin, J. Appl. Phys., 2013, 114(20), 204301.
41. I. Calizo, A. A. Balandin, W. Bao, F. Miao and C. N. Lau, Nano Lett., 2007, 7(9), 2645–2649.
42. K. F. Mak, M. Y. Sfeir, Y. Wu, C. H. Lui, J. A. Misewich and T. F. Heinz, Phys. Rev. Lett., 2008, 101(19), 196405.
43. I. Santoso, R. S. Singh, P. K. Gogoi, T. C. Asmara, D. Wei, W. Chen, A. T. Wee, V. M. Pereira and A. Rusydi, Phys. Rev. B: Condens. Matter, 2014, 89(7), 075134.
44. G. X. Ni, H. Z. Yang, W. Ji, S. J. Baeck, C. T. Toh, J. H. Ahn, V. M. Pereira and B. Özyilmaz, Adv. Mater., 2014, 26(7), 1081–1086.
45. D. Teweldebrhan and A. A. Balandin, Appl. Phys. Lett., 2009, 94(1), 013101.
46. G. Liu, D. Teweldebrhan and A. A. Balandin, IEEE Trans. Nanotechnol., 2011, 10(4), 865–870.
47. C. Y. Ho, R. W. Powell and P. E. Liley, J. Phys. Chem. Ref. Data, 1972, 1(2), 279–421.
48. L. G. Cançado, A. Jorio, E. M. Ferreira, F. Stavale, C. A. Achete, R. B. Capaz, M. V. O. Moutinho, A. Lombardo, T. S. Kulmala and A. C. Ferrari, Nano Lett., 2011, 11(8), 3190–3196.
49. M. M. Lucchese, F. Stavale, E. M. Ferreira, C. Vilani, M. V. O. Moutinho, R. B. Capaz, C. A. Achete and A. Jorio, Carbon, 2010, 48(5), 1592–1597.
50. C. A. Ratsifaritana and P. G. Klemens, Int. J. Thermophys., 1987, 8(6), 737–750.
51. J. Kotakoski, A. V. Krasheninnikov, U. Kaiser and J. C. Meyer, Phys. Rev. Lett., 2011, 106(10), 105505.
52. F. Banhart, J. Kotakoski and A. V. Krasheninnikov, ACS Nano, 2010, 5(1), 26–41.
53. T. Feng, X. Ruan, Z. Ye and B. Cao, Phys. Rev. B: Condens. Matter, 2015, 91(22), 224301.
54. A. Bagri, R. Grantab, N. V. Medhekar and V. B. Shenoy, J. Phys. Chem. C, 2010, 114(28), 12053–12061.
55. S. Zhou and A. Bongiorno, Sci. Rep., 2013, 3, 2484.
56. R. Larciprete, S. Fabris, T. Sun, P. Lacovig, A. Baraldi and S. Lizzit, J. Am. Chem. Soc., 2011, 133(43), 17315–17321.
57. H. Zhang, G. Lee and K. Cho, Phys. Rev. B: Condens. Matter, 2011, 84(11), 115460.
58. F. Hao, D. Fang and Z. Xu, Appl. Phys. Lett., 2011, 99(4), 041901.
59. B. Mortazavi and S. Ahzi, Carbon, 2013, 63, 460–470.
60. J. Haskins, A. Kınacı, C. Sevik, H. Sevinçli, G. Cuniberti and T. Çağın, ACS Nano, 2011, 5(5), 3779–3787.
61. A. Eckmann, A. Felten, A. Mishchenko, L. Britnell, R. Krupke, K. S. Novoselov and C. Casiraghi, Nano Lett., 2012, 12(8), 3925–3930.
62. P. Venezuela, M. Lazzeri and F. Mauri, Phys. Rev. B: Condens. Matter, 2011, 84(3), 035433.
63. L. G. Cancado, M. A. Pimenta, B. R. A. Neves, M. S. S. Dantas and A. Jorio, Phys. Rev. Lett., 2004, 93(24), 247401.
64. C. Casiraghi, A. Hartschuh, H. Qian, S. Piscanec, C. Georgi, A. Fasoli, K. S. Novoselov, D. M. Basko and A. C. Ferrari, Nano Lett., 2009, 9(4), 1433–1441.
65. B. Krauss, P. Nemes-Incze, V. Skakalova, L. P. Biro, K. V. Klitzing and J. H. Smet, Nano Lett., 2010, 10(11), 4544–4548.
66. A. K. Vallabhaneni, D. Singh, H. Bao, J. Murthy and X. Ruan, Phys. Rev. B: Condens. Matter, 2016, 93(12), 125432.
67. H. Cao, Q. Yu, L. A. Jauregui, J. Tian, W. Wu, Z. Liu, R. Jalilian, et al., 2009, arXiv preprint arXiv:0910.4329.
68. D. L. Nika, S. Ghosh, E. P. Pokatilov and A. A. Balandin, Appl. Phys. Lett., 2009, 94(20), 203103.
69. L. Lindsay and D. A. Broido, Phys. Rev. B: Condens. Matter, 2010, 81(20), 205441.
70. S. J. Plimpton, Comput. Phys., 1995, 117(1), 1–19.
71. F. Müller-Plathe, J. Chem. Phys., 1997, 106(14), 6082–6085.
72. G. Balasubramanian, I. K. Puri, M. C. Böhm and F. Leroy, Nanoscale, 2011, 3(9), 3714–3720.

---

Footnote

† Electronic supplementary information (ESI) available: Additional thermal conductivity measurements data. See DOI: 10.1039/c6nr03470e

---