Appl. Phys. A

DOI 10.1007/s00339-015-9489-1

Applied Physics A

Materials Science & Processing

Thermal conductivity reduction in graphene with silicon impurity

Byoung Seo Lee1 • Joon Sik Lee1

Received: 1 April 2015/Accepted: 11 September 2015

© The Author(s) 2015. This article is published with open access at Springerlink.com

Abstract We present a molecular dynamics investigation f Fraction

on the thermal conductivity of silicon-doped graphene and k Thermal conductivity (W/mK)

the resulting change in phonon properties. A significant m Mass

reduction in the thermal conductivity is observed in the L Length (nm)

presence of silicon impurity even at a small concentration l Length scale (nm)

of silicon atoms. Conductivity values continued to decrease q Heat flux (W/m2)

with an increase in silicon concentration. The increase in r Position vector

the scattering rate, which is measured by the reduction or r Interatomic distance (nm)

broadening of the peaks of the van Hove singularities, is T Temperature (K)

the most significant factor contributing to the large con- t Time (fs)

ductivity reduction. An analysis with scattering time v Speed or velocity (m/s)

models shows that the mass displaced by the silicon X, y, Z Space coordinates (nm)

impurity plays a significant role in reducing the conduc-

tivity, especially at a moderate concentration. The non- Greek symbols

mass effect, which comes from the change of the sp2 C-C c DOS calculated from MD

bonds to the sp3 Si-C bonds, is less strong or comparable D Increment

with the mass change effect. For high impurity concen- e Energy supplied (eV)

trations, the shape of the graphene is severely distorted and s Relaxation time

the irregularity of the ripples increases, which could con- / Pairwise potential

tribute to the reduction in conductivity. x Frequency

List of symbols

A Area (nm2) c Specific heat

cd Point defect concentration F Force vector

& Joon Sik Lee

jslee123@snu.ac.kr

Byoung Seo Lee elvis1@snu.ac.kr

1 School of Mechanical and Aerospace Engineering, Seoul

National University, Seoul 08826, Korea

Subscripts

d Dopant

eff Effective

i, j Indices for atoms

p Phonon

pd Point defect

oo Infinite

0 Host or pure reference

Abbreviations

DOS Density of state

EMD Equilibrium molecular dynamics

MD Molecular dynamics

NEMD Non-equilibrium molecular dynamics

Published online: 24 September 2015

a Springer

1 Introduction

Since graphene, a two-dimensional graphite, was fabricated by exfoliation from graphite [1], there has been significant research interest in graphene synthesis and graphene properties [2-4]. Numerous works have reported the remarkable properties of graphene, including extraordinary high electron mobility [5], thermal conductivity [6], stiffness and strength [7], as well as large surface area to volume ratio [8], and an unusual electronic structure [9]. However, some of the properties of graphene are not well suited to practical applications, which lead to intense research into new methods of functionalization or property tuning [10].

In applications where the thermal requirements are critical, it is important to regulate the change of the thermal properties. For graphene to be a high-performance thermoelectric (TE) candidate, it is desirable to lower its thermal conductivity as much as possible, which would increase its TE figure of merit [11]. However, it is often the case that the thermal transport properties of graphene are altered as a side effect of tuning other properties. Gra-phene's excellent heat transfer properties make it an ideal material for addressing challenges associated with the cooling loads of electric devices [6, 12]. Compared to carbon nanotubes (CNTs), graphene has lower contact resistance and superior thermal heat transfer properties. However, for practical applications, the zero bandgap of pristine graphene must be addressed with proper func-tionalization. The insertion of impurity or dopant atoms is the most effective method for addressing the bandgap concerns [10, 13]. On the other hand, it can negatively impact thermal conduction properties. Therefore, controlling the impact of impurity on heat transfer [14-16] becomes important for advancing graphene-based technology.

An understanding of the origin of graphene's high thermal conductivity is challenging because the models of two-dimensional and nanoscale physics are limited. The contribution of phonon scattering to the thermal conductivity is dominant over the electronic contribution [12]. Sadeghi et al. [17] have reviewed current challenges surrounding the thermal transport properties. As regards the heat transport of tuning graphene, investigations of isotope [18], nitrogen [16], and vacancy [19] have been performed. Recently, an elaborate model that corrects the Klemens' model, which was originally developed for general bulk materials and subsequently adopted for graphene, has been suggested [20].

The silicon-doped graphene has been recently synthesized [21], and density functional theory (DFT) calculation had shown the possibility of enhanced on/off efficiency [13]. It was proposed that bonding stability allows the

addition of different atoms bound to Si atoms [22]. In the present study, single-layered graphene with silicon as impurity has been investigated using non-equilibrium MD (NEMD) and equilibrium MD (EMD) to demonstrate the reduction in graphene's thermal conductivity and the relevant phonon behavior, respectively. MD has already been adopted in investigating the thermal conduction of gra-phene relating to rectification [14, 23] and thermal property management by controlling chirality, shape, and type of defect [14-16, 24-26]. To evaluate the scattering contributions in the Klemens' model from the MD results, Callaway type model calculations are conducted. Based on these calculations, additional consideration is given on the long-length behavior of the thermal conductivity.

2 Numerical methods

In classical MD, Newton's second law is applied to all atoms in the calculation domain, and the force Ft, acting on an atom i, is expressed as the gradient of potential energy according to

d/tj (r) drt '

where ri is the position vector of the atom i, and the intermolecular potential /jr) is given by Tersoff potential model [27, 28]. The potential is used for C-C, Si-Si, and C-Si interactions for silicon-doped graphene in a single form with multiple sets of parameters. All of the parameters used in the current study follow values previously tested and recommended in some applications by Tersoff [27]. The current programs have been validated through tests that reproduce previous results [29, 30]. The Newtonian equation is discretized in time using the Verlet algorithm with a fixed time step of At = 0.2 fs for both NEMD and EMD, and both methods have previously been sys-temically applied to silicon bonds with the Tersoff potential [31].

The procedure used by Wei et al. [32] is employed for NEMD calculations, where the half-sized domain, which is used for time-saving purposes, is the only difference from general NEMD calculations. The atomic arrangement of carbon and silicon is depicted in Fig. 1. The length of the model sheets ranges from 10.65 to 63.9 nm, and the width of all the simulated sheets is fixed at 3.9 nm. The silicon atoms are regularly distributed in convenience for the stability of the calculation and the consistent comparison of results for different sheet lengths. This artificial arrangement is not expected to have a significant impact on the results compared with using a random distribution. Five impurity concentrations, ranging between 0 and 10 %, are

Fig. 1 Atomic distribution for 5 % impurity with carbon atoms in white and silicon atoms in black

simulated. Periodic and free boundary conditions are applied in the y and z directions, respectively. The general procedures of NEMD (as described in [31]) are used in the simulations. The atoms start from their initial conditions and undergo the equilibration processes at a temperature of 300 K with 1,000,000 time steps. Then, the sheets are exposed to heat flux by adding and removing kinetic energy at the source and sink.

The thermal conductivity is calculated using Fourier's law from the resulting temperature gradient and the applied heat flux as

q = "k dX '

where heat flux, q, is applied as De/A not exceeding 23.0 x 10-5 eV/nm2 for each time step.

The same concentrations are used for both the EMD and NEMD calculations. The DOS calculated using EMD reveals the change of phonon property. The domain for the EMD calculations has x and y dimensions of 10.65 and 11.81 nm, respectively, and the periodic boundary conditions are extended to the x direction to simulate an infinite system. The Fourier transform over velocity autocorrelation function is given by

hP Vi(t)vi(0))

c(x) = •

<£ v2 (0)>

which is considered to be equivalent to the DOS [33].

3 Results and discussion

After a sufficient number of time steps to adapt to the heat flux, steady temperature profiles are established. Figure 2 presents a temperature profile obtained at a 63.9-nm-long pure graphene sheet subjected to a thermal energy of 15.3 x 10-5 eV/nm2. With the exception of the regions

Fig. 2 Temperature profile for the 63.9-nm-long pure graphene subjected to thermal energy of De/A = 15.3 x 10-5 eV/nm2

near the sink and source, a linear profile is generated, from which the thermal conductivity can be calculated.

The thermal conductivity values calculated for all conditions are shown in Fig. 3. The thermal conductivity increases with increasing sheet length, as expected because the phonon mean path depends on the sheet dimension. This can simply be explained using kinetic theory,

k = 3 CVP

1 1 2 l = 3 cv2

where c, vp, l, and s represent specific heat, phonon velocity, mean free path, and scattering time, respectively. As the impurity level of the graphene increases, the conductivity decreases.

Because of the extremely large computational cost, most of MD calculations employ an extrapolation scheme to estimate the bulk thermal conductivity based on smaller scale simulations, where Matthiessen's rule has been favored. In this case,

1/leff = 1/L + 1/li

where the free path of infinite size l? and the free path confined by distance between sink and source L contribute to the effective free path of finite length leff. According to the NEMD procedure, the effective conductivity approaches the bulk conductivity as 1/L approaches zero (i.e., at the infinite length). Thus, the inverse relationship shown in Fig. 3 facilitates the estimation of infinite conductivity k?. For the pristine graphene sheet, the extrapolation gives k? = 774 W/mK, which is in agreement with recent NEMD calculations on the same potential model [16, 32, 34]. The large deviation from linearity at the maximum concentration (10 % silicon atoms) is related to the

-1-1-r" impurity [%] -1-

—•—0.00

—0.63

—♦—1.25

_ —*—5.00

—*—10.0

0.00 0.02 0.04 0.06 0.08 0.10

ML [1/nm]

Fig. 3 Size-dependent thermal conductivities obtained from NEMD calculations for pure and impure graphenes

Fig. 4 Thermal conductivity ratios estimated by Eq. (5)

resolution of conductivity, i.e., the conductivity value is so small that the inverse value appears large.

The conductivity reduction ratios compared to pure infinite sheet, which is calculated from Eq. (5), are presented in Fig. 4. A drastic reduction in the thermal conductivity is observed. A 2 % impurity concentration results in more than an 80 % reduction, while a 10 % impurity concentration reduces the conductivity to <2 % of that in pure graphene. These reduction rates are larger than those obtained by simple mass substitution of isotopic atoms [18] or by nitrogen doping of graphene [16].

Using the NEMD model, an accurate fitting scheme is necessary to estimate the bulk thermal conductivity other than Eq. (5). There are three conjectures on the behavior of

the thermal conductivity: (1) As in MD simulation, the bulk thermal conductivity is usually estimated by the Matthiessen's rule preferred in general bulk (three dimensional) materials , which assumes a finite k at infinite length [31, 32]; (2) the k diverges in lower dimensions especially with logarithmic law (log dependence) for two dimensional lattices, which, in graphene, is known to be valid up to L < *9 im, and then diverges [35]; and (3) the k follows log dependence when L is smaller than about 100 im, and for longer L, it converges when the thermal transport reaches the diffusive regime [36].

Barbarino et al. [37] calculated the thermal conductivity for the graphene whose length is in the range of 0.83-100 im using approach-to-equilibrium MD (AEMD) with the reactive empirical bond order (REBO) potential, which can accommodate a large computational domain size with relatively shorter computing time. Their results support the third conjecture mentioned above.

We assumed Matthiessen's rule to estimate k?, which is turn out to be 774 W/mK. This is much smaller value compared to experimentally observed literature value of about 3000 W/mK [38, 39]. However, if we use the size dependence proposed by Barbarino et al. for fitting, instead of Matthiessen's rule, the thermal conductivity becomes about 2700 W/mK. The k/k (100 mm) at 60 nm is 0.1, and its value is 0.6 at about L = 3 mm where most experimental data are available; thus, 6 times of our k = 440 W/ mK at L = 60 nm yields about 2700 W/mK.

The DOS for the same configuration and concentration rate as in the NEMD cases is also obtained using EMD calculations. The calculated DOS of pure sheet is consistent with the dispersion curve obtained from Boltzmann transport calculations using the same original potential [40, 41]. Figure 5 presents the dispersion curve obtained from Boltzmann transport calculations and the DOS obtained from the EMD. Most of the singular points in DOS shown in Fig. 5 correspond to the frequencies of the symmetric points in the dispersion curve. Each of six phonon branches causes a high DOS peak at the M points. The nonzero and higher DOS at zero frequency (Fi_3 point) is related to the quadratic dispersion relation of out-of-plane phonon modes, which is a direct consequence of the two dimensionality of graphene [42]. At K1-2 where the optical and acoustic modes are contacted, the valley between the two M peaks is formed. This is analogous to the electronic band structure, where a zero bandgap is formed by the contact.

The DOS for the pristine and impure graphenes is presented in Fig. 6. To allow the DOS calculated for the different impurity concentrations to be compared, the frequency, x, has been corrected to adjust for the mass changes that result from the replacement carbon with silicon atoms [15], using

Xg/xd = л/fd (md/mG ) + (1 - fd ),

CO -Q JO

сл -Q

SL со О û

00 о Û

30 40 50 ^[GHz]

with subscripts 0 and d denoting the host and dopant fraction, respectively. Without the correction, the real confinement frequency is reduced from 73 GHz for the pure case to 70 GHz at 10 % impurity concentration.

The most remarkable change in the DOS distribution is that the increased concentration rate reduces or broadens all of the M point peaks of the pure sheet. In Adamyan's calculations [43] for sheets doped with aluminum atoms, which have a mass comparable to silicon atoms, the peaks are lower near the van Hove singularities. Generally, peak broadening indicates that the scattering rates increase [33, 44] for all phonons.

The contraction of the frequency regime, as shown in Eq. (6), suggests that the effect of the phonon velocity being reduced by the impurities can account for only a limited or small portion of the full conductivity reduction [see Eq. (4)]. Hence, the distribution of DOS indicates that the primary source of the reduction is the increased scattering rate. Both calculation models predict that the additional peaks arise at the same frequencies, K1-2 and K4-5. These additional peaks accompany the splitting of the intersected modes, which is analogous to the bandgap opening in the doped graphene [43, 45]. It has been suggested that the shapes of the additional peaks are the result of optical phonons that occupy a narrow band with relatively low phonon group velocity [45]. This conversion of peaks into modes that do not contribute to the heat transport could enhance the reduction in thermal conductivity. At the maximum concentration in Fig. 6e, the DOS is significantly different from the lower concentration models in that the peaks disappear and the gap separating the inplane acoustic modes from the other modes is bridged because of the growth of the additional distribution.

The cause of imperfect scattering due to defects in the crystalline structure can be analyzed with a few constituent elements in the Rayleigh scattering model, as shown by Klemens [46]. When adjusted to two-dimensional material, such as graphene [47], the phonon relaxation time is given by

24ncd x3 Spd = -

p, max

{sf + (S2 + S3 )2},

Fig. 6 DOS obtained from the present MD calculations with frequency corrected by Eq. (6). From top to bottom: 0.00, 0.63, 1.25, 5.00, and 10.0 % silicon concentrations

where cd is the concentration of defect atom. The elemental parameters S1, S2, and S3 are related to differences in mass, velocity, and radial spacing, respectively [46]. For example, S1 accounts for only mass changes, such as in isotopes. While the first term is the contribution of pure mass change, the latter two are commonly derived from the change of interatomic force. The scattering rate is larger

than the simple mass change when the change of the interatomic force involves different atomic species. Hence, the conductivity of graphene doped with heavy atoms such as silicon is lower than that of enriched isotope or of graphene doped with lighter atoms.

Figure 7 shows the instantaneous atomic arrangements. Even in the pure sheet, graphene has intrinsic ripples on the order of one nanometer, which is comparable to that of real suspended graphene [48]. The level of perturbation is getting higher than the natural ripple (Fig. 7a) (for pure graphene) with increasing silicon concentration. This perturbation occasionally results in irregular twisting of the smooth ripples. Non-planar shapes such as ripples and areas of high curvature can cause reduced thermal conductivity [49]; reductions by folding [50] and curvature [16] have been reported. While the large-scale irregular ripples occur when the impurity concentration is relatively high, as shown in Fig. 7, another source of strain induced by the change of bond type is present at even low impurity concentrations.

Figure 8 depicts tetrahedra comprised of four atoms. The tetrahedra comprised of four carbon atoms, i.e., those that do not include a Si impurity atom, are planar and show a propensity for sp2 bonding, as shown in Fig. 8a. The average bond angle and length are 120.0° and 0.146 nm, respectively, which are typical sp2 values. When the silicon impurity is at the center of the tetrahedra, with carbon neighbors, as shown in Fig. 8b, tetrahedra are no longer within two-dimensional plane, as the silicon atoms rise or sink away from the plane formed by the three adjacent carbon atoms. The average bond angle and length for every Si-C bond for all of the impure cases are 102.3° and 0.176 nm, respectively. These values indicate that the silicon impurity converts the two-dimensional sp2 bonds into sp3 bonds. The sp3 bond induces an increase in out-of-plane motion, and its disturbance to in-plane heat transfer is well known [51]. It is acknowledged that the addition of sp3 bonds reduces the in-plane thermal conductivity of graphene [52, 53]. Even a simple substitution of the bond type without the insertion of the impurity [52] can lead to a large reduction, of a similar order of magnitude to that which results from nitrogen doping [16], when sp2 orbitals are preserved in N-C bonds of the N-doped graphene. The thermal conductivity reduction due to either simple sp3 substitution or of lighter atomic substitution is commonly about 77 % at 5 % impurity concentration as in the references [16, 52]. The effect of sp replacement on pure strain (S3) is expected to be negligible for interstitial silicon compared with large mass substitution of a vacancy [54]. The second term, S2, which indirectly reflects changes in the intermolecular angle and distance due to the strain, for example, in Tersoff model, is expected to have a contribution comparable to the mass effect.

Fig. 7 Atomic arrangements of graphene with different concentrations. a 0.00 % Si, b 1.25 % Si, c 10.00 % Si

Callaway type calculations are performed using the same constituent models and calculation parameters as Alofi and Srivastav [47]. In Fig. 9, the solid symbol indicates the results obtained from the calculation and open symbol denotes the data obtained by the MD simulation. The conductivity values from the MD calculations are rescaled by half for comparison. The length dependence shown in Fig. 9a for pure graphene agrees with the recent reports that the thermal conductivity increases logarithmically but converges at a large (millimeter) scale [37]. The parameter S in the Klemens' imperfection model is

Fig. 8 Tetrahedra consisting of a central atom and its three carbon neighbor atoms. The central atom is a carbon and b silicon

determined from the thermal conductivity values of MD simulations. The mass parameter Si that results from the mass change is fixed because of the known mass of silicon, but non-mass terms, S23 = S2 + S3, are difficult to measure and are therefore calculated by fitting the thermal conductivity, which is obtained from the present MD calculation. In particular, at a moderate concentration of 0.63 %, the value of S23 = 0.15 gives the best fit to the thermal conductivity [a curve with open circle symbol in Fig. 9b]. While the mass parameter S1 is 0.38 at 0.63 % of silicon atoms, the S23 parameter is related to the change of intermolecular force and its resultant strain. The value

Fig. 9 Length dependence of thermal conductivity calculated using the Callaway type model (solid line). For impure graphene, S23 = 0.15 is used. a L < 100 im b L < 10 im with MD data (dashed line)

S23 = 0.15 is used for all of the impurity ratios. This choice of parameters is valid for the low concentrations, as shown in Fig. 9b, which demonstrates that the focused length scale is close to the range of MD calculation. At moderate concentrations, the mass effect (S1 = 0.38) is larger than or comparable to non-mass effect (S23 = 0.15). This is because the mass of the silicon atom is nearly two times larger than the mass of the host atom. In the Callaway type model calculation, the reduction ratios in long-length limit are increased compared with those for the small sizes. According to Eq. (5), the thermal conductivity at 0.63 % impurity is 28 % of the pristine graphene, as shown in Fig. 4; this ratio decreases to 22 % at 10 im in the Call-away model. The same tendency is found for the isotopic defect; some MD calculations [26, 55] based on the small

size have reported to underestimate the reduction in comparison with experimental and theoretical values [18, 36].

For higher concentrations, the value of S23 = 0.15 fitted for moderate concentration does not sufficiently account for the reduction in the thermal conductivity. The predicted thermal conductivity is smaller than that predicted using MD, and the deficiency in S is not negligible. A new kind of scattering model is needed to describe the effect of the irregular ripple, as shown in Fig. 7c. The strain evaluated for the sp3 bond follows the Klemens' scattering formula, as shown in Fig. 7b. In this case, every defect has the same strain rate and interatomic relation, so the defect density is proportional to the scattering rate, which the Klemens' model also implies. However, the ripple, which is not subject to a single point defect, requires an understanding of the principle of formation and the resultant effect on thermal transport.

4 Summary and conclusions

Molecular dynamics simulations are performed to investigate the effect of silicon impurity on the thermal conductivity of graphene sheets. The results obtained are in good agreement with previously published phonon property calculations. NEMD reveals a drastic reduction in thermal conductivity with even a small concentration of silicon impurity. The calculations show that only 0.63 % atomic replacement decreases the thermal conductivity to 30 % of that of pure graphene. Impurity concentrations of 3 % reduce the conductivity to less than one-tenth of the pure graphene conductivity. These results suggest that silicon impurity is more effective at reducing the conductivity than the isotopic dopant; experiments have shown that the thermal conductivity of graphene isotopi-cally enriched with 1.1 % of 13C is reduced to 63 % of pure graphene [18]. Further, our results suggest that silicon impurity is superior to nitrogen doping, which has been observed to reduce the thermal conductivity to 30 % of pure graphene at of 3 % nitrogen impurity [16]. For higher silicon concentrations, the present NEMD shows that the conductivity decreases to 1.7 % of the pure gra-phene conductivity at 10 % silicon. The reduction rate may increase when the logarithmic length dependence is applied, which has been recently supported by experiments and calculations.

The DOS obtained from EMD calculations shows the phonon properties expected in the presence of crystalline imperfectness. The peak broadening or reduction becomes more significant at higher concentration, which is related directly to the reduction in phonon heat transfer. The broadening implies an increase in the scattering rate, which results in a reduction in thermal conductivity. The increase

in scattering by imperfection is shared with all of the acoustic and optical branches of phonons.

From the scattering time model analysis, it can be seen that the contribution of the mass element (S1) of the silicon impurity plays a significant role in reducing the conductivity, especially at a moderate concentration. The non-mass effect comes from the change of the sp2 C-C bonds to the sp3 Si-C bonds. For high impurity concentrations, the shape of the graphene is severely distorted and the irregularity of the ripples increases, which could contribute to the reduction in conductivity. We propose that the ripples have a different scattering source than that which is normally modeled in Klemens' formula.

Acknowledgments This research work was supported partially by the World Class University (WCU) Program (R31-2008-000-10083-0) and partially by the Nano-Material Technology Development Program (R2011-003-2009), both through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT and Future Planning. Additionally, it was supported by the National Institute of Supercomputing and Network/Korea Institute of Science and Technology Information with supercomputing resources including technical support (KSC-2013-C01-034).

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

References

1. K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Electric field effect in atomically thin carbon films. Science 306(5696), 666-669 (2004). doi: 10.1126/science.1102896

2. C. Berger, Z. Song, X. Li, X. Wu, N. Brown, C. Naud, D. Mayou, T. Li, J. Hass, A.N. Marchenkov, E.H. Conrad, P.N. First, W.A. de Heer, Electronic confinement and coherence in patterned epitaxial graphene. Science 312(5777), 1191-1196 (2006). doi: 10.1126/science.1125925

3. K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. Khotke-vich, S.V. Morozov, A.K. Geim, Two-dimensional atomic crystals. Proc. Natl. Acad. Sci. USA 102(30), 10451-10453 (2005). doi:10.1073/pnas.0502848102

4. L. Jiao, L. Zhang, X. Wang, G. Diankov, H. Dai, Narrow gra-phene nanoribbons from carbon nanotubes. Nature 458(7240), 877-880 (2009). doi:10.1038/nature07919

5. K.I. Bolotin, K.J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, H.L. Stormer, Ultrahigh electron mobility in suspended graphene. Solid State Commun. 146(9-10), 351-355 (2008). doi:10.1016/j.ssc.2008.02.024

6. A.A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, C.N. Lau, Superior thermal conductivity of single-layer graphene. Nano Lett. 8(3), 902-907 (2008). doi:10.1021/ nl0731872

7. C. Lee, X. Wei, J.W. Kysar, J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science 321(5887), 385-388 (2008). doi:10.1126/science.1157996

8. M.D. Stoller, S. Park, Y. Zhu, J. An, R.S. Ruoff, Graphene-based ultracapacitors. Nano Lett. 8(10), 3498-3502 (2008). doi:10. 1021/nl802558y

9. Y. Zhang, Y.-W. Tan, H.L. Stormer, P. Kim, Experimental observation of the quantum hall effect and Berry's phase in graphene. Nature 438(7065), 201-204 (2005). doi:10.1038/nature04235

10. V. Georgakilas, M. Otyepka, A.B. Bourlinos, V. Chandra, N. Kim, K.C. Kemp, P. Hobza, R. Zboril, K.S. Kim, Functional-ization of graphene: covalent and non-covalent approaches, derivatives and applications. Chem. Rev. 112(11), 6156-6214 (2012). doi:10.1021/cr3000412

11. R. Verma, S. Bhattacharya, S. Mahapatra, Thermoelectric performance of a single-layer graphene sheet for energy harvesting. IEEE Trans. Electron Devices 60(6), 2064-2070 (2013)

12. S. Ghosh, I. Calizo, D. Teweldebrhan, E.P. Pokatilov, D.L. Nika, A.A. Balandin, W. Bao, F. Miao, C.N. Lau, Extremely high thermal conductivity of graphene: prospects for thermal management applications in nanoelectronic circuits. Appl. Phys. Lett. 92(15), 151911 (2008). doi:10.1063/1.2907977

13. M.S.S. Azadeh, A. Kokabi, M. Hosseini, M. Fardmanesh, Tunable bandgap opening in the proposed structure of silicon-doped graphene. Micro Nano Lett. 6(8), 582-585 (2011)

14. Y. Wang, S. Chen, X. Ruan, Tunable thermal rectification in graphene nanoribbons through defect engineering: a molecular dynamics study. Appl. Phys. Lett. 100(16), 163101 (2012). doi:10.1063/1.3703756

15. G. Balasubramanian, I. K. Puri, M. C. Bohm, and F. Leroy, Thermal conductivity reduction through isotope substitution in nanomaterials: predictions from an analytical classical model and nonequilibrium molecular dynamics simulations. Nanoscale. 3(9), 3714-3720, 045413 (2011). doi:10.1039/c1nr10421g

16. B. Mortazavia, A. Rajabpour, S. Ahzi, Y. Remonda, S.M.V. Allaeid, Nitrogen doping and curvature effects on thermal conductivity of graphene: a non-equilibrium molecular dynamics study. Solid State Commun. 152(4), 261-264 (2012). doi:10. 1016/j.ssc.2011.11.035

17. M.M. Sadeghi, M.T. Pettes, L. Shi, Thermal transport in gra-phene. Solid State Commun. 152(15), 1321-1330 (2012). doi:10. 1016/j.ssc.2012.04.022

18. S. Chen, W. Qingzhi, C. Mishra, J. Kang, H. Zhang, K. Cho, W. Cai, A.A. Balandin, R.S. Ruoff, Thermal conductivity of iso-topically modified graphene. Nat. Mater. 10(3), 203-207 (2012). doi:10.1038/nmat3207

19. H. Zhang, G. Lee, K. Cho, Thermal transport in graphene and effects of vacancy defects. Phys. Rev. B 84(11), 115460 (2011). doi:10.1103/PhysRevB.84.115460

20. G. Xie, Y. Shen, X. Wei, L. Yang, H. Xiao, J. Zhong, G. Zhang, A bond-order theory on the phonon scattering by vacancies in two-dimensional materials. Sci. Rep. 4, 5085 (2014). doi:10. 1038/srep05085

21. R. Lv, M.C. dos Santos, C. Antonelli, S. Feng, K. Fujisawa, A. Berkdemir, R. Cruz-Silva, A.L. Elias, N. Perea-Lopez, F. Lopez-Urias, H. Terrones, M. Terrones, Large-area Si-doped graphene: controllable synthesis and enhanced molecular sensing. Adv. Mater. 26(45), 7593-7599 (2014). doi:10.1002/adma.201403537

22. S.B. Fagan, A.J.R. da Silva, R. Mota, R.J. Baierle, A. Fazzio, Functionalization of carbon nanotubes through the chemical binding of atoms and molecules. Phys. Rev. B 67(3), 033405 (2003). doi: 10.1103/PhysRevB.67.033405

23. Q.-X. Pei, Y.-W. Zhang, Z.-D. Sha, V.B. Shenoy, Carbon isotope doping induced interfacial thermal resistance and thermal rectification in graphene. Appl. Phys. Lett. 100(10), 101901 (2012). doi:10.1063/1.3692173

24. S.-K. Chien, Y.-T. Yang, C.-K. Chen, The effects of vacancy defects and nitrogen doping on the thermal conductivity of armchair (10, 10) single-wall carbon nanotubes. Solid State

Commun. 151(14-15), 1004-1008 (2010). doi:10.1016/j.ssc. 2011.04.025

25. G. Zhang, B. Li, Thermal conductivity of nanotubes revisited: effects of chirality, isotope impurity, tube length, and temperature. J. Chem. Phys. 123(11), 114714 (2005). doi:10.1063/1. 2036967

26. J. Hu, S. Schiffli, A. Vallabhaneni, X. Ruan, Y.P. Chen, Tuning the thermal conductivity of graphene nanoribbons by edge passivation and isotope engineering: a molecular dynamics study. Appl. Phys. Lett. 97(13), 133107 (2011). doi:10.1063/1.3491267

27. J. Tersoff, Modeling solid-state chemistry: interatomic potentials for multicomponent systems. Phys. Rev. B 39(8), 5566-5568 (1989). doi:10.1103/PhysRevB.39.5566

28. J. Tersoff, Empirical interatomic potential for carbon, with applications to amorphous carbon. Phys. Rev. Lett. 61(25), 2879-2882 (1988). doi:10.1103/PhysRevLett.61.2879

29. T. Halicioglu, Comparative study on energy- and structure-related properties for the (100) surface of b-SiC. Phys. Rev. B 51(11), 7217-7223 (1995). doi:10.1103/PhysRevB.51.7217

30. M. Ishimaru, K. Yoshida, T. Kumamoto, T. Motooka, Molecular-dynamics study on atomistic structures of liquid silicon. Phys. Rev. B 54(7), 2641-4638 (1996). doi:10.1103/PhysRevB.54.4638

31. P.K. Schelling, S.R. Phillpot, P. Keblinski, Comparison of atomic-level simulation methods for computing thermal conductivity. Phys. Rev. B 65(14), 144306 (2002). doi:10.1103/ PhysRevB.65.144306

32. Z. Wei, Z. Ni, K. Bi, M. Chen, Y. Chen, In-plane lattice thermal conductivities of multilayer graphene films. Carbon 49(8), 2653-2658 (2011). doi:10.1016/j.carbon.2011.02.051

33. J.M. Dickey, A. Paskin, Computer simulation of the lattice dynamics of solids. Phys. Rev. 188(3), 1407-1418 (1969). doi:10. 1103/PhysRev.188.1407

34. A. Cao, Molecular dynamics simulation study on heat transport in monolayer graphene sheet with various geometries. J. Appl. Phys. 111(8), 083528 (2012). doi:10.1063/1.4705510

35. X. Xu, L.F.C. Pereira, Y. Wang, J. Wu, K. Zhang, X. Zhao, S. Bae, C. Tinh Bui, R. Xie, J.T.L. Thong, B.H. Hong, K.P. Loh, D. Donadio, B. Li, B. Ozyilmaz, Length-dependent thermal conductivity in suspended single-layer graphene. Nat. Commun. 5, 3689 (2014). doi: 10.1038/ncomms4689

36. G. Fugallo, A. Cepellotti, L. Paulatto, M. Lazzeri, N. Marzari, F. Mauri, Thermal conductivity of graphene and graphite: collective excitations and mean free paths. Nano Lett. 14(11), 6109-6114

(2014). doi:10.1021/nl502059f

37. G. Barbarino, C. Melis, L. Colombo, Intrinsic thermal conductivity in monolayer graphene is ultimately upper limited: a direct estimation by atomistic simulations. Phys. Rev. B 91(3), 035416

(2015). doi:10.1103/PhysRevB.91.035416

38. J.-U. Lee, D. Yoon, H. Kim, S.W. Lee, H. Cheong, Thermal conductivity of suspended pristine graphene measured by raman spectroscopy. Phys. Rev. B 83(8), 081419 (2011). doi:10.1103/ PhysRevB.83.081419

39. J.H. Seol, I. Jo, A.L. Moore, L. Lindsay, Z.H. Aitken, M.T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R.S. Ruoff, L. Shi, Two-dimensional phonon transport in supported graphene. Science 328(5975), 213-216 (2010). doi:10.1126/science. 1184014

40. L. Lindsay, D.A. Broido, Optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene. Phys. Rev. B 81(20), 205441 (2010). doi:10.1103/PhysRevB.81.205441

41. L. Lindsay, D.A. Broido, Erratum: optimized Tersoff and Brenner empirical potential parameters for lattice dynamics and phonon thermal transport in carbon nanotubes and graphene. Phys. Rev. B 82(20), 209903(E) (2010). doi:10.1103/PhysRevB. 82.209903

50. N. Yang, X. Ni, J.-W. Jiang, B. Li, K.S. Novoselov, T.J. Booth, S. Roth, How does folding modulate thermal conductivity of graphene? Appl. Phys. Lett. 100(9), 093107 (2012). doi:10.1063/1. 3690871

51. D. Yang, F. Ma, Y. Sun, T. Hu, K. Xu, Influence of typical defects on thermal conductivity of graphene nanoribbons: an equilibrium molecular dynamics simulation. Appl. Surf. Sci. 258(24), 9926-9931 (2012). doi:10.1016/j.apsusc.2012.06.052

52. A. Rajabpour, S.M.V. Allaei, Tuning thermal conductivity of bilayer graphene by inter-layer sp3 bonding: a molecular dynamics study. Appl. Phys. Lett. 101(5), 053115 (2012). doi:10. 1063/1.4740259

53. B. Liu, C.D. Reddy, J. Jiang, J.A. Baimova, S.V. Dmitriev, A.A. Nazarov, K. Zhou, Morphology and in-plane thermal conductivity of hybrid graphene sheets. Appl. Phys. Lett. 101(21), 211909 (2012). doi:10.1063/1.4767388

54. C.A. Ratsifaritana, P.G. Klemens, Scattering of phonons by vacancies. Int. J. Thermophys. 8(6), 737-750 (1987). doi:10. 1007/BF00500791

55. X. Li, J. Chen, C. Yu, G. Zhang, Comparison of isotope effects on thermal conductivity of graphene nanoribbons and carbon nan-otubes. Appl. Phys. Lett. 103(1), 013111 (2013). doi:10.1063/1. 4813111

doi:10.1103/PhysRevB.82.115427

42. R. Saito, G. Dresselhaus, and M. S. Dresselhaus, Physical Properties of Carbon Nanotubes, (Imperial College Press, Singapore, 1998), 093107, pp. 165-170

43. V. Adamyan, V. Zavalniuk, Phonons in graphene with point defects. J. Phys. Condens. Matter 23(1), 015402 (2011). doi:10. 1088/0953-8984/23/1/015402

44. H. Frase, B. Fultz, J.L. Robertson, Phonons in nanocrystalline Ni3Fe. Phys. Rev. B 57(2), 898-905 (1998). doi:10.1103/Phys RevB.57.898

45. S.K. Gupta, H.R. Soni, P.K. Jha, Electronic and phonon band-structures of pristine few layer and metal doped graphene using first principles calculations. AIP Adv. 3(3), 032117 (2013). doi:10.1063/1.4794949

46. P.G. Klemens, The scattering of low-frequency lattice waves by static imperfections. Proc. Phys. Soc. A 68(712), 1113-1128 (1955). doi:10.1088/0370-1298/68/12/303

47. A. Alofi, G.P. Srivastava, Thermal conductivity of graphene and graphite. Phys. Rev. 87(11), 115421 (2013). doi:10.1103/Phys RevB.87.115421

48. J.C. Meyer, A.K. Geim, M.I. Katsnelson, K.S. Novoselov, T.J. Booth, S. Roth, The structure of suspended graphene sheets. Nature 446(7131), 60-63 (2007). doi:10.1038/nature05545

49. L. Lindsay, D.A. Broido, N. Mingo, Flexural phonons and thermal transport in graphene. Phys. Rev. B 82(11), 115427 (2010).