Tailoring Thermal Transport Property of Graphene through Oxygen

Dec 30, 2013 - We compute thermal conductivity of graphene oxide at room temperature with molecular dynamics simulation. To validate our simulation mo...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

Tailoring Thermal Transport Property of Graphene through Oxygen Functionalization Hengji Zhang,† Alexandre F. Fonseca,‡,§ and Kyeongjae Cho*,† †

Department of Materials Science and Engineering and Department of Physics, The University of Texas at Dallas, Richardson, Texas 75080, United States ‡ Department of Physics, UNESP−Sao Paulo State University, Bauru, SP 17033-360, Brazil § Instituto de Física “Gleb Wataghin”, Universidade Estadual de Campinas, UNICAMP, Campinas, SP 13083-859, Brazil

ABSTRACT: We compute thermal conductivity of graphene oxide at room temperature with molecular dynamics simulation. To validate our simulation model, we have investigated phonon scattering in graphene due to crystal boundary length and isotope defect, both of which are able to diagnose the behavior of long wavelength and short wavelength phonon scattering. Our simulation shows that thermal conductivity of pristine graphene has logarithmic divergence for the boundary length up to 2 μm. As compared with pristine graphene, thermal conductivity of graphene oxide can be reduced by a factor of 25 at low oxygen defect concentration. Moreover, we find that not only the concentration but also the configuration of the oxygen functional groups (e.g., hydroxyl, epoxide, and ether) has significant influence on the thermal conductivity. Through phonon mode analysis, phonon defect scattering as well as phonon localization are mainly responsible for the conspicuous reduced thermal conductivity. The simulation results have provided fundamental insight on how to precisely control thermal property of graphene oxide for thermal management and thermoelectric applications.

1. INTRODUCTION Single layer graphene has attracted a great deal of attention because of its ultrahigh thermal conductivity,1,2 electron mobility,3 and mechanical strength.4 Despite all these attractive properties, zero band gap and limited-scale production of graphene are the major challenges that restrict its device applications. To overcome these challenges, extensive studies have been devoted on band gap opening through engineering its structure into graphene nanoribbon,5 hydrogenated graphene,6 or graphene oxide (GO).7−9 Meanwhile, mass production of graphene can be achieved by reduction of GO since it can be readily converted into graphene on large-scale through chemical or thermal reduction treatment.8,10,11 Oxygen functionalized graphene therefore has been recognized as an interesting chemically modified carbon material because it not only provides a possible solution to the aforementioned challenges in electronic application but also possesses thermal rectification property12 that may advance diverse thermal management such as electronics cooling, thermal diode, and thermal logic circuits.13−15 More than these, in recent © 2013 American Chemical Society

experimental work, oxygen functionalized graphene also has shown interesting thermoelectric properties.16 The efficiency of a thermoelectric material is determined by its figure of merit, ZT = S2σT/κ, where S is the Seebeck coefficient, σ is the electrical conductivity, T is the absolute temperature, and κ is the thermal conductivity consisting of phonon and electron heat conduction. According to recent experimental works on graphene, carbon nanotube composites, and graphite intercalation compounds, ZT values of these carbon materials are in the range of 10−3 to 10−2.16,17 In order to further enhance ZT value up to that of useful thermoelectric materials, the ratio of thermopower factor (S2σ) to thermal conductivity (κ) needs to be increased by 2−3 orders of magnitude. The recent experiment by Xiao et al.16 has demonstrated that thermopower factor of few-layer graphene after oxygen plasma treatment can increase by a factor of 15. In Received: September 27, 2013 Revised: December 26, 2013 Published: December 30, 2013 1436

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

Article

such an oxidative environment, oxidation of graphene would inevitably occur, and the oxygen defects can also decrease thermal conductivity due to the enhanced phonon scattering by defects. Therefore, intentionally introducing oxygen defects in graphene seems to be a promising route to tailor thermoelectric property of graphene-based materials. Considering all these potential applications of GO on thermal management and thermoelectric energy conversion, it is imperative to determine thermal conductivity of oxygen functionalized graphene. Thermal transport study usually can be carried out with lattice dynamics method,18,19 Boltzmann transport equation (BTE) theory,20−22 or classical molecular dynamics (MD) simulation.15,23 On the subject of graphene with defects, classical MD simulation has provided an effective approach to calculate thermal conductivity if a reliable and accurate interatomic potential is able to describe intrinsic and extrinsic phonon scatterings.24 Recently, the effects of various defects on the thermal transport property of graphene (e.g., vacancy,24,25 isotope,26 grain boundary,27,28 hydrogen,29−31 and fluorine32) have been systematically studied. Nevertheless, thermal transport modeling of GO with MD simulation is not attainable until an improved reactive empirical bond order potential of carbon/hydrogen/oxygen (REBO−CHO) has been developed.33 The REBO potential can simulate chemical bond breakage and formation by evaluating environmentally dependent bond orders.34 Compared with other well-known CHO potentials (e.g., ReaxFF-CHO35), the unique feature of reoptimized REBO−CHO is that it can describe atomic interactions for a micrometer sized GO flake with good accuracy and low computational expense. Thus, REBO−CHO is a more suitable candidate for thermal transport modeling. According to recent experiments, κ of GO can range from 2 to 1000 W/mK depending on the different oxygen reduction methods.36,37 To understand the impacts of oxygen functionalization, our simulation study would provide fundamental insight on phonon transport behavior for GO. In this article, we have performed nonequilibrium molecular dynamics simulations to investigate thermal transport property of oxygen functionalized graphene. In the following, we first calculated the κ of graphene as functions of crystal boundary length (L) and isotope defect. One notable finding is that κ of pristine graphene diverges and scales linearly as a logarithmic function of L (i.e., κ ∝ ln(L)). This result is consistent with the earlier studies on two-dimensional systems.38−41 As the oxygen functional groups (e.g., hydroxyl, epoxide, and ether) are introduced in graphene, κ is independent of crystal boundary length. Compared with pristine graphene of the same L, the reduction ratio of κ for GO can range from 3 to 25 at low oxygen concentration. Recalling the increase of S2σ by a factor of 15, it is promising that large enhancement for ZT value of graphene could be attainable through oxygen functionalization. In the last part, we conclude that phonon scattering and phonon localization effect due to oxygen functionalization are the major factors that contribute to the reduction of κ. The first order effect is due to phonon scattering, which strongly depends on the concentration of oxygen defects. The second order effect is due to phonon localization, which depends on the specific configuration of oxygen defects.

Figure 1. (a) G−OH. (b) G−Epoxide. (c) G−Ether. The carbon, hydrogen, and oxygen (a−c) are colored with black, blue, and red, respectively. (d) Schematic diagram of thermal conductivity calculation with nonequilibrium MD simulation. (e) Linear fit to temperature profile for pristine graphene with L = 300 nm at 300 K. The origin of the position axis in panel e corresponds to the midpoint position in panel d.

periodic boundary condition with W = 10 nm is applied in Y direction. We have confirmed that this size of W is large enough to have converged κ. Both the heat source and the heat sink have a size of L′ = 20 nm. The middle region between the heat source and the heat sink is defined as the boundary length of L. In order to run simulation for boundary length up to micrometer scales, a parallel version code for REBO−CHO potential is used in the study. The bond order term in REBO− CHO potential has been improved such that the binding energies for graphene with different oxygen functional groups are consistent with the results from density function theory.33 Meanwhile, for C−C interactions, we use a reoptimized REBO carbon potential because it can predict thermal conductivity of pristine graphene in close agreement with experiments due to the improved vibration frequencies for long wavelength phonons.24,42 In the simulation, the system is fully relaxed at 300 K using NVT ensemble. After reaching equilibrium, the system is switched to NVE ensemble. Then, according to the algorithm developed by Ikeshoji and Hafskjold,43 a constant amount of kinetic energy ΔE is injected into the heat source and extracted out of the heat sink in each simulation step Δt = 0.2 fs. Therefore, this can impose a constant heat current flow as J = ΔE/A/Δt in the middle region, which has a cross sectional area of A = W(3.4 Å). As the system has reached steady states, the temperature profile is linear as shown in Figure 1e, and a statistical average of temperature gradient can be used to calculate thermal conductivity with Fourier law expressed as k = J/Grad(T).

2. THEORETICAL METHODS MD Simulation Methods. The thermal transport model is illustrated by a schematic diagram in Figure 1d. As it shows, the atoms at both ends of graphene in X direction are fixed, and a 1437

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

Article

Phonon Localization Analysis. The phonon localization effect can be descripted with phonon participation ratio Pλ for each phonon mode, which can be expressed as44,45 N

Pλ−1 = N ∑ (∑ εi , α , λεi , α , λ)2 i=1

α

(1)

where εi,α,λ is the vibration amplitude of atom i along the Cartesian direction α for the mode λ. N is the total number of atoms. This mode-dependent participation ratio provides a way to compute the fraction of atoms participating in the mode λ, and its value is between O(1) for delocalized states and O(1/ N) for localized states. For two-dimensional graphene material, both the in-plane (i.e., longitudinal and transverse modes) and out-of-plane (i.e., flexural modes) acoustic modes are major heat carriers. Through calculating Pλ for the in-plane and outof-plane phonon modes, the phonon localization effect can be well characterized.

Figure 2. Thermal conductivity of graphene, as a function of boundary length, is calculated with MD simulation and Klemens model (KM). Solid square/triangle represents the data for pristine graphene. Empty square/triangle represents the data for graphene with 50% randomly mixed isotope C13. The inset figure shows the reference data using BTE theory.20 The coefficients of linear fit for different models are also given.

3. RESULTS AND DISCUSSION The oxygen functionalized graphene considered in the study are mainly categorized into three types of GO: graphene− hydroxyl (G−OH), graphene−epoxide (G−Epoxide), and graphene−ether (G−Ether). These major functional groups are hydroxyl and epoxide on the basal plane, and the ether group is attached to the edges.10,46,47 Because of nonstoichiometry and inhomogeneous nature of GO, these oxygen functional groups are randomly distributed on both sides of graphene surface. The top views of their atomic structures are shown in Figure 1a−c. Such a classification helps to understand how thermal conductivity of pristine graphene can be affected by different groups of oxygen defects. Recently, GO with high chemical homogeneity have been made in experiment through sequential removal of chemical functional groups.8 Also, thermally reversible oxidation of graphene can be made under ultrahigh vacuum conditions.9 Hence, the materials considered in the simulation can be practically synthesized. Before looking into the effect of different oxygen functional groups, we first investigated the boundary size and the isotope effects on the κ of graphene. In principle, boundary size plays an important role in long wavelength acoustic phonon scattering, whereas isotope defect has a relatively larger influence on short wavelength acoustic phonon scattering. Through studying these two effects, we can affirm the quality of the simulation model by comparing it with published theory and experiment results. The simulation result can also help to clarify the behavior of long wavelength and short wavelength phonons in pristine graphene, and it provides a useful reference before studying GO. On the boundary size effect, previous classical MD simulations on pristine graphene showed that κ is converged.23,24,27,48,49 BTE method is either converged50,51 or diverged.41,52 As we used reoptimized REBO carbon potential, we examined the boundary size effect on κ for graphene. Our simulation result in Figure 2 shows that κ of pristine graphene scales linearly as a function of ln(L). The similar linear logarithm relation has been predicted by earlier analytic and numerical studies for phonon transport in two-dimensional system.38−41 The data shown in the inset of Figure 2 also supports that linear logarithm relation with BTE method.20 Furthermore, we have calculated κ of graphene with Klemens model,53 and the same linear logarithm relation is demonstrated in Figure 2. Such comparison between classical MD

simulation and BTE method helps to clarify that our simulation model can accurately describe the phonon−phonon and phonon−boundary scattering processes for graphene. In addition, as graphene has high Debye temperature (∼2000 K) and quantum effects are not negligible at room temperature,18,50 it is not clear on whether classical MD simulation is able to accurately predict thermal conductivity. The BTE method predicts that κ of graphene is 3600 W/mK (10 μm in length). Using the reoptimized REBO potential, we can determine κ for graphene using the equation (κ = 526.8 ln(L) − 1590.9) shown in Figure 2. The κ calculated from our simulation is 3261 W/mK. The difference between quantum and classical MD method does not exceed 10%, which is in the range of calculation error bar of classical MD simulation. Therefore, classical MD simulation is a reasonable accurate method to predict κ for graphene at room temperature. On the isotope effect, we have calculated thermal conductivity of graphene with 50% C13 isotope with MD simulation and Klemens model. In the previous modeling work, we used Green−Kubo method to calculate thermal conductivity of graphene under periodic boundary conditions.54 The limitation of Green−Kubo method is the lack of long wavelength phonons and the boundary scattering effect. Hence, it can overestimate the reduction ratio for κ due to the lack of long wavelength phonons, which are mostly scattered at boundary rather than by isotope point defects. To address this problem, long wavelength phonon corrections based on Klemens model was applied in order to model sample length in experimental conditions.26 In this study, nonequilibrium MD simulation itself includes the boundary scattering effect, and we found that reduction ratio for κ predicted from our simulation is consistent with Klemens model as well as the experimental data.26 Moreover, both MD simulation and Klemens model imply that isotope impurities cannot effectively scatter long wavelength phonons, and κ keeps increasing as the boundary length increases. As for GO, hydroxyl and epoxide are commonly observed oxygen functional groups.8−11 During oxygen reduction process, these oxygen functional groups are mostly removed, whereas ether group as shown in Figure 1c can become a major 1438

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

Article

Figure 3. (a) Thermal conductivity of G−OH, G−Epoxide, and G−Ether as a function of O/C ratio at 300 K. The solid and dashed curve lines are fitted with the results from MD simulation and Klemens model, respectively. (b) Thermal conductivities of G−OH, G−Epoxide, and G−Ether for O/C ratio of 0.01 with different boundary length. (c) In-plane (κxy), out-of-plane (κz), and total (κxyz) thermal conductivities of G−Ether as a function of correlation time τ. κsum is the summation of κxy and κz. The O/C ratio for G−Ether is 0.04.

defect structure.47 Herein, we mainly focused on these three oxygen functional groups with random distribution on the graphene surface. With Klemens model and MD simulation, κ is computed as a function of O/C ratio for three types of GO as shown in Figure 3a. The boundary size of GO in this calculation is 300 nm. For the same O/C ratio, Klemens model predicts that both G−Epoxide and G−OH have almost the same κ, but G−Ether has a much lower κ. In MD simulation, κ(G−Ether) is still the lowest; however, in contrast to Klemens model, κ(G−Epoxide) is significantly lower than κ(G−OH). We will explain this notable difference between the two modeling methods in the following phonon mode analysis. It is worthwhile to mention that Klemens model has neglected the flexural phonons because of the zero group velocity in the long wavelength limit. According to recent studies, it is found that intermediate and short wavelength flexural phonons make a sizable contribution to κ of pristine graphene.24,52 In the case of GO, it would be interesting to know how the oxygen defects affect κ of flexural phonons. Using the freezing method,24 our calculated in-plane κxy and out-of-plane κz for G−Ether are shown in Figure 3c. The summation of κxy and κz is in good agreement with κxyz. This result indicates that freezing method is a reasonable approximation to calculate κxy and κz. It is, then, evident that the intermediate and short wavelength flexural phonons are strongly scattered by oxygen defects. Therefore, the flexural phonon contribution to κ is largely suppressed (κz = ∼5 W/mK), and Klemens model is still accurate to study heat transfer in graphene with defects.

As compared with pristine graphene of the same boundary length (κ = 1413 W/mK), we find that κ of oxygen functionalized graphene can be reduced by a factor ranging from 3 to 25, depending on the O/C ratio and the configuration of oxygen functional group. On the basis of the simulation result, hundreds times of ZT enhancement for oxygen functionalized graphene might be achieved because of the largely reduced κ as shown in the simulation and the 15 times enhancement for S2σ shown in earlier experiment.16 On thermal management application, we find that reduced GO with O/C ratio smaller than 0.01 can be applicable for device heat dissipation as κ ranges from 200 to 400 W/mK, which is as high as thermal conductivity of silver and copper. In Figure 3a, we also fit the MD simulation data with an inverse power law relation, which was used to predict κ at different point defect concentrations.55 According to this fitting curve, thermal conductivity of GO in a wide range of O/C ratios has been predicted. At the same O/C ratio of 0.10, κ(G−Epoxide) is ∼35.36 W/mK and κ(G−OH) is ∼85 W/mK. The hydroxyl group has shown greater impact on phonon transport than the epoxide group. A similar effect for electron transport in GO is also observed in experiment.56 In Figure 3b, we computed κ of GO as a function of boundary length. Since phonons are mainly scattered by oxygen functional groups, phonon transport in GO are more diffusive due to multiple scatterings with increased defect concentration. In contrast with isotope effect as shown in Figure 2, κ of GO is converged and independent of boundary length. This simulation demonstrates that oxygen defects in graphene are 1439

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

Article

additional phonon peaks at frequencies about 5 and 120 THz. The 5 THz peak corresponds to C−OH torsional vibrations, and the 120 THz peak corresponds to O−H bond vibrations. By calculating the partial DOS for all carbon atoms, we found that carbon atoms do not contribute for those vibration peaks. According to kinetic theory, κ is a function of heat capacity, phonon group velocity, and phonon mean free path. Our calculation suggests that there is a minor change on heat capacity and phonon group velocity because of the very little change on phonon peak position and phonon DOS value. Hence, the reduced κ is mainly caused by phonon scattering by defects, which result in the reduction of phonon mean free path. Using Klemens model, we can discuss about phonon scattering rate due to point defect. The phonon scattering rate (1/τ) due to point defect in Klemens model is given as 1/τ ∝ (ΔM/M)2, where M is the average atomic mass and ΔM is the mass difference of a single substitution defect in a crystal.58 In the case of C13 isotope defect, ΔM/M = 1/12. In an approximate way, the phonon scattering rate for G−Epoxide (G−OH) can be estimated as ΔM/M = 16/12 (17/12) assuming that both epoxide and hydroxyl are tightly bonded with carbon atoms. As the ether group in graphene is equivalent to vacancy, its phonon scattering rate can be determined by ΔM/M = 3.58 In principle, the higher the scattering rate, the smaller the thermal conductivity. With Klemens model, it can be predicted that the highest to lowest κ of graphene at the same defect concentration are κ(isotope) > κ(G−epoxide) > κ(G−OH) > κ(G−Ether). As shown in Figure 3a, MD simulation predicts that κ(isotope) > κ(G−OH) > κ(G−epoxide) > κ(G−Ether). As we have discussed earlier, MD simulation shows that κ(G− OH) is much larger than κ(G−Epoxide), whereas Klemens model suggests there is little difference between them. Through the comparison between Klemens model and MD simulation, we found that, in addition to the defect concentration, the configuration of oxygen defects have significant influence on phonon transport. To gain more insight on this effect, we have done phonon localization analysis to study the impact of oxygen defect configurations on phonon modes. In Figure 5, we have computed phonon participation ratios Pλ for pristine graphene, G−OH, G−Epoxide, and G−Ether with O/C ratio of 0.01. For different types of oxygen functional groups, we found that, when compared with Pλ of pristine

more effective than isotope defect in terms of scattering long wavelength phonons, and therefore, boundary scattering for phonons in GO become a secondary effect. On the basis of MD simulation results in Figure 3, it is clear that oxygen defect concentration plays a dominant role for thermal transport in GO. In addition, the specific configurations of oxygen functional groups in graphene (e.g., top site hydroxyl and bridge site epoxide) do have large influence over κ. In order to better understand the role of concentration and configuration of oxygen functional groups, we have done analysis on phonon density of states (DOS) and phonon localization effect. In MD simulation, we notice that κ reduction is very significant even for O:C ratio as low as 0.01. In order to understand the impact of oxygen functional group on phonon modes, we calculate phonon DOS by performing a fast Fourier transform of the velocity autocorrelation function.57 We compared phonon DOS between pristine graphene and GO with O/C ratio of 0.01. In Figure 4, we found G−OH has

Figure 4. Phonon DOS for pristine graphene (G−P), G−OH, G− Epoxide, and G−Ether.

Figure 5. Phonons participation ratio of the in-plane (i.e., LA, TA, LO, TO) and out-of-plane (i.e., ZA, ZO) vibration modes for pristine graphene (black), G−OH (blue), G−Epoxide (red), and G−Ether (green). The O/C ratio for GO is 0.01. The vertical axis value in each panel ranges from 0 to 1. 1440

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

■

graphene, ether group defects do not cause phonon localization for the in-plane and out-of-plane long wavelength acoustic phonons. However, epoxide and hydroxyl defects have demonstrated localization effect for the long wavelength phonons. In Figure 5, it is apparent that G−Epoxide has much stronger localization effect than G−OH. This result suggests that diffusive scattering is the primary scattering mechanism for G−Ether, but long wavelength phonon propagation in G−Epoxide is more spatially confined than that in G−OH. As a result, for the same O/C ratio in Figure 3a, MD simulation shows that κ(G−Epoxide) is always much lower than κ(G−OH), whereas Klemens model cannot distinguish their difference because of the lack of proper description on phonon localization effect in its model. The strong localization effect for G−Epoxide can arise from the oxygen binding configuration. In G−Epoxide, each oxygen atom is bonded with two carbon atoms at the bridge site. In G−OH, one OH radical is bonded with only one carbon atom at the top site. Therefore, epoxide functional group may apply stronger restriction on the vibration of carbon atoms than hydroxyl, and phonon localization for both in-plane and out-ofplane long wavelength phonon is apparent for G−Epoxide.

REFERENCES

(1) Balandin, A. A.; Ghosh, S.; Bao, W. Z.; Calizo, I.; Teweldebrhan, D.; Miao, F.; Lau, C. N. Superior Thermal Conductivity of SingleLayer Graphene. Nano Lett. 2008, 8, 902−907. (2) Balandin, A. A. Thermal Properties of Graphene and Nanostructured Carbon Materials. Nat. Mater. 2011, 10, 569−581. (3) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666−669. (4) Lee, C.; Wei, X.; Kysar, J. W.; Hone, J. Measurement of the Elastic Properties and Intrinsic Strength of Monolayer Graphene. Science 2008, 321, 385−388. (5) Yang, L.; Park, C. H.; Son, Y. W.; Cohen, M. L.; Louie, S. G. Quasiparticle Energies and Band Gaps in Graphene Nanoribbons. Phys. Rev. Lett. 2007, 99, 186801. (6) Haberer, D.; Vyalikh, D. V.; Taioli, S.; Dora, B.; Farjam, M.; Fink, J.; Marchenko, D.; Pichler, T.; Ziegler, K.; Simonucci, S.; et al. Tunable Band Gap in Hydrogenated Quasi-Free-Standing Graphene. Nano Lett. 2010, 10, 3360−3366. (7) Yan, J.-A.; Xian, L.; Chou, M. Y. Structural and Electronic Properties of Oxidized Graphene. Phys. Rev. Lett. 2009, 103, 086802. (8) Mathkar, A.; Tozier, D.; Cox, P.; Ong, P.; Galande, C.; Balakrishnan, K.; Reddy, A. L. M.; Ajayan, P. M. Controlled, Stepwise Reduction and Band Gap Manipulation of Graphene Oxide. J. Phys. Chem. Lett. 2012, 3, 986−991. (9) Hossain, M. Z.; Johns, J. E.; Bevan, K. H.; Karmel, H. J.; Liang, Y. T.; Yoshimoto, S.; Mukai, K.; Koitaya, T.; Yoshinobu, J.; Kawai, M.; et al. Chemically Homogeneous and Thermally Reversible Oxidation of Epitaxial Graphene. Nat. Chem. 2012, 4, 305−309. (10) Gao, W.; Alemany, L. B.; Ci, L.; Ajayan, P. M. New Insights into the Structure and Reduction of Graphite Oxide. Nat. Chem. 2009, 1, 403−408. (11) Park, S.; Ruoff, R. S. Chemical Methods for the Production of Graphenes. Nat. Nanotechnol. 2009, 4, 217−224. (12) Tian, H.; Xie, D.; Yang, Y.; Ren, T. L.; Zhang, G.; Wang, Y. F.; Zhou, C. J.; Peng, P. G.; Wang, L. G.; Liu, L. T. A Novel Solid-State Thermal Rectifier Based On Reduced Graphene Oxide. Sci. Rep. 2012, 2, 523. (13) Chang, C. W.; Okawa, D.; Majumdar, A.; Zettl, A. Solid-State Thermal Rectifier. Science 2006, 314, 1121−1124. (14) Li, B. W.; Wang, L.; Casati, G. Thermal Diode: Rectification of Heat Flux. Phys. Rev. Lett. 2004, 93, 184301. (15) Hu, J. N.; Ruan, X. L.; Chen, Y. P. Thermal Conductivity and Thermal Rectification in Graphene Nanoribbons: A Molecular Dynamics Study. Nano Lett. 2009, 9, 2730−2735. (16) Xiao, N.; Dong, X.; Song, L.; Liu, D.; Tay, Y.; Wu, S.; Li, L.-J.; Zhao, Y.; Yu, T.; Zhang, H.; et al. Enhanced Thermopower of Graphene Films with Oxygen Plasma Treatment. ACS Nano 2011, 5, 2749−2755. (17) Yao, Q.; Chen, L.; Zhang, W.; Liufu, S.; Chen, X. Enhanced Thermoelectric Performance of Single-Walled Carbon Nanotubes/ Polyaniline Hybrid Nanocomposites. ACS Nano 2010, 4, 2445−2451. (18) Turney, J. E.; McGaughey, A. J. H.; Amon, C. H. Assessing the Applicability of Quantum Corrections to Classical Thermal Conductivity Predictions. Phys. Rev. B 2009, 79, 224305. (19) Kong, B. D.; Paul, S.; Nardelli, M. B.; Kim, K. W. FirstPrinciples Analysis of Lattice Thermal Conductivity in Monolayer and Bilayer Graphene. Phys. Rev. B 2009, 80, 033406. (20) Lindsay, L.; Broido, D. A.; Mingo, N. Flexural Phonons and Thermal Transport in Multilayer Graphene and Graphite. Phys. Rev. B 2011, 83, 235428. (21) Nika, D. L.; Balandin, A. A. Two-Dimensional Phonon Transport in Graphene. J. Phys.: Condens. Matter 2012, 24, 233203. (22) Nika, D. L.; Askerov, A. S.; Balandin, A. A. Anomalous Size Dependence of the Thermal Conductivity of Graphene Ribbons. Nano Lett. 2012, 12, 3238−3244. (23) Ong, Z. Y.; Pop, E. Effect of Substrate Modes on Thermal Transport in Supported Graphene. Phys. Rev. B 2011, 84, 075471.

4. CONCLUSIONS We have performed MD simulations to investigate thermal transport in oxygen functionalized graphene. The interatomic interactions are determined by optimized REBO−CHO and REBO−carbon potentials. To validate the quality of the simulation, we have studied phonon scattering in graphene due to the effect of crystal boundary length and point defect, which are able to reflect the response of long and short wavelength phonon scattering, respectively. Our study shows that oxygen defect concentration is a major parameter for tuning the κ of GO. In addition, the configuration of oxygen defects (e.g., epoxide) can also significantly influence κ through phonon localization. Depending on the defect concentration and configuration considered in this study, κ of GO can be reduced by a factor ranging from 3 to 25. Through studying thermal transport property of GO, we hope that our simulation result can provide useful guidance for experiment to design oxidized graphene with precisely controlled thermal property for thermal management and thermoelectric application.

■

Article

AUTHOR INFORMATION

Corresponding Author

*(K.C.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by Nano Material Technology Development Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT, and Future Planning (2012M3A7B4049888). H.Z. and K.C. also acknowledge support from II−VI foundation. A.F.F. is a research fellow of the Brazilian agency ́ CNPq (Conselho Nacional de Desenvolvimento Cientifico e Tecnológico) and acknowledges grant #471116/2011-5 (CNPq) and grant #2012/10106-8, São Paulo Research Foundation (FAPESP). The authors would like to thank Dr. Susan B. Sinnott (University of Florida) for providing the parallel version of simulation code. 1441

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442

The Journal of Physical Chemistry C

Article

(24) Zhang, H.; Lee, G.; Cho, K. Thermal Transport in Graphene and Effects of Vacancy Defects. Phys. Rev. B 2011, 84, 115460. (25) Fthenakis, Z. G.; Tomanek, D. Computational Study of the Thermal Conductivity in Defective Carbon Nanostructures. Phys. Rev. B 2012, 86, 125418. (26) Chen, S. S.; Wu, Q. Z.; Mishra, C.; Kang, J. Y.; Zhang, H. J.; Cho, K. J.; Cai, W. W.; Balandin, A. A.; Ruoff, R. S. Thermal Conductivity of Isotopically Modified Graphene. Nat. Mater. 2012, 11, 203−207. (27) Bagri, A.; Kim, S.-P.; Ruoff, R. S.; Shenoy, V. B. Thermal Transport Across Twin Grain Boundaries in Polycrystalline Graphene from Nonequilibrium Molecular Dynamics Simulations. Nano Lett. 2011, 11, 3917−3921. (28) Serov, A. Y.; Ong, Z.-Y.; Pop, E. Effect of Grain Boundaries on Thermal Transport in Graphene. Appl. Phys. Lett. 2013, 102, 033104. (29) Chien, S. K.; Yang, Y. T.; Chen, C. K. Influence of Hydrogen Functionalization on Thermal Conductivity of Graphene: Nonequilibrium Molecular Dynamics Simulations. Appl. Phys. Lett. 2011, 98, 033107. (30) Kim, J. Y.; Lee, J.-H.; Grossman, J. C. Thermal Transport in Functionalized Graphene. ACS Nano 2012, 6, 9050−9057. (31) Liu, B.; Reddy, C. D.; Jiang, J. W.; Baimova, J. A.; Dmitriev, S. V.; Nazarov, A. A.; Zhou, K. Morphology and In-Plane Thermal Conductivity of Hybrid Graphene Sheets. Appl. Phys. Lett. 2012, 101, 211909. (32) Huang, W.; Pei, Q.-X.; Liu, Z.; Zhang, Y.-W. Thermal Conductivity of Fluorinated Graphene: A Non-Equilibrium Molecular Dynamics Study. Chem. Phys. Lett. 2012, 552, 97−101. (33) Fonseca, A. F.; Lee, G.; Borders, T. L.; Zhang, H. J.; Kemper, T. W.; Shan, T. R.; Sinnott, S. B.; Cho, K. Reparameterization of the REBO-CHO Potential for Graphene Oxide Molecular Dynamics Simulations. Phys. Rev. B 2011, 84, 075460. (34) Kemper, T. W.; Sinnott, S. B. Mechanisms of Ion-Beam Modification of Terthiophene Oligomers from Atomistic Simulations. J. Phys. Chem. C 2011, 115, 23936−23945. (35) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A., III. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (36) Schwamb, T.; Burg, B. R.; Schirmer, N. C.; Poulikakos, D. An Electrical Method for the Measurement of the Thermal and Electrical Conductivity of Reduced Graphene Oxide Nanostructures. Nanotechnology 2009, 20, 405704. (37) Mahanta, N. K.; Abramson, A. R. Thermal Conductivity of Graphene and Graphene Oxide Nanoplatelets. 13th IEEE InterSociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm), San Diego, CA, May 30−Jun 01, 2012; pp 1-6. (38) Lepri, S.; Livi, R.; Politi, A. On the Anomalous Thermal Conductivity of One-Dimensional Lattices. Europhys. Lett. 1998, 43, 271−276. (39) Wang, L.; Hu, B.; Li, B. Logarithmic Divergent Thermal Conductivity in Two-Dimensional Nonlinear Lattices. Phys. Rev. E 2012, 86, 040101. (40) Saito, K.; Dhar, A. Heat Conduction in A Three Dimensional Anharmonic Crystal. Phys. Rev. Lett. 2010, 104, 040601. (41) Nika, D. L.; Ghosh, S.; Pokatilov, E. P.; Balandin, A. A. Lattice Thermal Conductivity of Graphene Flakes: Comparison with Bulk Graphite. Appl. Phys. Lett. 2009, 94, 203103. (42) Lindsay, L.; Broido, D. A. Optimized Tersoff and Brenner Empirical Potential Parameters for Lattice Dynamics and Phonon Thermal Transport in Carbon Nanotubes and Graphene. Phys. Rev. B 2010, 81, 205441. (43) Ikeshoji, T.; Hafskjold, B. Nonequilibrium Molecular-Dynamics Calculation of Heat-Conduction in Liquid and Through Liquid-Gas Interface. Mol. Phys. 1994, 81, 251−261. (44) Bodapati, A.; Schelling, P. K.; Phillpot, S. R.; Keblinski, P. Vibrations and Thermal Transport in Nanocrystalline Silicon. Phys. Rev. B 2006, 74, 245207.

(45) Biswas, R.; Bouchard, A. M.; Kamitakahara, W. A.; Grest, G. S.; Soukoulis, C. M. Vibrational Localization in Amorphous-Silicon. Phys. Rev. Lett. 1988, 60, 2280−2283. (46) Mkhoyan, K. A.; Contryman, A. W.; Silcox, J.; Stewart, D. A.; Eda, G.; Mattevi, C.; Miller, S.; Chhowalla, M. Atomic and Electronic Structure of Graphene-Oxide. Nano Lett. 2009, 9, 1058−1063. (47) Acik, M.; Lee, G.; Mattevi, C.; Chhowalla, M.; Cho, K.; Chabal, Y. J. Unusual Infrared-Absorption Mechanism in Thermally Reduced Graphene Oxide. Nat. Mater. 2010, 9, 840−845. (48) Thomas, J. A.; Iutzi, R. M.; McGaughey, A. J. H. Thermal Conductivity and Phonon Transport in Empty and Water-Filled Carbon Nanotubes. Phys. Rev. B 2010, 81, 045413. (49) Pereira, L. F. C.; Donadio, D. Divergence of the Thermal Conductivity in Uniaxially Strained Graphene. Phys. Rev. B 2013, 87, 125424. (50) Singh, D.; Murthy, J. Y.; Fisher, T. S. On the Accuracy of Classical and Long Wavelength Approximations for Phonon Transport in Graphene. J. Appl. Phys. 2011, 110, 113510. (51) Bonini, N.; Garg, J.; Marzari, N. Acoustic Phonon Lifetimes and Thermal Transport in Free-Standing and Strained Graphene. Nano Lett. 2012, 12, 2673−2678. (52) Lindsay, L.; Broido, D. A.; Mingo, N. Flexural Phonons and Thermal Transport in Graphene. Phys. Rev. B 2010, 82, 115427. (53) Klemens, P. G. Theory of the A-Plane Thermal Conductivity of Graphite. J. Wide Bandgap Mater. 2000, 7, 332. (54) Zhang, H. J.; Lee, G.; Fonseca, A. F.; Borders, T. L.; Cho, K. Isotope Effect on the Thermal Conductivity of Graphene. J. Nanomater. 2010, 2010, 537657. (55) Che, J. W.; Cagin, T.; Goddard, W. A. Thermal Conductivity of Carbon Nanotubes. Nanotechnology 2000, 11, 65−69. (56) Xu, Z.; Bando, Y.; Liu, L.; Wang, W.; Bai, X.; Golberg, D. Electrical Conductivity, Chemistry, and Bonding Alternations under Graphene Oxide to Graphene Transition As Revealed by In Situ TEM. ACS Nano 2011, 5, 4401−4406. (57) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (58) Ratsifaritana, C. A.; Klemens, P. G. Scattering of Phonons by Vacancies. Int. J. Thermophys. 1987, 8, 737−750.

1442

dx.doi.org/10.1021/jp4096369 | J. Phys. Chem. C 2014, 118, 1436−1442