Equilibrium Molecular Dynamics ( MD ) Simulation Study of Thermal Conductivity of Graphene Nanoribbon : A Comparative Study on MD Potentials

The thermal conductivity of graphene nanoribbons (GNRs) has been investigated using equilibrium molecular dynamics (EMD) simulation based on Green-Kubo (GK) method to compare two interatomic potentials namely optimized Tersoff and 2nd generation Reactive Empirical Bond Order (REBO). Our comparative study includes the estimation of thermal conductivity as a function of temperature, length and width of GNR for both the potentials. The thermal conductivity of graphene nanoribbon decreases with the increase of temperature. Quantum correction has been introduced for thermal conductivity as a function of temperature to include quantum effect below Debye temperature. Our results show that for temperatures up to Debye temperature, thermal conductivity increases, attains its peak and then falls off monotonically. Thermal conductivity is found to decrease with the increasing length for optimized Tersoff potential. However, thermal conductivity has been reported to increase with length using 2nd generation REBO potential for the GNRs of same size. Thermal conductivity, for the specified range of width, demonstrates an increasing trend with the increase of width for both the concerned potentials. In comparison with 2nd generation REBO potential, optimized Tersoff potential demonstrates a better modeling of thermal Electronics 2015, 4 1110 conductivity as well as provides a more appropriate description of phonon thermal transport in graphene nanoribbon. Such comparative study would provide a good insight for the optimization of the thermal conductivity of graphene nanoribbons under diverse conditions.

conductivity as well as provides a more appropriate description of phonon thermal transport in graphene nanoribbon.Such comparative study would provide a good insight for the optimization of the thermal conductivity of graphene nanoribbons under diverse conditions.

Introduction
In recent years, high density integration and size minimization have taken a tremendous turn in device technology.As a result, further development of the silicon-based micro-electronic device urges for the search of a new type of high thermal conductivity material.In this context, graphene has drawn a lot of attention as it has a number of unique properties which make it interesting for both fundamental studies and future applications.Graphene nanoribbons (GNRs), the narrow strips of graphene with few nanometers width, are particularly considered as a significant element for future nano-electronics.GNRs exhibit several intriguing electronic [1], thermal [2] and mechanical [3] properties dominated by their geometry i.e., width or edge structure [4][5][6].The ultra-thin characteristics of graphene and GNR are extremely beneficial for the aforementioned high-density integration.For a device to perform reliably, it is very important to manage the device heat efficiently.Therefore, it is necessary to investigate the thermal transport properties of graphene as well as graphene nanoribbons deeply [7].Phonon vibration is the prime thermal energy transport mechanism for graphene or GNRs [2,8,9].The contribution of phonons to the thermal conductivity is about 50-100 times greater than that of electrons [10].Particularly in recent experiments, it has been found using micro-Raman spectroscopy that single layer graphene sheets show extremely large values of thermal conductivity [11] .With the use of micro-Raman spectroscopy, the room temperature thermal conductivity for a single layer graphene suspended across a trench was found in the range of 3080-5150 W/m-K by Balandin [12].Furthermore, significant thermal rectification for triangular-shaped GNR has been reported [13].The thermal conductivity of graphene nanoribbons found in molecular dynamics simulation is of the same order of magnitude (2000 W/m-K) [13] as that of the experimental value (5300 W/m-K) [8], but the magnitude changes considerably with graphene dimension.For isolated 6 nm × 6 nm single layer graphene sheet, Chen et al. [14] reported 1780 W/m-K thermal conductivity while a value of 500 W/m-K for a graphene flake of 200 nm × 2.1 nm has been found in [15].Furthermore, for 27 nm × 15 nm graphene nanoribbon, the study of [7] obtained a thermal conductivity value of 5200 W/m-k whereas Yu et al. [16] observed 550 W/m-K thermal conductivity for a 600 nm × 10.4 nm graphene nanoribbon.Because of the extraordinary high value of thermal conductivity of graphene or GNR, deeper investigation on their thermal transport is necessary for the enhancement of energy-efficient nano-electronics.For the purpose of analyzing the thermal conduction property of graphene or GNR, molecular dynamics (MD) simulation can be employed.In our study, we have performed equilibrium molecular dynamics (EMD) simulation to compute thermal conductivity of GNR.EMD is based on Green-Kubo (GK) method derived from linear response theory [17].The thermal conductivities extracted from non-equilibrium MD (NEMD) method are generally one order smaller in magnitude than the experimentally obtained results following the similar trend [7].This is due to the fact that NEMD approach might demonstrate non-linear effects because of applied temperature gradient.This non linear effect may be subjected to the strong scattering caused by the heat source or heat sink used in the NEMD approach [18].Furthermore, size effects are more prominent in NEMD simulation in comparison with EMD simulation [19].EMD gives provision of computing the thermal conductivity along all directions in one simulation while NEMD can calculate thermal conductivity in one direction only as it uses thermal gradient in a particular direction.As a result, EMD is highly applicable for the geometries involving periodic boundary conditions.On the other hand, EMD is computationally more expensive and simulation results are more susceptible to parameters.However, the EMD with GK method is advantageous over NEMD because of inclusion of the entire thermal conductivity tensor from one simulation along with some additional parameters like heat current autocorrelation function (HCACF) [20].This is why, in order to systematically observe the thermal conductivity of GNRs, EMD simulation has been carried out in this paper.In a recent study, Mahdizadeh et al. [20] showed the variation of thermal conductivity of graphene nanoribbon (10 nm × 2 nm) with temperature and discussed the size dependence of thermal conductivity only with varying lengths (2 nm to 30 nm) using optimized Tersoff potential.However, the reliability of MD simulations highly depends on the use of appropriate interatomic energies and forces that can be advantageously interpreted by interatomic classical potentials.So, in this study, we have taken into account two of the MD potential fields: optimized Tersoff potential and 2nd generation REBO potential with a view to interpreting the heat transport mechanism in the graphene nanoribbon and thereby discussing the comparative reliability of the two potentials to provide the more appropriate description of the thermal conductivity of graphene nanoribbon.Using these two potential fields, we have performed a comparative analysis for the size dependence, i.e. varying lengths (8 nm to 14 nm) and widths (1 nm to 2.5 nm ) and temperature dependence of thermal conductivity of graphene nanoribbons.

Interatomic MD Potentials
Theoretical and simulation based approaches of phonon thermal transport in graphene, nanotubes and nanoribbons employing MD simulation require accurate representation of interactions between atoms.Interatomic potential, one of the most dictating fundamentals that influence the accuracy of thermal conductivity [11] can conveniently represent these interactions in a realizable form.Literature reported that the original Tersoff [21] potential tends to overestimate the thermal conductivity of graphene structures [22,23] while original REBO potential [24], adaptive intermolecular REBO or AIREBO [25] are found to underestimate the experimental data on the same [23,26].2nd generation REBO potential produces a more reliable function for simultaneously predicting bond lengths, energies, and force constants than the original version of the REBO potential [27].On the other hand, optimized Tersoff and Brenner potential [23] parameters predict acoustic phonon (majority heat carriers in a material) velocity values to be in better agreement with the experimental data.In our work, we have done comparative computations using 2nd generation reactive empirical bond order potential (REBO) [27] and optimized Tersoff potential [23] because of their better accuracy in describing bond order characteristics and anharmonicity of both sp 2 and sp 3 carbon systems including different carbon structures like graphene, carbon nanotube and graphene nanoribbons.

Equilibrium Molecular Dynamics Simulation: Green-Kubo Method
The computation of dynamic properties such as thermal conductivity in EMD is based on the fluctuation dissipation and linear response theorem [17].In equilibrium, the heat flow in a system of particles fluctuates around zero and this fact is directly used in this method.The calculation of the heat flux vectors and their correlations is carried out throughout the simulation [11,19].In this method, thermal conductivity (K x ) is related to HCACF by the following equation: where V stands for the system volume defined as the product of ribbon planar area and the nominal graphite interplanar separation, i.e., van der Waals thickness (3.4 Å), K B is the Boltzmann constant, T is the system temperature and τ is the correlation time needed for reasonable decay of HCACF.J x stands for the x component of heat flux.The ensemble average of the heat flux autocorrelation function is presented by the J x (t) • J x (0) term.To use the above equation in EMD simulation, the integration of the equation is represented in the discrete form as the following equation including the time averaging: where ∆t is the MD simulation time step, M∆t corresponds to the correlation time τ of Equation (1) , J x (m + n) and J x (n) are the x th components of the heat current in x direction at MD time-steps (m + n) and n respectively.N represents the total number of simulation steps while the number of time steps required for the heat flux correlation vector is denoted by M. Therefore, for obtaining better average from statistical point of view, M should be less than N.The heat current J used in Equation ( 2) for pair potential is defined as the time derivative of the sum of the moments of the site energies as Here, r i (t) is the time-dependent position of atom i and ε i (t) is the site energy which is the summation of kinetic energy and potential energy of the i th atom.For two body interactions, the total potential energy content of a particle is defined in terms of the pair-wise interactions u 2 (r) as and the total site energy is expressed by and subsequently the corresponding heat current can be calculated through the following expression: where v i is the velocity of atom i, r ij = r i − r j and F ij is the interatomic force exerted by atom j on atom i [18].

Quantum Correction
Quantum correction to the classical MD calculations of temperature and thermal conductivity is necessary because of the fact that quantum effects in classical MD approach are neglected below Debye temperature ( T d ).As a result, thermal conductivity calculated by MD approach gives erroneous results for temperatures below T d .Classical energy per carbon atom should be equal to phonon energy per carbon atom at quantum temperature, T q .Therefore, where T M D is the temperature in MD simulation and n(ω, T q ) is the occupation number of phonons given by and the term 1 2 in Equation ( 7) represents the effect of zero point energy.Here, phonon density of states (DOS), D(ω) = 3ω 2πN v 2 , ω d = (4πN ) 1/2 v, N = number of carbon atoms in unit area of GNR and v is the effective phonon sound velocity satisfying where LA, TA, and ZA stand for in-plane longitudinal, in-plane transverse and out-of-plane acoustic phonons respectively.Phonons of three acoustic branches are considered with phonon sound velocities v LA = 19.5 km/s , v T A = 12.2 km/s and v ZA = 1.59 km/s [28].This process aims at mapping classically calculated results onto their quantum analogs at the same energy level.Quantum corrected temperatures can be obtained from the following equation: where T D is the Debye temperature which is 322 Kelvin.So, 2nd term of the equation becomes 107 K which represents that if MD simulation temperature is less than 107 K, there will be no quantum corrected temperature.The quantum corrected thermal conductivity K QC can be calculated by multiplying the uncorrected thermal conductivity K M D by dT M D dTq , From Figure 2, we can clearly see that quantum correction is dominant at low temperatures while at higher temperatures it is almost negligible.

Simulation Details
In this work, equilibrium molecular dynamics simulations were carried out using LAMMPS [29] with proper periodic boundary conditions applied in the plane, i.e., in the length and width directions of graphene nanoribbon to avoid the effect of fixed walls.We have considered a zigzag graphene nanoribbon (GNR) of 10 nm × 1 nm (length × width) at room temperature which is shown in Figure 3. Two potential fields, namely optimized Tersoff potential and 2nd generation REBO potential were applied to study their effect on thermal transport.The system was thermalized using NoseHoover thermostat for 3 × 10 5 time steps with a time step of 0.5 fs followed by a switching to NVE ensemble for 2 × 10 7 time steps.During MD simulations the equations of motion were integrated with a velocity-Verlet integrator.Energy minimization was done using steepest descent algorithm due to its robustness.The heat current data were recorded every 5 steps in the micro canonical ensemble average state.Heat flux autocorrelation values were calculated by averaging five obtained HCACFs.Variable autocorrelation lengths were used for different sizes of GNRs to ensure the reasonable decay of normalized HCACF values.

Potential Dependence of Thermal Conductivity
Our results show that 2nd generation REBO potential underestimates the thermal conductivity of graphene nanoribbon by a considerable margin in comparison with that of optimized Tersoff potential.In fact, REBO potential measures lower thermal conductivity than even the original Tersoff parameters [23,30,31].Lattice thermal conductivity depends strongly on the phonon dispersion energies and near-zone centre velocities.With the 2nd generation REBO potential parameters, the velocities of the transverse acoustic (TA) branch and longitudinal acoustic (LA) branch are found to be very low while dispersion of the out of plane branch is also underestimated.In fact, 2nd generation REBO potential does not appropriately measure the zone centre velocities for all the acoustic modes.Strong anharmonicity of 2nd generation REBO potential resulting in high phonon-phonon scattering rates may also contribute to lower thermal conductivity [23].On the other hand, in our study, thermal conductivity of GNR using optimized Tersoff potential gives much better estimation which is in accordance with the expectation of Lindsay et al. [23].Using optimized Tersoff potential, our extracted value of thermal conductivity for 10 nm × 1 nm GNR at room temperature is 3207 W/m-K while using 2nd generation REBO potential measured thermal conductivity is 1650 W/m-K at 300K.Evans et al. [22] also reported a very high thermal conductivity (∼3000 W/m-K at room temperature) of 10 nm × 1 nm sized graphene nanoribbon using Tersoff potential which might be due to the limited number of phonon-phonon scattering in the small systems [32].Overall improved estimation of ZA, TA and LA phonon branches along with a better fit of near-zone-centre velocities by the optimized Tersoff parameters can be demonstrated to provide a better modeling of thermal conductivity of graphene nanoribbons in contrast to 2nd generation REBO potential.

Temperature Dependence of Thermal Conductivity
Figure 4 suggests that thermal conductivity of GNR decreases as temperature increases for both optimized Tersoff and 2nd generation REBO potential fields.This phenomena is a consequence of thermal conductivity reduction by phonon-phonon scattering in pure (without defect) lattice structures [33].For optimized Tersoff potential, the thermal conductivities obtained in our study at 300 K and 400 K are ∼3000 W/m-K and ∼1650 W/m-K respectively.As temperature increases, the number of high frequency phonon increases.Hence phonon-phonon anharmonic interaction (Umklapp scattering) increases that makes the decay of HCACF profile faster resulting in a decrease of thermal conductivity in EMD method.Our results follow the similar trend of magnitude with the variation of temperature as reported by Mahdizadeh et al. [20] (∼2500 W/m-K and ∼1700 W/m-K at 300 K and 400 K, respectively) where 10 nm × 2 nm sized GNR is considered.Our results show that the estimated thermal conductivity using optimized Tersoff potential is higher approximately by a factor of 2 compared to that with 2nd generation REBO potential near room temperature and this factor decreases with the increase of temperature.At lower temperature, thermal conductivity for optimized Tersoff potential varies with T −1 but it deviates at higher temperature [32].This might be due to an increased non-linear thermal resistivity considering three phonon-phonon interactions and comparatively weaker but appreciable four phonon-phonon process at an elevated temperature [34].In addition, anharmonic interactions of two acoustic modes with an optical phonon mode might lead to the observed variation of thermal conductivity at sufficiently high temperatures [35] .
Figure 5a,b shows the quantum corrected and uncorrected thermal conductivity of GNR as a function of temperature for optimized Tersoff potential and 2nd generation REBO potential, respectively.At low temperature (upto Debye temperature), quantum corrected thermal conductivity increases almost linearly as temperature increases and reaches a maximum and then drops again.The Umklapp process is supposed to be inactive at low temperature and therefore not available to provide thermal resistance [32] which has been considered through quantum correction.At room temperature and above, Umklapp scattering becomes highly significant [36] and thermally excited high energy phonons dominate the thermal conductivity.As a result, thermal conductivity decreases with the increase of temperature.Figure 6 shows the variation of quantum corrected thermal conductivity as a function of temperature for both optimized Tersoff potential and 2nd generation REBO potentials.The figure shows that the peak thermal conductivity for optimized Tersoff potential is higher than that of 2nd generation REBO potential.

Length Dependence of Thermal Conductivity:
In Figure 7, the thermal conductivity of graphene nanoribbon as a function of length for both optimized Tersoff and 2nd generation REBO potential fields has been depicted.The figure shows that the length dependence of thermal conductivity for these two fields is opposite in nature.The thermal conductivity decays with the increase in length for optimized Tersoff potential while it shows an increasing trend with the increase of length for 2nd generation REBO potential field.According to Figure 7, the thermal conductivity of GNR decreases with the increase of length in a monotonic manner for optimized Tersoff potential field which can be firstly interpreted from the HCACF profile plotted in Figure 8.With the increase in GNR length, the number of phonon increases resulting in the rise of phonon-phonon interactions.This causes a faster decay rate for HCACF, i.e., HCACF decays to zero more quickly.Secondly, the contribution of flexural phonons (ZA) to thermal transport of graphene has to be taken into consideration.Using optimized Tersoff potential, Lindsay et al. [23] showed the significance of ZA modes as the majority heat carriers in describing thermal transport in contrary to the established hypothesis of ignoring the out of plane phonon contribution.In this context, dispersion characteristics of ZA mode phonons in graphene are accountable for the convergence of thermal conductivity [33].ZA modes play the dual role of providing a significant source of heat carriers and at the same time preventing the divergence of thermal conductivity through heat carrier scattering.Their group velocity approaches zero near Γ point and hence, the scattering of other heat carriers becomes their prime role in heat transport [23,33].Thus, ZA modes in the low frequency region of power spectrum lower the in-plane thermal conductivity of GNR as demonstrated by the convergence trends using optimized Tersoff potential.However, converged sampling of the low frequency acoustic flexural modes is required to accurately calculate the thermal conductivity of GNR in molecular dynamics simulation.In this case, our findings disagree with some literature specified in the micrometer range of length [37,38].This might be due to poor sampling of phase space and poor ergodicity issues, i.e., lack of correspondence between the system's statistical average with the ensemble average.Here, the term ensemble average represents the average taken over all the states present in a system whereas the system statistical average implies the time average of a particular sequence of events rather than all the events or states present [17].However, Nika et al. [38] considered that the out of plane acoustic phonons (ZA) have insignificant contribution to thermal transport.However, recent studies by Nissimagoudar et al. [39] found that, ZA branch has significant influence in the thermal properties of graphene nanoribbons due to its quadratic dispersion.At the same time, Sonvane et al. [37] believes that flexural phonons dominate the thermal conductivity of graphene with length and width smaller than one micron which is in agreement with our study.
On the other hand, the thermal conductivity rises with the increase in length when the 2nd generation REBO potential is involved.2nd generation REBO potential underestimates the dispersion of ZA branches in contrast to the optimized Tersoff potential [23].Moreover, 2nd generation REBO potential includes a dihedral bonding term which poorly fits the phonon frequencies for the ZA branch.In this context, the length dependence of thermal conductivity using 2nd generation REBO potential is mainly governed by the phonon boundary scattering mechanism instead of the dispersion of ZA modes.The length of the GNR controls the phonon mean free path (MFP).With the increase in length, the acoustic phonons with longer wavelengths become available for heat transfer [37].The phonon transport in graphene, for phonon cutoff frequencies including zero frequency, is two dimensional (2D).In graphene, the long-wavelength phonons weakly scatter in the three-phonon Umklapp processes.This is calculated in the first order and yields the divergent nature of the thermal conductivity.Klemens [40,41], while working with the thermal conductivity of graphite, assumed negligible contribution of long wavelength i.e. low frequency phonons to the in plane thermal conductivity.To avert the problem of long-wavelength phonons in case of 2D lattices, Klemens [40] explained the phenomena using the size-dependent low-bound cut-off frequency ω min,s .ω min,s is related to finite in plane size by equation L = τ (ω min,s )v s where L is the in-plane size, τ is the phonon relaxation time and v s is the phonon group velocity [42].Thus he related his study with other literature showing a logarithmic divergence of thermal conductivity as a function of length for 2D lattices.In fact, with the small GNR length compared to phonon MFP (775nm for graphene), Umklapp scattering among phonons becomes negligible and the phonon collision at the edge becomes the significant scattering mechanism.Therefore, the shorter the GNR is, the stronger the edge scattering becomes which abates the thermal conductivity.Conversely, longer GNR results in the weakening of edge scattering leading to the increase of the thermal conductivity [7].

Width Dependence of Thermal Conductivity
Figure 9 shows the dependence of thermal conductivity on GNR width changing from 1nm to 2.5 nm at a fixed length of 10 nm.The thermal conductivity increases monotonically in the mentioned width range.The impacts of edge localized phonon effect or boundary scattering effect and the phonon's Umklapp effect, both having negative effects on the thermal conductivity, might be considered.The impact of boundary scattering is weakened with the increase in GNR width.But as the GNR width increases, the probability of Umklapp scattering rises as a result of increased number of phonons and the reduced energy separation between the phonon modes.The dominance between these two mechanisms dictates the nature of thermal conductivity variation.In case of narrower GNRs, the reduced edge-localized phonon effect is dominant with the increase of width leading to the increased thermal conductivity which is reflected in our study range.However, for a large enough width, more and more phonons will be activated with remarkably significant phonon's Umklapp effect and as a consequence, the thermal conductivity decreases with the increase of width of graphene nanoribbon [22,37,43].

Conclusions
In this study, using EMD simulation based on GK method, we have compared two interatomic potentials, optimized Tersoff and 2nd generation REBO in terms of thermal conductivity of graphene nanoribbons as a function of temperature, length and width for periodic boundary conditions.Our results show that thermal conductivity decreases with the increase of temperature for both optimized Tersoff and 2nd generation REBO potentials because of increasing high frequency phonon-phonon scattering and faster decaying of HCACF values.Furthermore, by introducing quantum correction to include the quantum effect, thermal conductivity increases with temperature up to Debye temperature, achieves its peak value and then falls off monotonically with temperature for both the potentials.Our study of thermal conductivity as a function of length demonstrates that, thermal conductivity decreases with increasing length for optimized Tersoff potential while an increasing trend is observed for 2nd generation REBO potential.Significant contribution of the dispersion of flexural mode phonons which are dominant in optimized Tersoff potential is responsible for convergence of thermal conductivity with increasing length.On the contrary, 2nd generation REBO potential highly underestimates the dispersion and phonon frequencies for the ZA branch as well as TA and LA velocities.As a result, the divergence nature of thermal conductivity of 2nd generation REBO potential is supposed to be dictated by phonon boundary scattering mechanism.It is found that thermal conductivity increases with the increasing width for both the potentials due to the dominant boundary scattering effect which decreases with increasing width.Finally, the estimated thermal conductivity using optimized Tersoff potential is higher approximately by a factor of 2 compared to that of 2nd generation REBO potential near room temperature.Optimized Tersoff potential parameters provide a better fit with various phonon branches, particularly with near zone centre velocities of these branches without altering structural data for graphene nanoribbon in comparison with the 2nd generation REBO potential.Therefore, optimized Terosoff potential has resulted in a better estimation of thermal conductivity and better description of phonon thermal transport of graphene nanoribbons compared to those of 2nd generation REBO potential.

Figure 2 .
Figure 2. Rate of change of MD temperature with respect to Quantum Corrected temperature versus MD temperature for GNR.

Figure 4 .
Figure 4. Thermal Conductivity (uncorrected) of GNR (10 nm × 1 nm) as a function of temperature using two different interatomic potentials.

Figure 5 .
Figure 5.Quantum corrected and uncorrected thermal conductivity of GNR (10nm×1nm) as a function of temperature using (a) optimized Tersoff and (b) 2nd generation Reactive Empirical Bond Order (REBO) potentials.

Figure 6 .
Figure 6.Quantum corrected thermal conductivity (TC) of GNR (10 nm length, 1 nm width) as a function of temperature using optimized Tersoff and 2nd generation REBO potentials.

Figure 7 .
Figure 7. Thermal conductivity variation of GNR with length using optimized Tersoff and 2nd generation REBO potentials at 300 K temperature.The width of GNR is 1 nm.

Figure 8 .
Figure 8. Envelopes of decaying heat current autocorrelation function (HCACF) profiles as a function of correlation length for varying lengths of GNRs (1nm width) using optimized Tersoff potential.

Figure 9 .
Figure 9. Thermal conductivity variation of GNR with width using optimized Tersoff and 2nd generation REBO potentials at 300 K temperature.The length of GNR is 10 nm.