Vibrational stability of graphene

The mechanical stability of graphene as temperature rises is analyzed based on three different self-consistent phonon (SCP) models. Compared with three-dimensional (3-D) materials, the critical temperature Ti at which instability occurs for graphene is much closer to its melting temperature Tm obtained from Monte Carlo simulation (Ti ≃ 2Tm, K. V. Zakharchenko, A. Fasolino, J. H. Los, and M. I. Katsnelson, J. Phys. Condens. Matter 23, 202202). This suggests that thermal vibration plays a significant role in melting of graphene while melting for 3-D materials is often dominated by topologic defects. This peculiar property of graphene derives from its high structural anisotropy, which is characterized by the vibrational anisotropic coefficient (VAC), defined upon its Lindermann ratios in different directions. For any carbon based material with a graphene-like structure, the VAC value must be smaller than 5.4 to maintain its stability. It is also found that the high VAC value of graphene is responsible for it...

The experimental success in graphene 1 initiated a recent revival of interest in two-dimensional (2-D) solid state materials due to its peculiar physical properties and application prospect. 2,3 Yet the mechanism of thermo-stability for this type of new material is not understood. In fact, the existence and stability of low dimensional crystalline material has long been a fascinating scientific problem. Peierls and Landau 4, 5 proved long ago that 1-D and 2-D crystalline state cannot exist for materials made up of particles interacting with harmonic potentials. Their conclusion was later extended by Mermin 6 on a rigorous basis using the Bogoliubov's inequality from statistical physics. Mermin's work supported the result that long range crystalline order cannot persist for low-D solids with particles interacting with harmonic pair potentials. The dilemma between the aforementioned classic theoretical works (hereafter referred to as PLM theory) and the existence of graphene can be explained as follow: the materials considered in PLM theory is constrained in the low-D space where they are defined (e.g. for a 2-D solid discussed in PLM theory, the out-of-plane motion of particles is forbidden). Suspended graphene violates the assumption of PLM theory by two points: a). the out-of-plane motion is not forbidden. On the contrary, since there is only one layer of atoms, the outof-plane vibrational amplitude is much larger than the in-plane vibrational amplitude; b). the atoms in graphene are not strictly located in a plane. It is proved experimentally 7, 10 and theoretically 8,9 that intrinsic ripples exist for suspended graphene, even at low temperature. 11 Therefore, the PLM theory, although mathematically correct, is not appropriate for suspended graphene.
a Authors to whom correspondence should be addressed. According to Mermin, 6 the physical state of low-D materials is different from both crystalline phase and isotropic liquid phase. Thus melting of low-D materials to isotropic liquid phase is still possible. Kosterlitz, Thouless, Halperin, Nelson, and Young proposed a model (often referred to as the KTHNY theory) 12 aimed at two-dimensional melting. The KTHNY theory suggests that melting for 2-D solids is a two-stage process, with one initiated by dislocation pairs unbinding and the other by disclination pairs unbinding. For graphene, the first step of melting according to the KTHNY theory corresponds to unbinding of the Stone-Wales (SW) defect 21 to two separate 5-7 defects (or called (1,0) dislocations in Ref. 22). Although the SW defect is experimentally observed in graphene 23 and of relatively low formation energy (≈5 ev) 24 among possible topological defects, the unbinding is unlikely to happen. Since the formation energy of one isolated 5-7 defect is calculated as 7.5 ev, 22 replacement of one SW defect by two 5-7 defects is energetically unfavored. Monte Carlo (MC) simulation on melting of graphene 20 indicates a transition from 2-D structure to 1-D chains, where the melting process is initiated by clustering, rather than unbinding, of SW defects. The simulation results might be considered as a support of the grain-boundary based melting theory, 25 although it is questionable whether clusters of SW defects separating 2-D structure and 1-D chains can be called grain-boundaries.
Both theoretical models and simulations on graphene suggest that defects play a significant role in the melting process of 2-D materials, yet another basic question to be understood is the relation between melting and the vibrational stability of 2-D materials. In classic lattice dynamics, the atomic motion due to lattice vibration is assumed to be neglectable compared with the lattice constants, which leads to force constants independent of the vibration. As temperature rises, however, the atomic motion becomes gradually significant which affects the force constants and eventually destroys the mechanical stability of the structure. The temperature at which this instability occurs is called vibrational instability temperature, denoted by T i . For low-D materials such as graphene, the vibrational instability is especially important since the out-of-plane vibrational amplitude can be very large compared with the lattice constants even at low temperature due to its high structural anisotropy. The self-consistent phonon (SCP) theory 13-16 is very effective for analyzing lattice instability due to large amplitude atomic motion. It is established upon the Born-von Karman theory of lattice dynamics, where the potential energy between the atoms is written as a Taylor expansion of the atomic displacements. The basic idea of the SCP theory can be explained as follow: the Taylor series which represent the atomic potential energy (or more specifically the coefficients in the series) are related to the vibrational motion of the atoms in a self-consistent manner. In this letter we extended three self-consistent phonon (SCP) models and calculated the vibrational stability of graphene. The basic equations of the self-consistent harmonic approximation (SCH), which is the lowest order SCP theory, are briefly introduced as follow. The phonon frequency having wave vector q and branch λ is where ε αk (q, λ) is the phonon polarization vector of the kth atom in the unit cell in the α direction, is the dynamical matrix and R l k is the equilibrium position of the k th atom in the l th unit cell. The force constant is averaged over the relative vibrational displacement vector u = u kk (0l ) = u k (0) − u k (l ), where V 0l kk is the interaction potential between the kth atom in the 0th unit cell and the k th atom in the l th unit cell. The width of the Gaussian distribution is determined by Eqs.
(1)-(4) are solved iteratively at given temperature and equilibrium bond length (EBL) until a self-consistent solution of the phonon frequencies is obtained. The SCH model reduces to selfconsistent Einstein approximation (SCE) if the phonon frequency spectrum is replaced by the Einstein frequency (to achieve this, the other atoms are fixed at their equilibrium position when constructing the equations of motion for one atom). Both SCH and SCE models expand the potential to merely quadratic terms. The next-level approximation models where cubic terms are included in the expansion of the potential are called SCHC model and SCEC model. Consideration of the cubic terms in the expansion results in a shift of the phonon frequencies and the phonon life time. This is caused by interaction between phonons which is absent in a harmonic approximation. The phonon frequency in SCHC/ SCEC model is given by where ω qλ is the corresponding phonon frequency obtained in the SCH/SCE model and the expression of (q, λ; ω) can be found in Ref. 13.
In SCP models self-consistent solution may exist for a range of EBL and the free energy of the system can be used to determine the stable EBL at any given temperature. The Helmholtz free energy in the SCH and SCHC models is given by while the energy expression for the SCE and SCEC models is given by where ω i , (i = 1, 2, 3) corresponds to three vibrational frequency of one atom obtained in the SCE model due to material anisotropy and N refers to the number of atoms in the system. The 2-D atomic structure of graphene is illustrated in the x-y plane in Figure 1. There are two carbon atoms in one unit cell (enclosed by dotted lines) and the equilibrium distance between two nearest neighbors is denoted by a.
The C-C potential LCBOPII 17 is used to describe the interaction potential between the carbon atoms. The short range interaction in LCBOPII is a multi-body potential for its short range interaction: The potential between atom 1 and 2 in Figure 1 is related to the bond length r 12 and the bond angles α 312 , α 412 , α 521 , and α 621 . The expansion of the potential between atom 1 and 2 to cubic terms is given by 18 where the potential is a function of the bond length r i j and the cosine of the bond angle cos α i jk , and (X ) 0 denotes the value of X at EBL. r i j and cos α i jk can be expanded as functions of the atomic displacements. At the lowest order expansion, cos α i jk − cos α 0 is related to the out-of-plane displacement and r i j − r 0 is related to the in-plane displacement. For example where i u α denotes the αth component of the displacement of atom i (the related atomic number is labeled in Figure 1). Finally the atomic potential related to motion of atom 1 in Figure 1 takes the following form 1 = V sr (r 12 ) + V sr (r 13 ) + V sr (r 14 ) + V lr (r 15 ) + V lr (r 16 ) + V lr (r 17 ) + V lr (r 18 ) + V lr (r 19 ) + V lr (r 110 ) + V lr (r 111 ) + V lr (r 12 ) + V lr (r 113 ) where and V sr and V lr denote respectively the short range and long range potential. This expression takes into account up to third nearest neighbor interaction. Replacing the force constants in Eq. (11) with the averaged ones in Eq. (3), we applied the SCH, SCE and SCEC models to calculate the vibrational stability of graphene. The instability temperature T i predicted by the three models and the Lindemann ratio 19 at T i is listed in Table I. For isotropic solids, Moleko and Glyde 13 state without strict deduction that instability occur when certain phonon frequency calculated as a function of the lindermann ratio reaches its maximum. Here we show that this conclusion is just an approximation. Consider the SCE model for isotropic solids, Eqs. (1)-(3) is simplified as 13 Eq. (4) is simplified as Eqs. (12) and (13) are calculated iteratively until a pair of self-consistent Einstein frequency and Lindermann ratio is obtained. Assume a solution point (ω 0 , δ 0 ) satisfying Eqs. (12)-(13), then consider an arbitrarily small disturbance of the Lindermann ratio δ p (δ p > 0), the corresponding change of Einstein frequency ω p from Eq. (12) is where higher order terms of δ p and ω p are omitted. Substitution of Eq. (14) into (13) returns the new disturbance of the Lindermann ratio δ p1 after one cycle of iterative calculation: If δ p1 > δ p for any small δ p , solution point (ω 0 , δ 0 ) is completely unstable. This requires (16) Notice that the right hand side of Eq. (16) is always negative, thus strictly speaking the instability point is never at the maximum of f (δ) where d f (δ) dδ = 0. When δ p1 < δ p , solution point (ω 0 , δ 0 ) is conditionally stable which means that there exists a local domain near (ω 0 , δ 0 ) inside which self-consistent solution can be obtained. Generally speaking, the critical point where self-consistent solution becomes completely unstable is quite near the maximum, at which self-consistent calculation already converges very slowly. Thus it is still reasonable and more convenient to consider the maximum of f (δ) as the critical point where instability occurs.
SCP calculation for isotropic solids supports the Lindermann ratio as a criteria concerning the vibrational stability. Yet for graphene things are more complicated. Due to its 2-D atomic structure graphene is strongly anisotropic in that the amplitude of out-of-plane atomic motion is much larger than the amplitude of in-plane atomic motion. To describe this character, two Lindemann ratios are defined for graphene, one related to the in-plane atomic motion by δ 1 = u 2 1 /r 0 and the other related to the out-of-plane atomic motion by δ 3 = u 2 3 /r 0 . For classic materials, T i predicted by harmonic models (SCH, SCE) is often 20∼200 times higher than the melting temperature T m obtained in experiment. For graphene, MC simulation shows a melting temperature of 4900K, 20 which is approximately half of the temperature T i predicted by the three models here. This result suggests that the vibration-initiated instability is vital for graphene concerning its melting. This phenomenon cannot be explained by the Lindemann ratios which, according to the results in Table I, are highly model dependent. Define the vibrational anisotropy coefficient (VAC) for graphene as follow: For any low-D materials, the VAC can be defined in a similar manner. From Eq. (17) we learn that the VAC reflects the anisotropy of the material due to the atomic structure. The change of (δ 1 , δ 3 ) and VAC values during the calculation at different temperatures for three SCP models is illustrated in Figure 2. We found that the VAC values stay in a narrow range for all the temperature points calculated. And although defined from the model-dependent Lindermann ratios, the VAC relies weakly on the models. This means that temperature rising and the chosen model have a slight impact on the VAC once the structure is maintained. To further explore the relation between structural anisotropy and low T i for graphene, we assume the VAC to be an intrinsic parameter of a structure which does not change with temperature. In this case, the phonon frequencies are determined as functions of δ 1 (notice that δ 3 = χδ 1 ) from Eq. further simplified for graphene as: where ω p = ω 1 = ω 2 denotes the phonon frequency for in-plane motion. For graphene at constantχ , we plot the variation of ω 2 p /ω 2 p0 and ω 2 3 /ω 2 30 with δ 1 at different VAC values in Figure 3, where ω p0 and ω 30 denote respectively the phonon frequency for in-plane vibration and out-of-plane vibration at 0K. From the results in Table I, vibrational instability occurs at δ 1 = 0.149 using the SCE model (marked by dotted line in Figure 3(b)). This corresponds to a local maximum of ω 2 3 /ω 2 30 (with the VAC value equals to 3.91 calculated from Table I), but not to the global maximum at aroundδ 1 = 0.5. From Figure 3(b) we can also see that the local maximum of ω 2 3 /ω 2 30 depends highly on the VAC value, that it appears as the VAC value gradually increases. Therefore we come to an important conclusion that the vibrational instability for graphene is due to lost of self-consistent solution of the phonon frequency of the out-of-plane vibration, and the critical temperature and critical Lindermann ratios are closely related to the VAC value of the structure. In Figure 4, it is shown that at different VAC values, the presence, shape and position of the local maximum change dramatically. Further calculation indicates that when χ ≥ 3.8 the local maximum appears and when χ ≥ 5.4 the value of δ 1 corresponding to the local maximum becomes smaller than its value for zero-point vibration (marked by dotted line in Figure 4). This means that for such a graphene-like structure made up of carbon atoms, its VAC value cannot be larger than 5.4.
By calculating the minimum of the Helmholtz free energy of graphene, the EBL at different temperatures predicted by the three models is plotted in Figure 5. It is shown that the predicted EBL is highly dependent on the model. The well-known negative thermal expansion coefficient of graphene at low temperature range is predicted clearly by using the SCH model, yet not by the FIG. 6. Self-consistent in-plane force constant (a) C 1 , (b) C α , predicted by SCH, SCE and SCEC models for graphene. SCE and SCEC models. In the SCE and SCEC models, the phonon spectrum is simplified by two Einstein frequencies in Eq. (18). This simplification underestimates the vibrational amplitude of the atoms (proved by the large discrepancy of the Lindermann ratios listed in Table I obtained from the SCH model and from the SCE, SECS models), and also underestimates the VAC value of the structure (see Figure 2(b)). The physical reason for the second underestimation is that the phonon branch for the out-of-plane vibration is low-lying 26 and has a different dispersion law from the one for in-plane vibration. 15 Thus the high VAC value of graphene is also responsible for the negative thermal expansion coefficient of graphene at low temperature range. Notice that the negative thermal expansion coefficient changes sign at T ≈ 900 K according to Ref. 9 while our SCH model predicts a critical point of T ≈ 5000 K . This difference is contributed to the neglect of nonlinear effects in the SCH model, which tends to accelerate the thermal expansion.
The self-consistent force constants C 1 and C a predicted by SCH, SCE and SCEC models are plotted in Figure 6. When temperature rises, the self-consistent force constants are affected by two factors, lattice expansion and increase of vibrational amplitude, in a rather complicated manner. Generally speaking, the self-consistent force constants decrease with lattice expansion. For graphene, as the vibrational amplitude increases, C 1 first decreases and then increases, while C a first increases, then decreases, and finally increases. From Figure 6 we can see that the self-consistent force constants are dominantly affected by the change of lattice constants.
In this work, we extended three SCP models to deal with the vibrational stability of graphene as temperature rises. A new concept, the vibrational anisotropic coefficient (VAC) is introduced which characterized the peculiar vibrational properties of graphene, as well as its negative thermal expansion coefficient at low temperature range. Moreover, provided the atomic interaction and structure, analyzing the VAC value of the structure gives a criterion regarding the stability of low-D materials. The SCP models used here did not take into account the effects of any possible topologic defects, which tend to decrease the stability of the material. Thus the results obtained from the VAC analysis may be considered as an upper bound for the stability of any low-D materials.