0Impact of vacancies on the thermal conductivity of graphene nanoribbons: A molecular dynamics simulation study
Maliha Noshin, Asir Intisar Khan, Ishtiaque Ahmed Navid, H. M. Ahsan Uddin, and Samia Subrina
Citation: AIP Advances 7, 015112 (2017); doi: 10.1063/1.4974996 View online: http://dx.doi.org/10.1063/1.4974996 View Table of Contents: http://aip.scitation.org/toc/adv/7/1 Published by the American Institute of Physics
AVE YOU HEARD?
Employers hiring scientists and engineers trust
PHYSICS TODAY I JOBS
www.physicstoday.org/jobs
(■) CrossMark
VjPHf <-click for updates
Impact of vacancies on the thermal conductivity of graphene nanoribbons: A molecular dynamics simulation study
Maliha Noshin,1,a Asir Intisar Khan,1,a Ishtiaque Ahmed Navid,1 H. M. Ahsan Uddin,1,2 and Samia Subrina1,b
1 Department of Electrical and Electronic Engineering, Bangladesh University of Engineering and Technology, Dhaka 1205, Bangladesh
2Aix-Marseille Université, Jardin du Pharo, 58 Boulevard Charles Livon, 13284 Marseille, France
(Received 28 October 2016; accepted 16 January 2017; published online 25 January 2017)
Equilibrium molecular dynamics simulation using 2nd generation Reactive Bond Order interatomic potential has been performed to model the thermal transport of nanometer sized zigzag defected graphene nanoribbons (GNRs) containing several types of vacancies. We have investigated the thermal conductivity of defected GNRs as a function of vacancy concentration within a range of 0.5% to 5% and temperature ranging from 300K to 600K, along with a comparative analysis of those for pristine GNRs. We find that, a vacancy concentration of 0.5% leads to over 90% reduction in the thermal conductivity of GNRs. At low defect concentration, the decay rate is faster but ceases gradually at higher defect concentration. With the increasing temperature, thermal conductivity of defected GNRs decreases but shows less variation in comparison with that of pristine GNRs at higher temperatures. Such comprehensive study on several vacancy type defects in GNRs can provide further insight to tune up the thermal transport characteristics of low dimensional carbon nanostructures. This eventually would encourage the characterization of more stable thermal properties in thermal devices at an elevated temperature as well as the potential applicability of GNRs as thermoelectrics. © 2017 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). (http://dx.doi.org/10.1063/L4974996]
I. INTRODUCTION
The successful synthesis of graphene, a two-dimensional, single layer arrangement of carbon atoms in a honeycomb lattice configuration with sp2 bonds, triggered numerous technological advancements,1 owing to its remarkable electronic and thermal properties.2-4 As obtained from experimental results,5,6 it has been demonstrated that the single-layer graphene (SLG) exhibits strikingly high thermal conductivity in the range of 4840 W/m-k to 5300 W/m-K which is much greater than the thermal conductivity of all other known materials. Such a finding affirms the fact that graphene can very well be used as heat transmit units for electronic devices ranging from nano to micro scale taking into account the thermal management concerns. Experimentally fabricated graphene sheets are usually not perfect and the existence of several types of defects is inevitable. Defects in materials have profound impact on their physical and chemical properties.7 Graphene sheets can have defects in the form of vacancies, dislocations or impurities. Due to the highly durable and stable nature of graphene sheets, no noticeable concentration of point vacancy at temperatures below the melting point is observed. Actually, during the crystal growth process with the chemical vapor
bCorresponding Author: E-Mail: samiasubrina@eee.buet.ac.bd; ssubrina@gmail.com. Tel.: +880-19-3795-9083; +88-02-9668054; Fax: +88-02-9668054
l@ <D I
2158-3226/2017/7(1 )/015112/7 7,015112-1 © Author(s) 2017 t^KI^
deposition (CVD) technique, the defects in graphene are created.8 It is very significant to fully comprehend the intriguing characteristics of graphene and also to analyze and study the impact of several defects on graphene properties for its effective application purpose. In this context, simulation approaches for nano-scale material characterization technique are commonly used instead of experimental characterization technique because of its predicaments and complications in the set-up process.
Using NEMD, Mortazavi et at.1 showed the effects of randomly distributed defects on the thermal conductivity and tensile response of single layer graphene sheets. Yang et at.9 studied the thermal transport properties of hydrogenated and carbon isotope doped graphene by NEMD. They found that carbon isotope and hydrogen atoms will also cause severe reduction of thermal conductivity because of increased defect scattering. Classical molecular dynamics with the AIREBO potential is used to investigate the thermal conductivity of both zigzag and armchair graphene nanoribbons possessing different densities of Stone-Thrower-Wales (STW) defects by Ng et at.10 Their results indicate that the presence of the defects can decrease thermal conductivity by more than 50%. Using Tersoff type potential, the influence of defects on the thermal and electronic transport properties of graphene nanoribbons (GNRs) is explored by molecular dynamics and non-equilibrium Green's function methods in the study of Haskins etat.n NEMD simulation using AIREBO potential is used to investigate and compare the thermal conductivity of both zigzag and armchair graphene nanoribbons possessing various densities of defects, within a temperature range of 100K-600K in the study of Yeo et at.12 Their results indicate that the presence of both these kinds of defects can decrease the thermal conductivity by more than 80%.
Defects, especially in the form of vacancies can be of several types such as point vacancy, bi-vacancy, edge vacancy or edge roughness and a mixture of all these types of vacancies. A detailed study of the impact of all these types of vacancies on the thermal conductivity of graphene nanoribbons could be very effective to tune up their thermal transport characteristics for proper thermal applications. Hence, in this work, using EMD simulation with 2nd generation REBO13 potential, we have analyzed the change in thermal conductivity of GNR by varying the percentage of different vacancies and compared the obtained values with the thermal conductivity of pristine GNR. Subsequently, the study on the effect of varying temperature on the thermal conductivity of defected GNR has been carried out to have further insight on the properties of GNR with several types of defects.
II. SIMULATION DETAILS
In this study, equilibrium molecular dynamics simulations have been carried out using LAMMPS
(Large-scale Atomic/Molecular Massively Parallel Simulator)14 for the purpose of calculating the thermal conductivity of defected GNRs. To simulate the interatomic interactions in the MD simulations and to model the bonding interactions between carbon atoms in GNRs, we have employed 2nd generation REBO potential parameterized by Brenner et at.13 The computation of thermal conductivity in EMD, based on linear response theorem is related to the ensemble average of heat current
auto correlation (HCACF) function with the following Green-Kubo formulation
where Kx is the thermal conductivity in x direction, KB is the Boltzmann constant, T is the system temperature, V is the system volume defined as the area of GNR multiplied with the van der Waals thickness (3.3 Â) and t is the correlation time required for the reasonable decay of HCACF. The ensemble averaging term, for a two body potential, (Jx(t).Jx(0)) can be calculated from heat current equation,
where is the total site energy, v; is the velocity of atom i, r j = r - r;- (ri is the time-dependent position of atom i) and Fj is the force exerted by atom j on atom i.
FIG. 1. Schematic representation of a zigzag graphene nanoribbon including a) point vacancy b) bi-vacancy c) edge vacancy d) mixed vacancy (combination of point vacancy, bi-vacancy and edge vacancy).
We have considered 10 nm x 2 nm sized GNR structures with different types of vacancies as depicted in Figure 1. Point vacancy results from the removal of a single atom randomly from the lattice as shown in Figure 1(a) while a special type of point vacancy defined as edge vacancy in this study involves discarding atoms solely from the lattice boundary pictured in Figure 1(c). On the other hand, bi-vacancy originates from either the lumping of two consecutive point vacancies or simply from the random elimination of a pair of bonding atoms depicted in Figure 1(b). Figure 1(d) represents a mixture of all three types of vacancies i.e point vacancy, bi-vacancy and edge vacancy. In this study, defect concentrations are defined as the percentage ratio of the number of vacancies to the total number of atoms considered.
In our study, periodic boundary condition is applied in the plane i.e. in the length and width direction to avoid the effect of fixed walls or the boundary effects caused by the long mean free path. Alongside, to avoid interactions between the GNR in the simulation domain and its image, a vacuum region is incorporated. Equations of atomic motion were integrated using a velocity-Verlet integrator with a time step of 0.5 fs. The equilibration of the system was achieved using Nose-Hoover thermostat for 3 x 105 time steps followed by NVE ensemble for 105 time steps. For minimizing the system energy steepest decent algorithm has been used. The in-plane heat current data were recorded in every 5 steps and heat flux autocorrelation values were calculated by averaging the five obtained HCACFs. Using equation (1), the thermal conductivity has been calculated as the average of two orthogonal planar components.
In order to evaluate the phonon density of states (PDOS), Fix Phonon command15 of LAMMPS has been used to directly calculate the dynamical matrices from MD simulation based on fluctuation-dissipation theory using the motion of atoms in the equilibrated system. With the dynamical matrices obtained, PDOS was then calculated using an auxiliary post-processing code 'phana'. In this study, q (wave vector) points are generated uniformly and a tricubic16 interpolation method has been considered for the calculation of PDOS.
III. RESULTS AND DISCUSSIONS
Figure 2 represents the effect of defect on the thermal conductivity of GNRs at room temperature using 2nd generation REBO potential. As the figure suggests, for all types of vacancies, thermal conductivity decreases monotonically with the increase in the defect percentage. Bi-vacancy imposes a less dramatic reduction of thermal conductivity while point vacancy and edge vacancy i.e. point vacancy at the edge are the most destructive in terms of reduction in thermal conductivity which conforms to the findings of the previous studies.11,17,18 In case of point vacancy, the breaking of sp2 characteristic of the local lattice generates less stable two-coordinated atoms which causes higher degree of vibrations and scattering thereby reducing the thermal conductivity at a higher degree.11,19 Furthermore, the dominance of scattering due to the vacancies at the edge accompanied by edge roughness causes a higher percentage of decrease in thermal conductivity. Mixed defects including random point, edge and bi-vacancies are found to reduce the thermal conductivity most dominantly as can be seen from Figure 2. The degree of inelastic scattering around the vacancy centers and also at the distance from the vacancy centers may vary20 based on the type of vacancy. Hence, the impact
FIG. 2. Thermal conductivity as a function of defect percentage for 10 nm x 2 nm zigzag graphene nanoribbon with various types of vacancies at room temperature (300K).
of randomly mixed vacancy on thermal conductivity degradation is also dependent on the relative composition of the constituent vacancies. As a result, the overall impact of mixed vacancy on the thermal conductivity of GNR varies based on which type of vacancy is dominant in terms of percent composition in the mixed vacancy configuration.
Furthermore, for all types of vacancies, thermal conductivity profile shows a fast decrease at low defect concentrations (0.5% to 2.5%) and then the decay rate decreases gradually for higher defect concentrations upto 5%, as depicted in Figure 2. This phenomenon can be attributed to the reduction in phonon mean free path due to the scattering of phonons at the defect centers. At low defect concentration, both phonon-phonon scattering length and phonon-defect scattering length play a vital role to reduce the overall phonon mean free path. Moreover, when the defect concentration is low, defects serve as localized phonon scattering centers while at higher defect concentration, the scattering centers get spread throughout the whole nanoribbon leading to a delocalized scattering of phonons due to the interaction of different scattering centers.10,21
Figure 3 interprets the drastic reduction in thermal conductivity of GNRs with different types of vacancies in comparison with the same sized pristine GNR. For a defect free 10 nm x 2 nm GNR, the thermal conductivity was calculated to be -1850 W/m-K by Khan et al.22 using 2nd generation REBO potential. From Figure 3 it can be seen that, at a defect concentration of 0.5%, thermal conductivity experiences a reduction by nearly 92-95% for various defect structures while a reduction over 98%
FIG. 3. Percentage decrease of thermal conductivity in comparison with 10 nm x 2 nm pristine GNR for same sized GNRs with different types of vacancies with respect to variation in defect percentage.
at higher defect percentage (~5%) is observed. However, our considered small domain size for calculating thermal conductivity of graphene nanoribbons might predict somewhat higher thermal conductivity values with a possible reason of limited number of phonon-phonon scattering in the small systems in comparison with the large systems. This is also in line with the study of Mahdizadeh et al.23 and Evans et al.24 where similar domain size has been considered to simulate the thermal conductivity of pristine graphene nanoribbon. For a 10 nm x 2 nm pristine zigzag GNR, the study of Mahdizadeh et al.23 calculated a thermal conductivity value of around 2500 W/m-K whereas for 10 nm x 1 nm sized GNR, Evans et al.24 reported a high thermal conductivity value of ~3000 W/m-K using EMD.
Literature7,12,18,21 reported similar exponential decay of thermal conductivity in the presence of various types of defects. The decrease in thermal conductivity in the presence of vacancy can firstly be attributed to the localization of long wave length phonons around a vacancy in GNR.12 These long wavelength i.e low frequency phonons are the majority heat carriers in graphene and in the presence of vacancy the localization of these important heat carriers shrinks the thermal transport capability of GNRs. Moreover, the high frequency optical phonon modes also experience strong inelastic scattering at around the vacancy centers and also at a distance from the vacancy centers, thereby reducing the thermal conductivity of GNRs in the presence of different types of defects such as vacancies.20 Secondly, phonon density of states of the defected graphene with different vacancy concentrations can be taken into consideration as reported in the previous literature.7,12,18
Figure 4(a) shows the phonon density of states of a pristine GNR while Figure 4(b) presents the phonon density of states of the same sized GNR structure with 2% random mixed vacancies at room temperature. As can be seen from the figure, with the inclusion of vacancy, the PDOS peaks in the high frequency region get softened as opposed to those of pristine graphene. This phenomenon forces a reduction of the lifetime hence, mean free path of the corresponding modes. Moreover, some finite peaks appear in the low-frequency regions of the PDOS which may also cause a further reduction in the relaxation time and mean free path thereby reducing the overall thermal conductivity. The shift of the PDOS towards the low frequency region can be attributed to the increase in the number of the unsaturated C (Carbon) atoms with some dangling bonds in GNRs in the presence of different types of vacancies.
Figure 5(a) shows a decaying thermal conductivity trend with respect to the increasing temperature for defected GNRs with defect concentration of 0.5%. The decaying trend is similar for all three types of vacancies as well as for mixed vacancies. Apart from phonon scattering due to defects or vacancies, temperature dependance of thermal conductivity can be explained from the view point of phonon phonon anharmonic interaction also known as Umklapp scattering mechanism. As temperature increases, the number of high frequency phonons increases. Hence phonon-phonon anharmonic interaction (Umklapp scattering) increases.25,26 As a result, because of the increased scattering, thermal conductivity decreases for increasing temperatures. Secondly, as shown in
FIG. 4. Phonon density of states, PDOS at room temperature for 10 nm x 2 nm zigzag GNR (a) with pristine structure and (b) structure with random mixed vacancy (2% defect concentration).
FIG. 5. a) Thermal conductivity for 10 nm x 2 nm GNRs with different types of vacancies (defect concentration 0.5%) with respect to variation in temperature, b) Envelopes of decaying heat current autocorrelation function (HCACF) profiles as a function of correlation time for GNRs (10 nm x 2 nm) with point vacancies (0.5% defect concentration) for varying temperatures.
Figure 5(b), decaying HCACF profiles as a function of correlation time for increasing temperatures can be taken into consideration. As the temperature increases, HCACF decays to zero more quickly and thermal conductivity decreases consequently for defected GNRs.
But, with the increase in temperature, the variation in thermal conductivity for defected GNRs is less than that of similar sized pristine GNRs using the same potential as reported by Khan et al.22 With the increasing temperature, thermal conductivity shows relatively flat curves as opposed to that of pristine GNR case.
For a fixed defect concentration at a specific temperature, percentage decrease of thermal conductivity is defined here by
TCpristine,T' — TCdefected ,T'
TCpristine,T'
where TCpristine,T' is the thermal conductivity of pristine GNR at a certain temperature (T') and TCdefected,T > is the thermal conductivity of a defected GNR (0.5% vacancy concentration in this case) at that temperature.
751-■-1-1-
300 400 500 600 700 Temperature (K)
FIG. 6. Percent decrease of thermal conductivity of GNRs (10 nm x 2 nm) with various types of vacancies with respect to variation in temperature. At each temperature, percentage decrease in thermal conductivity for defected GNRs (0.5% defect concentration) has been evaluated with respect to the calculated value of thermal conductivity of pristine GNR at that temperature.
From Figure 6, we observe that, percentage decrease of thermal conductivity decreases with increasing temperatures. This is due to the fact that at higher temperatures, phonon-phonon scattering becomes more dominant than scattering of phonons by defects in the GNR structure. As a result, percentage decrease of thermal conductivity at higher temperatures is comparatively lower in comparison with that at lower temperatures.
IV. CONCLUSIONS
Using equilibrium molecular dynamics simulations, we have demonstrated the drastic impact of defects on the thermal conductivity of GNRs containing several types of random vacancies namely point vacancy, bi-vacancy, edge vacancy and mixed vacancy. Thermal conductivity experiences a reduction by nearly 92-95% for various defect structures at a defect concentration of 0.5%, while a reduction over 98% occurs at a higher defect percentage. Localization of long wavelength phonons around the vacancy centers and shifting of the PDOS peaks towards the low frequency region as well as the softening of phonon modes at high frequency regions can be attributed to the drastic reduction of thermal conductivity of GNRs in the presence of vacancy. Point vacancies and edge vacancies are found to be the most destructive ones in terms of reduction in thermal conductivity. For all types of vacancies, thermal conductivity decreases sharply at low defect concentration due to strong phonon-phonon and phonon-defect scattering. But at higher defect concentrations, because of the delocalized scattering of phonons due to the interaction of different scattering centers, the decay rate gradually ceases. Further, with the increase of temperature, thermal conductivity of defected GNR decreases. This can be explained by the increased phonon-phonon scattering by vacancies, dominant Umklapp scattering and faster decay of HCACF profile with the increasing temperature. But, thermal conductivity profile of defected GNRs experience less variation in comparison with that of pristine GNRs. In fact, at higher temperatures, phonon-phonon scattering dominates over phonon-defect scattering which can be attributed to the smaller decrease in percentage of thermal conductivity of defected GNRs with respect to pristine GNRs at elevated temperatures. Apart from providing a better understanding of the thermal transport characteristic of defected GNRs including several types of vacancies, our results could help in finding an efficient way to tune the thermal conductivity of GNRs for the future practical design of devices based on low dimensional carbon nanostructures.
1 A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007).
2Z. Yuanbo, T. Yan-Wen, H. L. Stormer, and P. Kim, Nature 438, 201 (2005).
3 S. Stankovich, D. A. Dikin, G. H. B. Dommett, K. M. Kohlhaas, E. J. Zimney, E. A. Stach, R. D. Piner, S. T. Nguyen, and R. S. Ruoff, Nature 442, 282 (2006).
4 A. K. Geim and P. Kim, Sci. Am. 298, 90 (2008).
5 A. A. Balandin, S. Ghosh, W. Bao, I. Calizo, D. Teweldebrhan, F. Miao, and C. N. Lau, Nano Lett. 8, 902 (2008).
6 S. Ghosh, W. Bao, D. L. Nika, S. Subrina, E. P. Pokatilov, C. N. Lau, and A. A. Balandin, Nat. Mater. 9, 555 (2010).
7 B. Mortazavi and S. Ahzi, Carbon N. Y. 63, 460 (2013).
8 D. Zhang, B. Hu, D. Guan, and Z. Luo, Catal. Commun. 76, 7 (2016).
9 D. Yang, F. Ma, Y. Sun, T. Hu, and K. Xu, Appl. Surf. Sci. 258, 9926 (2012).
10 T. Y. Ng, J. J. Yeo, and Z. S. Liu, Carbon N. Y. 50, 4887 (2012).
11 J. Haskins, A. Kinaci, C. Sevik, H. Sevincli, G. Cuniberti, and T. Cagin, ACS Nano 5, 3779 (2011).
12 J. J. Yeo, Z. Liu, and T. Y. Ng, Nanotechnology 23, 385702 (2012).
13 D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J. Phys. Condens. Matter 14, 783 (2002).
14S. Plimpton, J. Comput. Phys. 117, 1 (1995).
15 L. T. Kong, Comput. Phys. Commun. 182, 2201 (2011).
16 F. Lekien and J. Marsden, Int. J. Numer. Methods Eng. 63, 455 (2005).
17 M. H. Gass, U. Bangert, A. L. Bleloch, P. Wang, R. R. Nair, and A. K. Geim, Nat. Nanotechnol. 3, 676 (2008).
18 H. Zhang, G. Lee, and K. Cho, Phys. Rev. B - Condens. Matter Mater. Phys. 84, 1 (2011).
19 W. Zhao, Y. Wang, Z. Wu, W. Wang, K. Bi, Z. Liang, J. Yang, Y. Chen, Z. Xu, and Z. Ni, Sci. Rep. 5, 11962 (2015).
20 G. C. Loh, E. H. T. Teo, and B. K. Tay, Diam. Relat. Mater. 23, 88 (2012). 21F. Hao, D. Fang, and Z. Xu, Appl. Phys. Lett. 99, 41901 (2011).
22 A. Khan, I. Navid, M. Noshin, H. Uddin, F. Hossain, and S. Subrina, Electronics 4, 1109 (2015).
23 S. J. Mahdizadeh and E. K. Goharshadi, J. Nanoparticle Res. 16, 2553 (2014). 24W. J. Evans, L. Hu, and P. Keblinski, Appl. Phys. Lett. 96, 203112 (2010).
25 P. Pichanusakorn and P. Bandaru, Mater. Sci. Eng. R Reports 67, 19 (2010).
26 J. M. Ziman, Electrons and Phonons: The Theory of Transport Phenomena in Solids (Oxford University Press, Amen House, London, UK, 1960).
