Electronics 2015, 4, 1109-1124; doi:10.3390/electronics4041109
OPEN ACCESS
electronics
ISSN 2079-9292
www.mdpi.com/journal/electronics
Article
Equilibrium Molecular Dynamics (MD) Simulation Study of
Thermal Conductivity of Graphene Nanoribbon: A
Comparative Study on MD Potentials
Asir Intisar Khan † , Ishtiaque Ahmed Navid † , Maliha Noshin, H. M. Ahsan Uddin,
Fahim Ferdous Hossain and Samia Subrina *
Department of Electrical and Electronic Engineering, Bangladesh University of Engineering and
Technology (BUET), Dhaka 1205, Bangladesh; E-Mails: asir.sakib@gmail.com (A.I.K.);
ishtiaque3nav@gmail.com (I.A.N.); riditapixy@gmail.com (M.N.);
anit.ahsan2011@gmail.com (H.M.A.U.); 010fahim@gmail.com (F.F.H.)
†
These authors contributed equally to this work.
* Author to whom correspondence should be addressed; E-Mail: samiasubrina@eee.buet.ac.bd;
Tel.: +880-19-3795-9083; +88-02-9668054; Fax: +88-02-9668054.
Academic Editor: Frank Schwierz
Received: 9 October 2015 / Accepted: 15 December 2015 / Published: 21 December 2015
Abstract: 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.
Keywords: thermal conductivity; graphene nanoribbon; equilibrium molecular dynamics;
Green-Kubo; optimized Tersoff potential; 2nd generation REBO potential
1. 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–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
Electronics 2015, 4
1111
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.
2. Theory and Simulation
2.1. 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
Electronics 2015, 4
1112
optimized Tersoff potential [23] because of their better accuracy in describing bond order characteristics
and anharmonicity of both sp2 and sp3 carbon systems including different carbon structures like graphene,
carbon nanotube and graphene nanoribbons.
2.2. 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 (Kx ) is related to HCACF by the following equation:
1
Kx =
V KB T 2
Zτ
hJx (t) · Jx (0)idt
(1)
0
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 Å), KB is the Boltzmann constant, T is
the system temperature and τ is the correlation time needed for reasonable decay of HCACF. Jx stands for
the x component of heat flux. The ensemble average of the heat flux autocorrelation function is presented
by the hJx (t) · Jx (0)i 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:
M
N
−m
X
∆t X
1
Kx =
Jx (m + n)Jx (n)
V KB T 2 m=1 N − m n=1
(2)
where ∆t is the MD simulation time step, M∆t corresponds to the correlation time τ of Equation (1) ,
Jx (m + n) and Jx (n) are the xth 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
J=
d X
ri (t)εi (t)
dt i
(3)
Here, ri (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 ith atom. For two body interactions, the total potential energy
content of a particle is defined in terms of the pair-wise interactions u2 (r) as
Epot =
1 X
u2 (rij )
2 ij,i6=j
(4)
and the total site energy is expressed by
1X
1
u2 (rij )
εi = m i v i 2 +
2
2 j
(5)
Electronics 2015, 4
1113
and subsequently the corresponding heat current can be calculated through the following expression:
J(t) =
X
εi v i +
i
1 X
(Fij · vi )rij ,
2 ij,i6=j
(6)
where vi is the velocity of atom i, rij = ri − rj and Fij is the interatomic force exerted by atom j on
atom i [18].
2.3. 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 ( Td ). As a result, thermal conductivity calculated by MD approach gives erroneous results
for temperatures below Td . Classical energy per carbon atom should be equal to phonon energy per carbon
atom at quantum temperature, Tq . Therefore,
Zωd
1
D (ω) [n(ω, Tq ) + ]hωdω = 3Kb TM D
2
(7)
0
where TM D is the temperature in MD simulation and n(ω, Tq ) is the occupation number of phonons
given by
1
n(ω, Tq ) =
e
h̄ω
k b Tq
(8)
−1
and the term 12 in Equation (7) represents the effect of zero point energy. Here, phonon density of states
3ω
(DOS), D(ω) = 2πN
, ωd = (4πN )1/2 v, N = number of carbon atoms in unit area of GNR and v is the
v2
effective phonon sound velocity satisfying
3
1
1
1
=
+
+
(9)
2
2
2
v
vLA
vT A
vZA 2
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
vLA = 19.5 km/s , vT A = 12.2 km/s and vZA = 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:
TM D
Tq3
=2 2
TD
TZ
D /Tq
x2
1
dx
+
TD ,
ex − 1
3
(10)
0
d
where TD = h̄ω
. TD is the Debye temperature which is 322 Kelvin. So, 2nd term of the equation
kb
becomes 107 K which represents that if MD simulation temperature is less than 107 K, there will be no
quantum corrected temperature.
Electronics 2015, 4
1114
900
Quantum Corrected Temperature with zero point energy
Molecular Dynamics Temperature
Quantum Corrected Temperature without zero point energy
800
Temperature (K)
700
600
500
400
300
200
100
0
0
100
200
300
400
500
600
700
800
900
Temperature (K)
Figure 1. Molecular Dynamics (MD) temperature versus Quantum Corrected temperature
for graphene nanoribbon (GNR).
The quantum corrected thermal conductivity KQC can be calculated by multiplying the uncorrected
thermal conductivity KM D by dTdTMqD ,
KQC =
dTM D
× KM D
dTq
(11)
From Figure 2, we can clearly see that quantum correction is dominant at low temperatures while at
higher temperatures it is almost negligible.
1
dTMD / dTq
0.8
0.6
0.4
0.2
0
0
200
400
600
Temperature (K)
800
Figure 2. Rate of change of MD temperature with respect to Quantum Corrected temperature
versus MD temperature for GNR.
Electronics 2015, 4
1115
2.4. 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 × 105 time steps with a time step of 0.5 fs followed by a switching to NVE ensemble for 2 × 107 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.
Figure 3. Schematic atomic structure of 10 nm×1 nm zigzag GNR.
3. Results and Discussion
3.1. 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
Electronics 2015, 4
1116
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.
3.2. 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.
Figure 4. Thermal Conductivity (uncorrected) of GNR (10 nm × 1 nm) as a function of
temperature using two different interatomic potentials.
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] .
Electronics 2015, 4
1117
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.
Thermal Conductivity (W/m−K)
8000
Optimized Tersoff (Uncorrected)
Optimized Tersoff (Corrected)
6000
4000
2000
0
100
200
300
400
Temperature(K)
500
(a) Optimized Tersoff
Thermal Conductivity (W/m−K)
3000
Second Generation REBO (Corrected)
Second Generation REBO (Uncorrected)
2500
2000
1500
1000
500
100
200
300
Temperature (K)
400
500
(b) 2nd Generation REBO
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.
Electronics 2015, 4
1118
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.
Quantum Corrected TC (W/m−K)
5000
Optimized Tersoff
Second Generation REBO
4500
4000
3500
3000
2500
2000
1500
1000
500
100
200
300
400
Quantum Corrected Temperature (K)
500
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.
3.3. 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.
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.
Electronics 2015, 4
1119
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.
Envelope of Normalized HCACF
1
8nm
10nm
12nm
13nm
14nm
0.8
0.6
0.4
0.2
0
0
2
4
6
Correlation Time (ps)
8
10
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.
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,
Electronics 2015, 4
1120
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 )vs where L is the in-plane size,
τ is the phonon relaxation time and vs 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].
3.4. 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].
Electronics 2015, 4
1121
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.
4. 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
Electronics 2015, 4
1122
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.
Author Contributions
Asir Intisar Khan; Ishtiaque Ahmed Navid; Maliha Noshin; H. M. Ahsan Uddin and
Fahim Ferdous Hossain jointly designed and performed the simulation under the supervision of Samia
Subrina All authors discussed and analyzed the data and interpretation, and contributed during the writing
of the manuscript. All authors have given approval to the final version of the manuscript.
Conflicts of Interest
The authors declare no conflict of interest.
References
1. Novoselov, K.; Geim, A.; Morozov, S.; Jian, D.; Katsnelson, M.; Grigorieva, I.; Dubonos, S.; Firsov, A.
Two-dimensional gas of massless Dirac fermions in graphene. Nature 2005, 438, 197–200.
2. Seol, J.H.; Jo, I.; Moore, A.L.; Lindsay, L.; Aitken, Z.H.; Pettes, M.T.; Li, X.; Yao, Z.;
Huang, R.; Broido, D.; et al. Two-Dimensional Phonon Transport in Supported Graphene. Science
2010, 328, 213–216.
3. 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.
4. Nakada, K.; Fujita, M. Edge state in graphene ribbons: Nanometer size effect and edge shape
dependence. Phys. Rev. B 1996, 54, 954–960.
5. Son, Y.W.; Cohen, M.L.; Louie, S.G. Energy Gaps in Graphene Nanoribbons. Phys. Rev. Lett.
2006, 97, 216803.
6. Han, M.Y.; Ozyilmaz, B.; Zhang, Y.; Kim, P. Energy Band-Gap Engineering of Graphene
Nanoribbons. Phys. Rev. Lett. 2007, 98, 206805.
7. Yang, D.; Ma, F.; Sun, Y.; Hu, T.; Xu, K. Influence of typical defects on thermal conductivity of
graphene nanoribbons: An equilibrium molecular dynamics simulation. App. Surf. Sci. 2012,
258, 9926–9931.
8. Balandin, A.A.; Ghosh, S.; Bao, W.; Calizo, I.; Teweldebrhan, D.; Miao, F.; Lau, C.N. Superior
Thermal Conductivity of Single-Layer Graphene. Nano Lett. 2008, 8, 902–907.
9. Singh, V.; Joung, D.; Zhai, L.; Das, S.; Khondaker, S.I.; Seal, S. Graphene based materials: Past,
present and future. Prog. Mater. Sci. 2011, 56, 1178–1271.
10. Hone, J.; Whitney, M.; Piskoti, C.; Zettl, A. Thermal conductivity of single-walled carbon nanotubes.
Phys. Rev. B 1999, 59, 2514–2516.
11. Cao, A. Molecular dynamics simulation study on heat transport in monolayer graphene sheet with
various geometries. App. Phys. Lett. 2012, 111, 083528.
Electronics 2015, 4
1123
12. Ghosh, S.; Calizo, I.; Teweldebrhan, D.; Pokatilov, E.P.; Nika, D.L.; Balandin, A.A.; Bao, W.;
Miao, F.; Lau, C.N. Extremely high thermal conductivity of graphene: Prospects for thermal
management applications in nanoelectronic circuits. App. Phys. Lett. 2008, 92, 151911.
13. Hu, J.; Ruan, X.; Chen, Y.P. Thermal Conductivity and Thermal Rectification in Graphene
Nanoribbons: A Molecular Dynamics Study. Nano Lett. 2009, 9, 2730–2735.
14. Chen, L.; Kumar, S. Thermal transport in graphene supported on copper. J. Appl. Phys. 2012,
112, 043502.
15. Zhang, Y.; Cheng, Y.; Pei, Q.; Wang, C.; Xiang, Y. Thermal conductivity of defective graphene.
Phys. Lett. A 2012, 376, 3668–3672.
16. Yu, C.; Zhang, G. Impacts of length and geometry deformation on thermal conductivity of graphene
nanoribbons. J. App. Phys. 2013, 113, doi:10.1063/1.4788813.
17. Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications;
Academic Press: San Diego, CA, USA, 2002.
18. Schelling, P.K.; Phillpot, S.R.; Keblinski, P. Comparison of atomic-level simulation methods for
computing thermal conductivity. Phys. Rev. B 2001, 65, 144306.
19. Khadem, M.H.; Wemhoff, A.P. Comparison of Green-Kubo and NEMD heat flux formulations for
thermal conductivity prediction using the Tersoff potential. Comp. Mater. Sci. 2013, 69, 428–434.
20. Mahdizadeh, S.J.; Goharshadi, E.K. Thermal conductivity and heat transport properties of graphene
nanoribbons. Springer 2014, 16, 1–12.
21. Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon.
Phys. Rev. Lett. 1988, 61, 2879–2882.
22. Evans, W.J.; Hu, L.; Keblinski, P. Thermal conductivity of graphene ribbons from equilibrium
molecular dynamics: Effect of ribbon width, edge roughness, and hydrogen termination.
Appl. Phys. Lett. 2010, 96, 203112.
23. 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.
24. Brenner, D.W. Empirical potential for hydrocarbons for use in simulating the chemical vapor
deposition of diamond films. Phys. Rev. B 1990, 42, 9458–9471.
25. Stuart, S.J.; Tutein, A.B.; Harrison, J.A. A reactive potential for hydrocarbons with intermolecular
interactions. J. Chem. Phys. 2000, 112, 6472–6486.
26. Zhang, H.; Lee, G.; Fonseca, A.F.; Borders, T.L.; Cho, K. Isotope Effect on the Thermal Conductivity
of Graphene. J. Nanomater 2010, 2010, 537657.
27. Brenner, D.W.; Shenderova, O.A.; Harrison, J.A.; Stuart, S.J.; Ni, B.; Sinnott, S.B.
A second-generation reactive empirical bond order (REBO) potential energy expression for
hydrocarbons. J. Phys. Condens. Matter 2002, 14, 783–802.
28. Lee, Y.H.; Biswas, R.; Soukoulis, C.M.; Wang, C.Z.; Chan, C.T.; Ho, K.M. Molecular-Dynamics
Simulation of Thermal Conductivity in Amorphous Silicon. Phys. Rev. B 1991, 43, 6573–6580.
29. Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys.
1995, 117, 1–19.
Electronics 2015, 4
1124
30. 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.
31. Yu, C.; Shi, L.; Yao, Z.; Li, D.; Majumdar, A. Thermal Conductance and Thermopower of an
Individual Single-Wall Carbon Nanotube. Nano Lett. 2005, 5, 1842–1846.
32. ZIMAN, J.M. Electrons and Phonons: The Theory of Transport Phenomena in Solids; Oxford
University Press: Amen House, London, UK, 1960.
33. Pereira, L.F.C.; Donadio, D. Divergence of the thermal conductivity in uniaxially strained graphene.
Phys. Rev. B 2013, 87, 125424.
34. Ecsedy, D.J.; Klemens, P.G. Thermal resistivity of dielectric crystals due to four-phonon processes
and optical modes. Phys. Rev. B 1977, 15, 5957–5962.
35. Steigmeier, E.F.; Kudman, I. Acoustical-Optical Phonon Scattering in Ge, Si, and III-V Compounds.
J. Appl. Phys. 1966, 141, 767–774.
36. Pichanusakorn, P.; Bandaru, P. Nanostructured thermoelectrics. Mater. Sci. Eng. R 2010, 67, 19–63.
37. Sonvane, Y.; Gupta, S.K.; Raval, P.; Lukacevic, I.; Thakor, P. Length, Width and Roughness
Dependent Thermal Conductivity of Graphene Nanoribbons. Chem. Phys. Lett. 2015, 634,
doi:10.1016/j.cplett.2015.05.036.
38. Nika, D.; Ghosh, S.; Pokatilov, E.; Balandin, A. Lattice Thermal conductivity of graphene flakes:
Comparison with bulk graphite. App. Phys. Lett. 2009, 94, 203103.
39. Nissimagoudar, A.S.; Sankeshwar, N.S. Significant reduction of lattice thermal conductivity due to
phonon confinement in graphene nanoribbons. Phys. Rev. B 2014, 89, 235422.
40. Klemens, P.G. Theory of the A-Plane Thermal Conductivity of Graphite. J. Wide Bandgap Mater.
2000, 7, 332–339.
41. Klemens, P.G.; Pedraza, D.F. Thermal Conductivity of Graphite in the Basal Plane. Carbon 1994,
32, 735–741.
42. 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.
43. Cao, H.Y.; Guo, Z.X.; Xiang, H.; Gong, X.G. Layer and size dependence of thermal conductivity in
multilayer graphene nanoribbons. Phys. Lett. A 2012, 376, 525–528.
© 2015 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article
distributed under the terms and conditions of the Creative Commons Attribution license
(http://creativecommons.org/licenses/by/4.0/).