Table Of Content

Vibrational properties of nano-graphene Sandeep Kumar Singh and F. M. Peeters Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerpen, Belgium (Dated: January 29, 2013) Theeigenmodesandthevibrationaldensityofstatesofthegroundstateconfigurationofgraphene clusters are calculated using atomistic simulations. The modified Brenner potential is used to describe the carbon-carbon interaction and carbon-hydrogen interaction in case of H-passivated edges. ForagivenconfigurationoftheC-atomstheeigenvectorsandeigenfrequenciesofthenormal modes are obtained after diagonalisation of the dynamical matrix whose elements are the second derivative of the potential energy. The compressional and shear properties are obtained from the divergence and rotation of the velocity field. For symmetric and defective clusters with pentagon 3 arrangement on the edge, the highest frequency modes are shear modes. The specific heat of the 1 clustersisalsocalculatedwithintheharmonicapproximationandtheconvergencetotheresultfor 0 bulk graphene is investigated. 2 n PACSnumbers: 63.22.Kn,63.22.Rc,65.80.Ck a J 8 I. INTRODUCTION calandthermodynamicalpropertiesofcarbonallotropes 2 were computed by Mounet et al.25 using a combination ofdensity-functionaltheorytotal-energycalculationsand During the past few years considerable attention has ] density-functional perturbation theory lattice dynamics i been paid to the study of the spectral properties such c in the generalized gradient approximation. Very good as the energy spectrum, the eigenmodes (i.e. phonon s agreement was found for the structural properties and - spectrum), and the phonon density of states of car- l the phonon dispersion. r bonallotropessuchasgraphene1–4,carbonnanotubes5,6, t m fullerenes7–9, graphene flakes10 and carbon clusters11,12. High-level ab initio calculations were carried out to Ramanspectroscopyisawidelyusedtooltostudythe . reexamine the relative stability of bowl, cage, and ring at vibrational modes of carbon based materials13–15. Be- isomersofC andC− byAnetal.26 Thetotalelectronic m cause the normal modes of molecules are unique, they 20 20 energies of the three isomers showed a different energy have their own set of Raman frequencies. Similarly for - ordering which depends on the used hybrid functional. d graphene clusters their phonon modes will give us infor- n mationone.g. theirstructuralanddynamicalproperties. Finite size planar structures which are close to the o Lattice dynamics calculations of a single graphene ground state, and in particular nanographene-like struc- c [ sheet were done by Rao et al.16 based on C-C force tureswerestudiedbyKosimovetal.27 usingenergymini- constants17,18 which were obtained from a fit to the mizationwiththeconjugategradient(CG)method. The 1 two dimensional, experimental phonon dispersion. Good lowest energy, i.e. the ground state27,28 configurations v agreement with Raman spectra of single wall car- were determined for up to N=55 carbon atoms. 7 bon nanotubes was found. Tommasini et al.19 pre- 3 5 sented a semiempirical approach for modeling the off- Here we will investigate the normal modes for the 6 resonance Raman scattering of molecular models of con- ground state configuration of flat carbon clusters, also . fined graphene and compared with the results from den- called nanographene as function of the number of C 1 0 sity functional theory calculations. Pimenta et al.20 re- atoms in the cluster. The content of the paper is as 3 viewed the defect-induced Raman spectra of graphitic follows. First (see Sec. 2), we present the main compu- 1 materialsandpresentedRamanresultsonnanographites tational approach that is used to calculate the normal : and graphenes. Ryabenko et al.21 analyzed the spectral modes. The normal modes for linear (N ≤ 5) and ring v i properties of single-wall carbon nanotubes encapsulat- (6 ≤ N ≤ 18) clusters which are exactly one and two X ing fullerenes and found that fullerene alters the aver- dimensional, respectively are discussed in Sec. 3. Next r agediameteroftheelectroncloudaroundthesinglewall (see Sec. 4), the phonon spectrum of two dimensional a carbon nanotube. Michel et al.22 presented a unified nanographene (19≤ N ≤ 55) and trigonal-shaped nan- theory of the phonon dispersion and elastic properties odisk29,30 as well as the effect of hydrogenation of the of graphene, graphite and graphene multilayer systems. edge carbon atoms on the phonon spectrum are investi- Theystartedfromafifth-nearest-neighborforce-constant gated. InSec.5,wecalculatethephonondensityofstates model derived from the full in-plane phonon dispersions ofvariousclusters. ByincreasingthenumberofCatoms of graphite23, and extended this model with interplanar in the cluster we will show how the phonon spectrum interactions, and found that the graphite eigenfrequency of graphene appears. We also calculate the temperature ω ≈ 127 cm−1 is reached within a few percent for dependence of the specific heat within the harmonic ap- B2g1 N=10 layers which agrees with similar results obtained proximation in Sec. 6. Summary (see Sec. 7) close the for the electron spectrum24. The structural, dynami- paper. 2 II. SIMULATION METHOD Using molecular dynamic (MD) simulations with the Berendsenaftercorrectorthermostat,weobtaintheequi- TheBrennersecondgenerationreactiveempiricalbond librium configuration of stable quasi-planer clusters. We order (REBO) potential function between carbon atoms found that about 104-5×105 MD steps are needed to is used in the present work. The values for all the pa- obtain an accuracy of 10−5-10−6 eV in the energy per rametersusedinourcalculationfortheBrennerpotential particle. We used free boundary conditions in our simu- canbefoundinRef.31andarethereforenotlistedhere. lation. The Brenner potential (REBO) energy is given by: Here, we calculate the dynamical matrix numerically usingthefinite-displacementmethod. Wedisplaceparti- (cid:88) (cid:88) cle j in the direction β by a small amount ±u, and eval- Eb = [VR(rij)−bijVA(rij)]. (1) uate the forces on every particle in the system F± (we α,i i j(>i) took typically u=0.003 ˚A). Then, we compute numerically the derivative using the central-difference formula HereE isthetotalbindingenergy,VR(r )andVA(r ) b ij ij are a repulsive and an attractive term, respectively, whererij isthedistancebetweenatomsiandj,givenby dF F+ −F− d2E α,i = lim α,i α,i = b . (7) duβ,j u→0 2uβ,j duα,iduβ,j VR(r)=fC(r)(1+Q/r)Ae−αr, (2) Diagonalization of the dynamical matrix yields the eigenvalues and the eigenvectors from which we can cal- 3 VA(r)=fC(r)(cid:88)B e−βnr, (3) culate the vibrational modes. The eigenfrequencies are n the square root of the eigenvalues. The stability of the n=1 configurationcanbeverifiedbecausealleigenvalueshave wherethecut-offfunctionfC(r)istakenfromtheswitch- to be real. ing cutoff scheme 1 r <Dmin fC(r)=(cid:2)1+cos(cid:0) (r−Dimjin) (cid:1)(cid:3)/2 Dmini<j r <Dmax ij (Dmax−Dmin) ij ij 0 ij ij r >Dmax ij (4) where Dmax −Dmin defines the distance over which the ij ij function goes from one to zero and A,Q,α,B ,β ,(1 ≤ n n n≤3) are parameters for the carbon-carbon pair terms. Here n is the type of bond (i.e. single, double and triple bonds). Theempiricalbondorderfunctionusedinthisworkis written as a sum of terms: 1 b = [bσ−π+bσ−π]+bπ, (5) ij 2 ij ji ij where the functions bσ−π and bσ−π depend on the local ij ji positionandbondanglesdeterminedfromtheirarrange- ment around each atom (i and j, respectively) and is governed by the hybridization of the orbitals around the FIG.1: (Coloronline)Spectrumofnormalmodesasafunc- tionofthenumberofparticlesintheclusterforlinearN≤5 atom. The first term in Eq. (5) is given by, andring6≤N ≤18clusters. Thesolidcurvesaretheanalyt- icalresultsfromasimplechain(N≤5)andring(6≤N ≤18) bσ−π =(cid:2)1+ (cid:88) fC(r )G(cos(θ ))(cid:3)−1/2. (6) model. The dashed curve indicates the breathing mode. The ij ik ik ijk symbols are the numerical results from the present simula- k((cid:54)=i,j) tionwherethesolid(hollow)symbolsarethefrequencieswith correspondingeigenvectorswhichareinplane(outofplane). Here the angular function G(cos(θ )) modulates the ijk The frequencies that decrease with the number of particles contributionofeachnearestneighbourandisdetermined are fitted to a 1/N dependence (solid curves). bythecosineoftheangleofthebondsbetweentheatoms i−j−k. 3 direction to each other along the chain as depicted in Fig. 3. While the cluster N=3 with mode k=9 and the cluster N=5 with k=15 have alternating particles moving in opposite direction which makes an angle with the chain direction (see Fig. 3). The breathing mode for the linear chain is shown by the dotted line in Fig. 1. For N=2 the excitation corresponding to the breathing mode is shown in Fig. 3(a). But for N=3, this mode is similar for N=2 but where thedirectionofoscillationmakesananglewiththechain direction and the middle particle does not move due to conservationofmomentum. Similarly, forN=5wefound that the breathing mode is similar to the one for N=4 except that the deflection makes an angle with the chain directionandthemiddleparticleisatrest(seeFigs.4(a, b)). For mode k=6 and 7 with N=3,4,5 particles, the particles on the boundaries of the clusters move in the FIG. 2: Fitting parameter f of linear chain (N ≤ 5) and same direction while the remaining particles move in the ring (6 ≤ N ≤ 18) model of particles connected by identical opposite direction with different amplitudes to conserve springs as a function of the number of particles. For N ≥ 6 totalangularmomentum(seeFig.4(c))forbothin-plane theresultsarefittedbyasimpleanalyticalexpressionforeven andout-of-planedirection. Thecorrespondingfrequency (full curve) and odd (dashed curve) clusters. of these modes can be reasonably well fitted by a+b N (solid curve in left panel of Fig. 1) where the fitting parameters are (a,b)= (1809 (±103), -234 (±27)) cm−1. III. LINEAR AND RING CLUSTERS Andformodek=8and9withN=4,5particles,adjacent particlesmoveinoppositedirectionwithdifferentampli- The energy of the normal modes corresponding to the tudes where the central particle is at rest in the case of ground-state configuration for linear and ring clusters is N=5 particles (see Fig. 4(d)) for both in plane and out shown in Fig. 1 as function of the number of C atoms in of plane direction and they can be fitted by the same theclusters. ThenumberofnormalmodesforN-particles function with (a,b)= (1783 (±103), 0.0) cm−1. moving in 3-dimensions is 3N. Some of these modes will be degenerate. There are 5 modes with ω = 0 which correspond to 3 uniform translations and 2 independent rotations of the whole cluster, due to the translational and rotational invariance of the cluster. Whenclustersareexactlyonedimensional,thephonon modes correspond to oscillations along the chain or or- thogonal to it. The normal modes of such a chain of N particlesconnectedbyspringswithspringconstantκand mass m can be calculated analytically32 and is given by ω =f×sin(φn/2) wheref =2(cid:112)κ,andφ =(n−1)π/N sinφN m n (n = 1,2,··· ,N). For each value of N, we determined f by taking the maximum frequency obtained from our model equal to our numerically found result. The ob- tainedfrequenciesofourlinearspringchainmodel(solid curves in Fig. 1) agree remarkably well with our numerical results. For N ≤ 5, the obtained spring constant is plotted in Fig. 2 and is almost constant (less than 4% variation). As shown in Ref.28 all the C-C bonds in FIG. 3: Eigenvectors of the highest frequency modes for the chain are not identical and therefore all the spring the clusters with N=2,3,4 and 5 particles for different mode constants between the C-atoms are expected not to be number with frequency (a) ω6 = 1660 cm−1, (b) ω9 = 2092 identical, whichisthereasonwhyoursimplemodeldoes cm−1, (c) ω12 =2188 cm−1 and (d) ω15 =2233 cm−1. not give perfect agreement with our numerical data. Notice that the low frequency eigenmodes are not de- The ground state configuration for C (N=6-18) con- N scribed by this simple model. Some of them correspond sists of particles arranged in a ring. For clusters that to out-of-plane motion. For cluster N=2 with mode k=6 are exactly two dimensional, as is the case of these ring (k counts the eigenvalues in increasing order) and for structures, the phonon modes correspond to motion in N=4 with mode k=12, the particles move in opposite the plane or vertical to the plane. There is no cou- 4 FIG. 4: Eigenvectors for the clusters with N=4 and 5 particles for different mode number with frequency (a) ω =905 10 cm−1, (b) ω = 725 cm−1, (c) ω = 150 cm−1 and (d) 12 6 ω =466 cm−1. 9 FIG. 5: Eigenvectors for a cluster with N=18 particles for pling between the in-plane and the out-of-plane motion. different mode number with frequency (a) ω = 398 cm−1, 23 We can construct a simple model of equidistant par- (b)ω =550cm−1,(c)ω =2224cm−1 and(d)ω =2259 38 53 54 ticles arranged in a ring, where neighbor particles are cm−1. connected by a spring. This model can be viewed as the linear chain model with periodic boundaries. For fixed radius of the ring this results in the following spondstonearest-neighboratomsoscillatingoutofphase eigenfrequencies33 ω = f ×sin(k/2), where k = 2πn/N, with different amplitudes except for the two neighbor (n = 1,2,··· ,N). As before we determined f such that atoms that oscillate in phase (see Fig. 6(b) k=23 (471 the maximum frequency obtained from the model equals cm−1))whichresultsinalowerfrequencybyabout20% the numerical result. This fitted value is shown in Fig. 2 as compared to the corresponding even N ring clusters. which can be described approximately by the function A breathing mode exists for all the ring clusters as f(N)=a×(1+bN), where: a=3200 cm−1 and i) for even shown in Fig. 6. The frequencies for N=8, 12, 13, 16, (1+cN) Nclustersb=-0.295(±0.001),andc=-0.395(±0.002); and 18 are 837 cm−1, 585 cm−1, 543 cm−1, 447 cm−1, and ii) for odd N clusters b = -0.276 (± 0.001), and c = and 398.25 cm−1, respectively, and have a clear depen- -0.369(±0.001). Severalofthenumericalfrequenciesare denceonN.Thesimplespringmodel34 predictsthatthis reasonably well described by this simple model. Notice breathing mode has the frequency ω=(π ×f)/N which that the eigenfrequencies exhibit a clear discontinuous remarkably agrees with our numerical results. behaviorwhenaclusterchangesfromalinearintoaring The eigenmodes corresponding to ω have similar max configuration. As an example we show in Fig. 5 some of vibrationalmodesforevennumberofparticles. Weplot- the interesting eigenmodes of a cluster with N=18 parti- ted the normal modes of ω of cyclic structures for max cles which are arranged in a circle and is therefore pure N=12 (ω =2183 cm−1) as shown in Fig. 6(e) and they 36 two dimensional. The mode k=23 is the breathing mode correspond to dipole-type of oscillations between nearest while mode k=38 corresponds to radial out of phase ra- neighbor particles while for odd number of ring struc- dial oscillation of nearest neighbor particles. The modes tures i.e. N=13 (ω = 2187 cm−1 in Fig. 6(f)), similar 39 k=53, and 54 correspond to pure angular oscillations of dipole-type of oscillations between nearest neighbor par- the particles. ticles are found but with decreasing magnitude towards Itisinterestingtonotethatforanevennumberofpar- opposite sides. ticles,thereisanearlycommonfrequency(i.e. N=8(545 Nowweinvestigatetheeigenmodesforout-of-planevi- cm−1),N=10(543cm−1),N=12(546cm−1),N=14(547 brationsforringclusters. ForN=10andmodek=10,the cm−1), N=16 (549 cm−1), N=18 (550 cm−1)) which is normal mode corresponds to a bending mode while for almost independent of N (see Fig. 1). This mode corre- modek=16,thenormalmodecorrespondstoasinusoidal sponds to out of phase oscillations of nearest-neighbor type of motion. The higher modes k=17 and k=19 have atoms (see Fig. 5(b) the k=38 mode for N=18 and nearest-neighbor atoms oscillating in opposite direction Fig. 6(a) the k=24 mode for N=12). But for odd ring perpendiculartothering-planewithdifferentwavelength clustersi.e. N=13,thecorrespondingnormalmodecorre- along the ring as seen in Fig. 7. For even ring structures 5 FIG. 7: Eigenvectors for out-of-plane motion for a cluster withN=10particlesfordifferentmodenumberwithfrequency (a)ω =165cm−1,(b)ω =482cm−1,(c)ω =482cm−1 10 16 17 and (d) ω =536 cm−1. 19 FIG. 6: Normal modes for the out-of-phase (top panels), breathing (middle panels) and ω (bottom panels) for the max clusters N=12 (left panel) and N=13 (right panel). TABLEI:Fittingparametersfortheeigenfrequenciesofout- of-planemotionshownbythesolidcurvesinFig.1forω<550 cm−1: ω= a +b. FIG. 8: 1/DC−C (dashed curve) and ωmax (solid curve for N even and odd N) as a function of the number of particles. a cm−1 b cm−1 2782 (±48) -107 (±5) FromFig.1,itisinterestingtonotethatthemaximum 4258 (±120) -97 (±11) frequency in the excitation spectrum for ring clusters, 5405 (±251) -62 (±20) on the average, slowly increases with N. This can be ex- plainedbycalculatingthephononspectrumofaninfinite 5686 (±59) 25 (±5) system. InFig.8weshowthattheminimalinterparticle 5814 (±499) 102 (±33) distance decreases slowly with the number of C-atoms 4942 (±589) 220 (±36) in the cluster. As a consequence, the maximum value of the wave vector k ≈π/l (l is the mean distance be- 0 0 tweentheparticles)andalsothecorrespondingfrequency the highest mode for out-of-plane vibration is similar to will increase weakly with cluster size. The inverse of the the one shown in Fig. 7(d). The out-of-plane vibration inter-carbon distance could be fitted to (dashed curve in can be reasonably well fitted by a +b (solid curves in Fig. 8) N Fig. 1 for ω<550 cm−1) where the fitting parameters a and b are listed in Table I. The eigenmodes belonging to 1 1 1+aN = (8) each solid curve correspond to the same type of mode. D D 1+bN C−C 0 6 with the fitting parameters: 1 =7.9 nm−1, a = -0.316 D0 (± 0.002), and b = -0.331 (± 0.002). We are able to fit a curve through the maximum frequency of these ring clusters using (solid curves in Fig. 8) 1+aN ω =ω , (9) 01+bN where ω =3200 cm−1, and: i) for even N clusters a = 0 -0.295 (± 0.001), and b = -0.395 (± 0.002); and ii) for odd N clusters a = -0.276 (± 0.001), and b = -0.369 (± 0.001). IV. NANOGRAPHENE (N >18) The clusters with N >18 have an inner structure and their configurations have been investigated in detail in FIG. 9: Eigenvectors for a cluster with N=20 particles for Ref. 27. Their ground state configuration can be clas- different mode number with frequency (a) ω = 124 cm−1, sified in three groups: 1) nanographene clusters consist- 7 (b) ω =161 cm−1, (c) ω =233 cm−1 and (d) ω =1772 ing of only hexagons, 2) clusters with pentagon on the 9 11 60 cm−1. boundary, and 3) bowl shaped configurations that have typically pentagons in the inner part of the cluster. For the latter one < z2 >(cid:54)= 0, where z is the position coordinate of the cluster along the out-of-plane direction and they are found for N=20, 28, 38 and 44 which are buckled-likestructures. Intheseclustersonepentagonis surrounded by five hexagons. The normal mode oscillations in-plane and out-of-plane are now coupled, i.e. the normal mode oscillations are 3-dimensional and some of the interesting ones are shown in Fig. 9 for N=20 and mode k=7, 9, 11 and 60. For mode k=7, the opposite atomsontheboundaryoftheclustervibrateinthesame directionwhereasfork=9, thenormalmodecorresponds to a bending mode of the cluster. But, for the higher mode k=11, the corner atoms show mixed type of oscillations while for large value of k, the normal modes correspond to in-plane oscillations of the C-atoms. ThenormalmodefrequenciesforN ≥19areplottedin Fig.10. ForN=22and39,aheptagonisontheboundary surrounded with pentagon and hexagon. In such a case local modes are found where only particles close to the FIG.10: (Coloronline)Normalmodespectrumasafunction defects (pentagon and heptagon) oscillate. If there is of the number of C atoms for 19≤N ≤55. The modes with one pentagon on the boundary i.e. N=49, and 51, the eigenvectors in-plane are shown in black, those with out-of- local modes with larger amplitudes are on the opposite planeeigenvectorsinredandmixedtypeeigenvectorsinblue. side of the clusters while for N=53 where a pentagon is surroundedbyfourhexagons,thelocalmodeswithlarger amplitudes are found only near the defect. For clusters atoms arranged in hexagons near the 6 corners of the which are symmetric and without defect i.e. N=54, we hexagonal disk. The next mode k=59 is similar but now foundthatthelocalmodesaresituatedontheboundary hexagons are displaced by an angle of 30◦ and are cen- as shown in Fig. 11. tered in the middle of the sides of the hexagonal disk. A special subgroup of clusters consists of hexagonal Notice that for k=58 the 6 rotational oscillations are in and trigonal shaped nanographene. The clusters with the same direction while for k=59 they alternate, i.e. size N=24 and N=54 are planar graphene structures vortex/anti-vortex like arrangements. For higher modes which are hexagonal shaped with zigzag edges. These k=139, onlytheinnerparticlesparticipateinthenormal clusters have a close-packed structure, which consist mode oscillations. For k=148, a dipole type of oscilla- purely of hexagons. Let us consider the cluster with tions is found of nearest neighbor C-atoms arranged in N=54 (see Fig. 12), the mode k=34 corresponds to the a shell around the middle between the center and the breathing mode. Mode k=58 exhibits circular motion of perimeter. For higher frequency only the outer particles 7 FIG. 12: Eigenvectors for a cluster with N=54 particles for different mode number with frequency (a) ω = 313 cm−1, 34 (b)ω =531cm−1,(c)ω =534cm−1,(d)ω =593cm−1, 58 59 70 (e)ω =1589cm−1,(f)ω =1624cm−1,(g)ω =1628 139 148 149 cm−1, (h) ω =1638 cm−1 and (i) ω =1914 cm−1. 150 162 situated close to the three corners of the trigonal nanodisk, but notice that the rotation direction and the po- sitionofthecenteroftherotationisnotalwaysthesame for both modes . FIG. 11: Eigenvectors for highest mode number for clusters We plotted the average distance between the C-C with N=22, 39, 49,51, 53, 54 particles with frequency (a) ω = 1944 cm−1, (b) ω = 1926 cm−1, (c) ω = 1907 atoms as a function of the number of particles N. The 66 117 147 cm−1, (d) ω = 1903 cm−1, (e) ω = 1952 cm−1 and (f) averageradiusincreaseslinearlywiththechainlengthfor 153 159 ω162 =1914 cm−1. (N=3-5) which could be fitted to DC−C = D0 +a×N whereD =0.13nmanda=0.00055(±0.00009)nm(see 0 solid curve in Fig. 15). For ring clusters it decreases ex- oscillate while the inner particles exhibit only very small ponentially as D = D +a×exp(b×N) (see solid C−C 0 displacements as shown in Fig. 12. curve in Fig. 15) where the parameter D = 0.134 nm, 0 FromFig.10wenoticethatthereisaregionalongthe a = 0.039(±0.001) nm and b = −0.337(±0.007) as the frequency axis with a higher density of normal modes. numberofparticlesincreases. ForN >18theaverageC- Examples of the displacements of such modes are shown CdistancefluctuatesasfunctionofNaroundtheaverage in Fig. 13 for N=54. Mode k=27 corresponds to the ex- value0.141nmwhichcompareswith0.142nmfortheC- citationofavortex/antivortexpair. Modek=33consists Cdistanceinbulkgraphene. Theclusterswithpentagon of rotational oscillations around the hexagonal corners and Stone-Wales defects inside the clusters show larger of the nanodisk, while for k=35 these rotational motions averageC-Cdistancethanthepurehexagonalstructures are centered around the middle of the hexagonal sides. while the clusters with pentagons on the boundary (i.e. Notice that mode k=45 consists of a central large vor- N=26, 31, 34, 36, 41, 46, 51 and 53) show the highest tex motion surrounded by 6 anti-vortex type of motions average distance as shown in Fig. 15. situated closer to the edge of the nanodisk. For completeness we also plotted (in Fig. 16) the en- The nanodisk cluster with N=46 carbon atoms has ergy spectrum of the normal modes for ground-state H- a metastable configuration consisting of closed zigzag passivated C (N=11-55) clusters27. Notice that dif- N edges and a trigonal-shaped structure27. As shown in ferent from Fig. 10 (i.e. non H-passivated graphene) Fig. 14 the mode k=18 consists of an asymmetric vor- there is a region with enhanced density of modes, i.e. tex/antivortexpairwhilemodek=29isaclearbreathing around 2900 cm−1. The corresponding modes are con- mode. The modes k=24 and k=32 show three rotations nected to oscillations of the C-H atoms. This is illus- 8 FIG. 13: Eigenvectors for the cluster with N=54 particles FIG. 14: Eigenvectors for a cluster with N=46 particles for forfourdifferentvaluesofthemodenumberkwithfrequency different mode number with frequency (a) ω = 185 cm−1, 18 (a)ω =244cm−1,(b)ω =302cm−1,(c)ω =314cm−1 (b) ω = 244 cm−1, (c) ω = 288 cm−1 and (d) ω = 322 27 33 35 24 29 32 and (d) ω =408 cm−1. cm−1. 45 trated in Fig. 17 for H-passivated nanographene clusters with N=53 and 54. In the optical region which is a region along the frequency axis with higher density of nor- malmodes,forN=53withk=172thecorrespondinglocal normalmodeswithlargeramplitudesareontheopposite side of the defect while for N=54 with k=187, only the inner C-atoms oscillate with larger amplitude (Figs. 17 (a, b)). The highest frequency mode of the N=53 and 54 clusters (Figs. 17 (c, d)) is a local mode near the lo- caldefectwheremainlyouterC-Hatomsofthearmchair C-atoms oscillate. V. PHONON DENSITY OF STATES In order to compare the normal modes with those of graphene we calculate the phonon density of states FIG. 15: The average distance between C-C atoms as a (PDOS). Furthermore, we calculated the PDOS of ex- function of the number of particles in the graphene cluster. actly trigonal and hexagonal two dimensional and de- TheerrorbarsindicatetherangeofC-Cdistanceswithineach fective clusters and analyze what is the effect of defects cluster. on the density of states. Because of the discreteness of the normal modes frequency spectrum we introduced a Gaussian broadening in the cluster. We notice that the introduction of a pentagon in the N=53 cluster as compared to the perfect 3N (cid:88) hexagon lattice in N=54 introduces more high frequency ρ(ω)= exp(−(ω−ω )2/δω2), (10) i modesasshowninFig.18whichareduetomodeslocal- i=1 ized around the defect. where the summation is over all normal modes, ω is the Fig. 19 shows the phonon density of states (PDOS) of i normal mode frequency of the ith mode, broadening is armchair hexagonal and zigzag trigonal clusters. Notice chosenδω =30cm−1 andNisthetotalnumberofatoms that several of the peaks in the PDOS of the hexagonal 9 FIG. 18: (Color online) The phonon density of states for clusters with N= 53 and 54 carbon atoms. FIG.16: (Coloronline)Normalmodespectrumasafunction of the number of C particles for 11 ≤ N ≤ 55 H-passivated clusters. The modes with eigenvectors in plane are shown in black, those having out-of plane eigenvectors in red. FIG. 19: (Color online) The phonon density of states for clusterswithN=42(hexagonalwitharmchair)and46(trigonal) carbon atoms. FIG. 17: Eigenvectors for a cluster with N=53, and 54 with H-passivation for k=172 for N=53 and k=187 for N=54, and for ω corresponding to frequency (a) ω = 1693 cm−1, The effect of H-passivation of the edge atoms of the max 172 (b) ω =1698 cm−1, (c) ω =2925 cm−1 and (d) ω = cluster is shown in Fig. 21 for N=53 and 55. The C- 187 210 222 2915 cm−1. H bonds introduce high frequency modes in the N=53 and N=55 clusters as compared to the non H-passivated N=53 and N=55 clusters. The H-passivated clusters ex- clustersaresplitintotwopeaksinthePDOSofthetrig- hibit pronounced peaks around 1700 and 2900 cm−1 onalshapedcluster. Thisisaconsequenceofthereduced whichcorrespondtotheopticalregionandtheC-Hatom rotational symmetry of the trigonal cluster. frequency region, respectively. We also notice that the We compare in Fig. 20 the phonon density of states of introduction of a pentagon in the N=53 cluster as com- clusters N=55 and 1600 with those for graphene35 which paredtotheperfecthexagonlatticeforN=55introduces shows nicely the convergence of the phonon spectrum of morehighfrequencymodesasshowninFig.21similarto nanographene to bulk graphene. For N=1600 we almost thecasefornonH-passivatedclusters(seeFig.18)which recoverthetheoreticalresultsofgraphene35. Noticethat are due to modes localized around the defect. for small clusters there is a pronounced PDOS for ω =0 The modes can have a shear-like or a compression-like which is due to the relative importance of the ω = 0 character. The compressional and shear properties can translational and rotational motion. be extracted from the divergence and rotation of the ve- 10 N 1 (cid:88) Ψ (k)= Ψ2 (k), (12) r N r,i i=1 wherethevaluesofΨ (k)andΨ (k)fortheithparticle d,i r,i are, respectively, given by J Ψ (k)= 1 (cid:88)((cid:126)r −(cid:126)r ).[A(cid:126) (k)−A(cid:126) (k)]/|(cid:126)r −(cid:126)r |2, (13) d,i J i j i j i j j=1 J Ψ (k)= 1 (cid:88)((cid:126)r −(cid:126)r )×[A(cid:126) (k)−A(cid:126) (k)]/|(cid:126)r −(cid:126)r |2. (14) r,i J i j i j i j j=1 Here, j and J denote the index and the number of neighboring particles of particle i, respectively. (cid:126)r is the posi- j tional coordinate of a neighboring particle and A(cid:126) (k) is i the eigenvector of particle i for mode k. Note also that we calculate the squared average over all the particles because the simple spatial average is of course zero. We plot Ψ (k) and Ψ (k) as a function of the mode FIG. 20: (Color online) The phonon density of states for d r number k for cluster N=39 which contains a 5-7 defect, clusterswithN=55,1600andperfectgraphene(bottompanel is taken from Ref. 35). cluster N=49 with pentagon defect on the boundary, the pure symmetric hexagonal cluster with N=54 and for a large cluster having N=398 particles. In general, the lower eigenfrequency spectrum corresponds to rotational type of excitations which are vortex-antivortex like excitations for large clusters (see Fig. 13). Vortex excitation is only expected for sufficiently large clusters, because the velocity of dissipation of the vortex energy is inverselyproportionaltoR2,whereRistheradius,which increases with increasing cluster size. In the second half ofthespectrumthedivergenceΨ (k),whichcorresponds d to compression-like modes, is appreciably different from zeroandwehavemixedmodesthathavebothashear-like and a compression-like component. For symmetric clusters and defective clusters with pentagon on the boundary, the highest modes consist of only shear-like modes (seeFigs.11(c,d,f)). ForclusterN=39witha5-7defect, compression-like modes appear for the highest frequency modes (see Fig. 11(b)). The heat content of a system is directly related to the FIG. 21: (Color online) The phonon density of states for specific heat whose lattice contribution is determined by H-passivated clusters with N=53 and 55. the phonons i.e. the eigenmodes of the system, while electronic contribution to it can be neglected even at a few Kelvins37. The specific heat C per unit mass, at P locity field. Here, we will associate a single number for temperature T is within the harmonic approximation38 the shear-like and compression-like character of the dif- given by the following expression: ferent modes by calculating the squared average of the divergence∇(cid:126).(cid:126)υ andthevorticity(∇(cid:126) ×(cid:126)υ)z ofthevelocity C = kB (cid:88)3N (cid:18) (cid:126)ωk (cid:19)2csch2(cid:18) (cid:126)ωk (cid:19), (15) field, following the approach of Ref. 36. P MN 2k T 2k T B B k=1 The z component of the rotation Ψ (k) = (cid:126)e .rotΨ(k) r z and the divergence Ψd(k) = divΨ(k) of the field of the where the summation is over all normal modes, kB is eigenvectors of mode k are the Boltzmann constant, M is the mass of each carbon atom and N is the total number of atoms in the cluster. From Eq. (15), the specific heat depends sensitively on N Ψ (k)= 1 (cid:88)Ψ2 (k), (11) the characteristics of the phonon spectrum and on its d N d,i vibrational density of states. i=1