phys. stat. sol. (b) 245, No. 4, 695–700 (2008) / DOI 10.1002/pssb.200743461 pss solidi physica status b www.pss-b.com basic solid state physics Stability analysis of graphene nanoribbons by molecular dynamics simulations N. Dugan and Ş. Erkoç* Department of Physics, Middle East Technical University, 06531 Ankara, Turkey Received 7 November 2007, revised 14 December 2007, accepted 14 December 2007 Published online 1 February 2008 PACS 31.15.xv, 61.48.De, 81.05.Uw * Corresponding author: e-mail erkoc@erkoc.physics.metu.edu.tr, Phone: +90 312 210 32 85, Fax: +90 312 210 50 99 In this work, stability of graphene nanoribbons are investi- gated using molecular dynamics. Simulations include heating armchair and zigzag-edged nanoribbons of widths varying be- tween one and nine hexagonal rings until the bonds between carbon atoms start to break. Breaking temperatures and bind- ing energies per atom for different widths are presented for both armchair and zigzag-edged cases. A nontrivial relation between stability and width is observed and discussed. Armchair-edged graphene nanoribbon at 1 K and 3100 K. © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim 1 Introduction Graphene, the two-dimensional coun- zigzag-edged GNRs are always metallic whereas armchair- terpart of three-dimensional graphite has recently drawn edged GNRs can be metallic or insulating depending on much attention due to its electronic properties. Some of the their widths [14–16]. As an additional effect of the width theoretical interest is focused on the use of relativistic on the electronic properties of GNRs, the band gap is Dirac equation for explaining electronic properties and on found to decrease with increasing width for all kinds of the observation of quantum Hall effect in this two- GNRs [17]. In a recent Monte Carlo study, structural dimensional material [1–4]. The value of minimum con- changes in relatively larger graphene sheets has been simu- ductivity [5] and the reason for the linear dependence of lated [18]. conductivity on the electron density are some of the open Stability is an important issue for a material which is questions related to graphene [6]. Practical interest, on the intended to be used as a building block of electrical circuits. other hand, lies in the possible applications in the field of Even though perfect two-dimensional crystals are proven carbon based electronics [7], where most of the studies are to be unstable, graphene is found to be stabilized by corru- concentrated on carbon nanotubes [8, 9]. Production of gations in the third dimension [19]. Understanding the ef- nanotubes however is not an easy task and the confinement fect of GNR width on the stability is therefore also central effects due to the finite size of tube thickness is not a prop- issue for the possible applications in the near future. In this erty restricted to the tube geometry. Planar graphene study, structural stability of armchair and zigzag-edged nanoribbons (GNRs) also display confinement properties GNRs of widths varying between one and nine hexagonal as mentioned by Walt de Heer from Georgia Institute of rings were investigated using molecular dynamics (MD). Technology [10]. Graphene has drawn attention after it Details of the method are given in the next section together was first prepared in 2005 by micromechanical cleavage with the potential energy function used in the calculations. method [11]. Effects of edges on electronic properties of GNRs have been investigated recently [3, 12, 13]. Due to 2 Method of calculation The Tersoff empirical po- the localized edge states that lie close to the Fermi level, tential energy function developed for carbon [20] is used in © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim  pss b solidi physica status 696 N. Dugan and Ş. Erkoç: Stability analysis of graphene nanoribbons the calculations. Total interaction energy of the system of Table 1 Final temperatures and binding energy per particle val- particles is taken to be sum of total two-body and total ues at 1 K, 300 K and final temperature of armchair edged gra- three-body contributions [21] phene ribbons. Width specifies number of hexagonal rings across the finite dimension. Energy values are given in electronvolt units. Φ = φ2 + φ3 . (1) width final T (K) – E/par – E/par – E/par Total two-body and three-body energies are expressed, re- – (1 K) – (300 K) – (final T) spectively, as 1 3300 – 5.48 – 5.25 – 5.18 N 2 3900 – 5.85 – 5.77 – 5.62 φ2 = A  U (1) ij , (2) 3 3400 – 6.03 – 5.83 – 5.56 i< j 4 3500 – 6.17 – 6.18 – 5.83 -1/ 2 n 5 3600 – 6.26 – 6.06 – 5.87 N È Ê N ˆ n˘ Í ˜ ˙˙ 6 3700 – 6.33 – 5.90 – 5.90 φ3 = - B  U (2) Í ij Í 1+ β nÁ Á ÂW ijk ˜˜ ˙ . (3) 7 3500 – 6.38 – 6.12 – 6.06 Í ÁÁ ˜¯ ˙ i< j Ë k πi , j ÎÍ ˙˚ 8 3300 – 6.43 – 6.17 – 5.94 9 3100 – 6.46 – 6.12 – 6.11 Here U ij and Wijk represent the two-body and three-body interactions, respectively: 3 Results and discussion MD calculations de- U ij(1) = f c (rij ) exp ( -λ1rij ) , scribed in the previous section were carried out for arm- chair and zigzag-edged thin GNRs with periodic boundary U ij(2) = f c (rij ) exp ( -λ2 rij ) , (4) conditions applied in one dimension. Width of the finite dimension was varied from one hexagonal ring to nine Wijk = f c ( rik ) g (θ ijk ) , hexagonal rings. Images are presented in Figs. 1 and 2 re- spectively for both armchair and zigzag edged GNRs. Two where different views at 1 K, at 300 K and at the final breaking c2 c2 temperature for widths of 1, 2, 5 and 9 hexagonal rings are g (θ ijk ) = 1 + 2 - 2 , (5) given in these figures. Final temperature values are dis- d d + (h - cos θ ijk ) 2 played in Tables 1 and 2 together with binding energies per atom. Slight corrugations are observed in the third dimen- Ï1 , r <R-D sion at 300 K for all configurations of GNRs studied. Am- ÔÔ 1 1 π r - R˘ f c (r ) = Ì - sin ÈÍ , R-D<r <R+D . plitudes of corrugations are observed to increase with tem- Ô 2 2 Î 2 D ˚˙ perature and the GNRs become highly distorted at the final temperature. In our study this final temperature at which ÓÔ0 , r >R+D the GNR distorts with broken bonds is regarded as a meas- (6) ure of stability. The parameters of the PEF for carbon are as follows Final temperature values in Tables 1 and 2 which [20]: A = 1393.6 eV, B = 346.74 eV, λ1 = 3.4879 Å, can also be seen in Figs. 3 and 4 do not indicate a simple λ2 = 2.2119 Å, β = 1.5724 ¥ 10 -7 , n = 0.72751 , c = 38049, relation between width and stability. Stability increases d = 4.3484 , h = -0.57058, R = 1.95 Å, and D = 0.15 Å. by a significant amount when width increases from one Verlet integration algorithm [22] is used in the canoni- hexagonal ring to two hexagonal rings for both types of cal ensemble (NVT) MD simulations [23]. The tempera- ture of the system is related to the time average of the ki- Table 2 Final temperatures and binding energy per particle val- netic energy to scale the velocities. This type of scaling ues at 1 K, 300 K and final temperature of zigzag edged graphene procedure is usually called as the Woodcock thermostat ribbons. Width specifies number of hexagonal rings across the fi- [24]. Temperature is kept constant by applying temperature nite dimension. Energy values are given in electronvolt units. scaling at every other time step which is chosen to be 10-16 seconds. Initial velocities of the particles are deter- width final T (K) – E/par – E/par – E/par mined according to the Maxwell–Boltzmann distribution. – (1 K) – (300 K) – (final T) Simulations start from a temperature of 1 K and 100 K 1 3600 – 5.90 – 5.91 – 5.10 temperature increment is applied after each 40000 MD 2 4800 – 6.16 – 6.12 – 5.69 steps. Temperature is increased until bonds between car- 3 4300 – 6.30 – 5.88 – 5.88 bon atoms start to break, causing the GNR to deform. Ini- 4 4100 – 6.38 – 6.36 – 6.00 tialization of bond breaking is defined based on the three- 5 4300 – 6.43 – 5.75 – 6.05 body term in the potential energy function. It is assumed 6 3200 – 6.47 – 6.49 – 6.11 that a bond breaking occurs if two neighbors for an atom 7 3500 – 6.50 – 5.83 – 6.15 can not be found in the cut off distance determined by the 8 3000 – 6.52 – 5.80 – 6.10 parameters R and D. 9 3200 – 6.54 – 5.93 – 6.02 © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.pss-b.com  Original Paper phys. stat. sol. (b) 245, No. 4 (2008) 697 Figure 1 Armchair edged graphene ribbons. Left column shows top view and right column shows tilted view. Ribbons of width 1, 2, 5 and 9 hexagonal rings are shown. 1 K, 300 K and final temperature images are given for each width. GNRs, being more significant for zigzag edged GNRs. nine rings. In a recent density functional theory (DFT) stu- However, increase in stability does not continue monotoni- dy on the zero temperature stability of GNRs up to 3 nm cally throughout the range of widths studied. For armchair- width [25], only an increase in the stability has been ob- edged GNRs (Fig. 3), omitting the instance of width two, served with increasing width. Despite the range of 3 nm stability increases almost linearly up to a width of six rings width is being far beyond the critical width of 0.861 nm at after which it decreases again linearly up to the width of which stability of armchair GNRs start to decrease in our www.pss-b.com © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim  pss b solidi physica status 698 N. Dugan and Ş. Erkoç: Stability analysis of graphene nanoribbons Figure 2 Zigzag edged graphene ribbons. Left column shows top view and right column shows tilted view. Ribbons of width 1, 2, 5 and 9 hexagonal rings are shown. 1 K, 300 K and final temperature images are given for each width. © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.pss-b.com  Original Paper phys. stat. sol. (b) 245, No. 4 (2008) 699 4000 Table 3 Distance ( Dh ) between the atoms which are most sepa- rated from average GNR plane at two sides for armchair-edged GNRs of different widths. 3800 width final T Dh (1 K) Dh (300 K) Dh (final T) final temp (K) 3600 1 3300 0.0004 2.0244 3.0874 2 3900 0.0013 2.4611 3.3935 3400 3 3400 0.0039 1.8748 3.2370 4 3500 0.0094 1.0121 3.3275 5 3600 0.0099 1.9879 4.5590 3200 6 3700 0.0201 3.3341 4.7173 7 3500 0.0142 1.9025 5.8760 8 3300 0.0093 2.2187 5.0753 3000 9 3100 0.0117 2.4507 5.8555 0 2 4 6 8 10 width (hexagonal rings) be more important compared to the symmetry argument Figure 3 Width versus final temperature graph for armchair described above and avoids appearance of even- edged GNRs. ness/oddness effects for armchair-edged GNRs. Another interesting result obtained is that GNRs of a MD study, there is no evidence of stability decrease for width of two hexagonal rings are the most stable for both wider GNRs in the mentioned DFT study. In the current types of edges. This can be explained by the rotational high temperature MD simulations dynamical effects such freedom about the infinite axial dimension of the ribbons. as rotational freedom are also being considered which may GNRs of a width of a single ring also have this rotational be the reason of stability decrease at higher widths. freedom but their stability is reduced due to their extremely For zigzag-edged GNRs (Fig. 4), the stability displays thin finite dimension. a more complicated behavior when width increases. An in- The increase followed by a decrease in stability with crease followed by a decrease in stability can also be pro- GNR width may thus be explained as follows: As the posed for this case but it is not as sharp as the armchair width increases, rotational freedom disappears gradually edged case. Evenness/oddness of the width also seems to and this situation causes a decrease in stability. However be important for the stability of zigzag edged GNRs. Omit- being wider is also a property that increases stability if ro- ting once again the instance of width two, it is observed tational freedom is not considered. Combination of these that zigzag-edged GNRs that have an odd number of hex- two effects result in the complicated stability width rela- agonal rings are more stable compared to their even- tion as seen in Figs. 3 and 4. Absence of rotational freedom numbered neighbors. Symmetry with respect to the central becomes more important for wider GNRs and therefore axis of finite dimension possessed by the GNRs with a their stability are lower compared to thinner ones. Beside width of odd-numbered rings may be responsible for the being non-monotoning, there is an increase in binding higher stability observed in this case. This symmetry ar- energy per atom values with GNR width, as seen in Ta- gument is valid also for the armchair-edged GNRs. How- bles 1 and 2. Therefore there is no clue about the lower sta- ever, the fact that their edges being highly curved seems to bility for larger widths in the binding energy per atom val- ues. 4800 We calculated height fluctuation for each GNR as the distance ( Dh ) between the atoms which are most separated 4500 from average GNR plane at two sides. Value of this dis- tance is approximately zero for 1 K and increases with 4200 temperature. These Dh values are given in Tables 3 and 4 final temp (K) for armchair and zigzag edged GNRs respectively. In a re- 3900 cent Monte Carlo study of larger graphene sheets [18], height fluctuations with a typical size of 0.7 Å have been 3600 observed at room temperature. In the present study room temperature Dh values for GNRs are always larger than 3300 this value and varying approximately from 1.0 to 4.0 Å. 3000 Finite size effects may be responsible for these higher val- ues of hight fluctuations of GNRs. 0 2 4 6 8 10 In conclusion, we have observed a complicated relation width (hexagonal rings) between stability and width of GNRs by MD simulations, Figure 4 Width versus final temperature graph for zigzag edged rotational freedom being the stabilizing factor for thin GNRs. GNRs. GNRs of a width of two hexagonal rings turn out to www.pss-b.com © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim  pss b solidi physica status 700 N. Dugan and Ş. Erkoç: Stability analysis of graphene nanoribbons Table 4 Distance ( Dh ) between the atoms which are most sepa- [4] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 rated from average GNR plane at two sides for zigzag-edged (2005). GNRs of different widths. [5] M. I. Katsnelson, Eur. Phys. J. B 51, 157 (2006). [6] C. L. Kane, Nature 438, 168 (2005). width final T Dh (1 K) Dh (300 K) Dh (final T) [7] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, 1 3600 0.0006 1.0485 3.5386 Science 306, 666 (2004). 2 4800 0.0040 1.3488 2.6138 [8] P. L. McEuen, Nature 393, 15 (1998). 3 4300 0.0197 2.0695 1.9402 [9] Z. Yao, H. W. Ch. Postma, L. Balents, and C. Dekker, Na- 4 4100 0.0151 0.9255 3.6811 ture 402, 273 (1999). 5 4300 0.0267 2.3399 4.8865 [10] Graphene’s Unique Properties Offer Much Potential, APS 6 3200 0.0305 0.8413 3.1611 News 15(5), 1 (2006). 7 3500 0.0426 2.2247 4.4047 [11] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. 8 3000 0.0458 2.7973 5.1321 Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Natl. 9 3200 0.0524 3.9193 3.5454 Acad. Sci. 102, 10451 (2005). [12] K. Kusakabe and M. Maruyama, Phys. Rev. B 67, 092406 be the most stable whereas absence of rotational freedom (2003). makes the wider GNRs to be less stable for both armchair [13] L. Brey and H. A. Fertig, Phys. Rev. B 73, 195408 and zigzag-edged cases. Evenness/oddness of the width (2006). causes an alternating stability behavior for zigzag-edged [14] K. Nakada, M. Fujita, G. Dresselhaus, and M. S. Dressel- GNRs. haus, Phys. Rev. B 54, 17954 (1996). [15] K. Wakayabashi, M. Fujita, H. Ajiki, and M. Sigrist, Phys. Acknowledgements The authors would like to thank Rev. B 59, 8271 (1999). METU and TUBITAK for partial support through the projects [16] M. Ezawa, Phys. Rev. B 73, 045432 (2006). METU-BAP-2006-07-02-00-01 and TUBITAK-TBAG-107T142, [17] Y. Son, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 97, respectively. One of the authors (N.D.) would like to thank Turk- 216803 (2006). ish Petroleum Foundation for partial financial support. The au- [18] A. Fasolino, J. H. Los, and M. I. Katsnelson, Nature Mater. thors would like to thank Dr. Hande Üstünel for critically reading 6, 858 (2007). the manuscript. [19] J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature 446, 60 (2005). References [20] J. Tersoff, Phys. Rev. Lett. 61, 2879 (1988). [21] J. Tersoff, Phys. Rev. B 38, 9902 (1998). [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. [22] L. Verlet, Phys. Rev. 159, 98 (1967). Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Fir- [23] D. W. Heermann, Computer Simulation Methods in Theore- sov, Nature 438, 197 (2005). tical Physics (Springer Verlag, Berlin, Heidelberg, 1990), [2] Y. Zhang, Y. Tan, H. L. Stormer, and P. Kim, Nature 438, p. 35. 201 (2005). [24] L. V. Woodcock, Chem. Phys. Lett. 10, 257 (1971). [3] N. M. R. Peres, F. Guinea, and A. H. Castro Neto, Ann. [25] V. Barone, O. Hod, and G. E. Scuseria, Nano Lett. 6, 2748 Phys. 321, 1559 (2006). (2006). © 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim www.pss-b.com