- Email: [email protected]
Out-of-plane shear and out-of plane Young’s modulus of double-layer graphene
Chemical Physics Letters 564 (2013) 37–40
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Out-of-plane shear and out-of plane Young’s modulus of double-layer graphene Balázs Hajgató a,⇑, Songül Güryel a, Yves Dauphin b, Jean-Marie Blairon b, Hans E. Miltner b, Gregory Van Lier a, Frank De Proft a, Paul Geerlings a a b
Free University of Brussels – Vrije Universiteit Brussel (VUB), Research Group General Chemistry (ALGC), Pleinlaan 2, B-1050 Brussels, Belgium SOLVAY S.A., Innovation Center, Ransbeekstraat 310, B-1120 Brussels, Belgium
a r t i c l e
i n f o
Article history: Received 27 December 2012 In final form 12 February 2013 Available online 19 February 2013
a b s t r a c t In this Letter we predict the out-of-plane Young’s modulus (perpendicular to the basal plane) and the out-of-plane shear modulus of double-layer graphene using density functional theory calculations with periodic boundary conditions using the GAUSSIAN 09 program package. These values are discussed in the context of the corresponding values of graphite. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Until the 1985, only two ordered forms of carbon were known to scientists, namely diamond, with its perfect crystal structure, and graphite, also crystalline but black and flaky and not at all transparent. Besides those ordered forms, also coal, coke, soot, lampblack, and the many kinds of charcoal were known. The graphite structure determines its properties, since it is made up of sheets of carbon atoms arranged in a hexagonal lattice, like a honeycomb of fused benzene rings, and with weak bonding between adjacent sheets. This means that graphite easily forms flakes where the sheets can slide over each other, providing use of graphite as a lubricant. In 2007, researchers in Manchester found a way to mechanically peel single two-dimensional sheets from threedimensional graphite crystals [1]. Graphene is the name given to this flat monolayer of carbon atoms tightly packed into a twodimensional honeycomb lattice. Since the first experimental analysis [1], graphene has recently gained significant attention. In particular, its excellent mechanical properties are an important advantage for the practical applications of graphene [2]. These mechanical properties have been extensively investigated, and in particular, the most important elastic property, namely the (in-plane) Young’s modulus has been studied using a wide range of experimental and theoretical approaches (see references in [3]). Recently we also investigated the in-plane Young’s modulus and flexural moduli of double-layer graphene, including a short review about the currently available experimental and theoretical predictions of the (in-plane) Young’s modulus and flexural modulus [3,4]. If one runs through the available experimental an theoretical results of in-plane Young’s modulus, it is apparent that the results have a rather big deviation, ranging from 0.762 to 5.189 TPa (See references in [3]). The origin ⇑ Corresponding author. E-mail address: [email protected] (B. Hajgató). 0009-2614/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2013.02.018
of these big differences lies in the different interpretation of the thickness of atomic platelets, which is not well defined and greatly influences the results. The in-plane shear modulus of single and multiple-layer graphenes have been measured to be 280 and 53 GPa, respectively [5–7]. The out-of-plane shear modulus has been theoretically determined to be in the range of 0.482– 16.1 GPa using molecular structural mechanics and a mixed atomistic continuum finite element technique, respectively [8,9]. There are some experimental studies investigating the transverse elastic properties of suspended graphene sheets, but only the in-plane Young’s modulus was extracted from the measurements [10]. According to our knowledge, up to now, no theoretical or experimental predictions exist for the out-of-plane Young’s modulus of double layer graphene. The classical engineering definitions of these quantities are not unequivocally transferable to atomic scale materials, especially in the case of single-layer graphene, but it can be done in the case of multiple-layer graphene with no more complications than in the case of in-plane Young’s modulus. The available theoretical results for the out-of-plane Young’s modulus of graphite are in the range of 0.3–45.4 GPa [11–18], and the outof-plane shear modulus of graphite is in the range of 3.9–5.0 GPa [13,16,18]. In the case of the out-of-plane Young’s modulus the differences between the obtained results are rather big. If one analyses more deeply those results it turns out that the differences have a methodological origin, since many theoretical methods are not capable to describe correctly the interlayer interactions, and not a consequence of some fundamental difference as in the case of thickness definition for the in-plane Young’s modulus. For multi-layered structures it is a rather good assumption that the weak interlayer interactions do not influence the in-plane Young’s modulus. This assumption is not so evident in the cases of out-ofplane Young’s and shear moduli, since the reacting forces are significantly lower, and most probably a weak interaction can be of influence much more than in the case of the in-plane Young’s modulus.
38
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
The goal of this Letter is to predict the out-of-plane Young’s modulus and shear modulus of double-layer graphene, with special attention to the comparison with the corresponding values of graphite. The in-plane Young’s modulus of double-layer graphene can be approximated through the in-plane Young’s modulus of few-layer graphene and graphite; we investigate whether this approximation is especially valid in the case of the out-of-plane Young’s modulus and the out-of-plane shear modulus. 2. Computational methods and theory All quantum-chemical calculations have been performed using Density Functional Theory (DFT) using the pure LDA [19–23], and PBE [24,25], the short-range hybrid HSE06H [26–32] developed for solids (sometimes referred as HSEh1PBE), and the local M06L [33] functionals, in conjunction with 6-31G⁄ [34] split-valence basis sets using Periodic Boundary Conditions (PBC) implemented in the GAUSSIAN 09 [35] program package. We used the smallest possible unit cells (2 carbons in graphite and 4 carbons in double-layer graphene) in our PBC calculations. Based on our experience, the mechanical properties of single- and double-layer graphene can be calculated accurately enough using the minimum unit cells [3]. To minimize errors arising from numerical instabilities, the convergence in the self consistent field procedure was set to 1010 a.u. on the energy and 108 a.u. on the electron density, ultrafine integration grid (99 radial shells and 590 angular points per radial shells), and tight optimization criteria (15 la.u. on maximum force, 10 la.u. on RMS force, 60 la.u. on predicted maximum displacement, and 40 la.u. on RMS predicted displacement) were used. The Young’s modulus (E) was calculated according to [36]
E¼
tensile stress r AF FL0 ¼ ¼ ¼ ; tensile strain e DL L ADL
ð1Þ
0
where F is the force, A is the area where the force is applied, L0 is the relaxed length, and DL is the elongation/compression. The shear modulus (G) was calculated according to
G¼
tensile stress s AF Fl0 ¼ ¼ ¼ ; tensile strain c Dl z ADz
ð2Þ
0
where F is the force, A is the area where the force is applied, l0 is the relaxed interlayer distance, and Dz is the transverse displacement. (See Figure 1) We note that the methodology applied here is different from the one in Ref. [18], since we directly calculate forces from distorted structures, as we did in our previous studies [3,4,36]. We made use of the engineering approximation in the case of our Young’s and shear moduli calculations, so the geometrical changes arising from the elongation/compression or distortion have been taken into account (in other words, the area in the Eqs. (1) and (2) always corresponds to the actual elongation/compression or distortion geometry, respectively, not the zero strain situation). It is also important to stress that in the case of the in-plane Young’s modulus the thickness of the graphene sheet or the interlayer
Δz
A
F
Figure 1. Definition of the quantities of relevance for the determination of the shear modulus.
distance of graphite can greatly change the predicted value, meanwhile in the case of the out-of-plane Young’s modulus the thickness of a (double or more) layer graphene or the interlayer distance of graphite has almost no influence. This difference is arising from the fact that the area where the force is acting in the case of the out-of-plane Young’s modulus can always be pragmatically defined and can be determined with a relatively high precision (only the inplane atomic distances are needed), while this is not the case for the in-plane Young modulus, because the area cannot always be defined pragmatically and the uncertainty of the determination is also higher (one needs the out-of-plane atomic distances, and/or one should define a platelet thickness of an atomic scale material). We have to emphasize that for the zero-strain point the equilibrium interlayer distance should always be used during the calculation of the out-of-plane Young’s modulus from quantummechanical equations; otherwise the corresponding stress would not be zero. In all our calculations, the area in Eqs. (1) and (2) was calculated purely from the theoretical (optimized) geometries. In addition, in the case of the in-plane Young’s modulus of graphite we also used the experimental interlayer distance (3.35 Å) [37] for comparative purposes. Finally we have to mention that our methodology, due to its nature, is not suitable to determine the out-ofplane Young’s and shear moduli of single-layer graphene. During the out-of-plane Young’s modulus calculations we defined two hypothetical planes (for the two layers of the doublelayer graphene), and the carbon atoms were forced to remain in these planes. The out-of-plane Young’s modulus was calculated using the average force acting on the atoms at 1.25%, 2.50%, and 3.75% elongation and 1.25% compression with respect to the equilibrium interlayer distance. The only further constraint was the increased or decreased distance of the two hypothetical layers, and the atoms were able to freely move within the two defined planes during the geometry optimizations of the elongated/compressed structures. In the case of the out-of-plane shear modulus calculations, two hypothetical planes were defined representing the two layers of the double-layer graphene and one atom of each layer was forced to be in the corresponding plane. Geometry optimizations were performed with the only further constraint being the distance between one of these clamped atoms and the projection of the other clamped atom on the same hypothetical plane (the transverse displacement). Therefore, the transverse movement corresponds to the minimum energy path. The shear modulus was calculated using the average force acting on atoms in 0.07 Å steps in the transverse displacement. (In total 19 steps in between AA and AB stacking, using the minimum energy path). In the case of graphite the Young’s modulus was calculated at 0.5% compression and elongation with respect to the corresponding crystallographic parameter; the atoms in the unit cell and the other crystallographic parameters (lengths and angles) were fully optimized. The corresponding forces were calculated on the elongated or compressed crystallographic parameter. The shear modulus was calculated similarly, but one of the crystallographic angles was fixed, and all the other parameters, including the atomic positions in the unit cell, were optimized. The corresponding shear stress was 0.2% and 0.6% (0.5% and 1.5% change in the interlayer-shift).
3. Results and discussion In order to validate the quality of different functionals and the default parameters in our PBC calculations, we compare the calculated results of mechanical properties with accurate and widely accepted experimental results [38]. Unfortunately, mechanical properties determined by experiments are incomplete for graphene. As an alternative, we validated our methodology for graphite [38] The results are summarized in Table 1. The in-plane Young’s
39
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
modulus is described almost perfect with the LDA and M06L functionals, and calculations with the two other functionals considered yields less accurate results. The origin of this error is that the interlayer distance is overestimated by PBE and HSE06H functionals. The in-plane forces are still well reproduced, since we can see that the in-plane Young’s modulus is perfectly described with all the functionals under consideration if the experimental interlayer distance is used. If we run through the data for the out-of-plane Young’s modulus we see that none of the functionals describe well the forces in the out-of-plane direction. The PBE and HSE06H functionals seem to describe incorrectly the dispersion forces, therefore, it is not surprising that the interlayer distance is too large, and, consequently, the out-of-plane Young’s modulus is underestimated. The M06L functional overestimates significantly the out-of-plane Young’s modulus, but the predicted interlayer distance is almost perfect. The LDA functional is also overestimating the out-of-plane Young’s modulus, with the predicted interlayer distance being too short. In other words, the PBE and HSE06H functionals give too shallow potential energy surfaces with the minimum being too far, and the M06L functional describes the minimum correctly on the too deep potential energy surface. In the case of the LDA functional, the depth of the potential energy surface seems to be accurate, but the interlayer-distances are underestimated. Generally speaking, all the four functionals have some difficulties to describe the dispersion forces perfectly in this case. The shear modulus of graphite is surprisingly well described by the LDA, PBE, and HSE06H functionals, despite the non-correct description of dispersion forces. On the other hand, the M06L functional significantly overestimate the shear forces. Based on the comparison on interlayer distance in graphite and double-layer graphene, we can conclude that the use of graphite experimental interlayer distance in the case of double-layer graphene is a reasonable choice, since the calculated interlayer distances are practically matching. The Young’s modulus of double-layer graphene has been calculated following the same methodology described in our previous paper [3]. The only difference with this contribution is that the elongation and compression were now made perpendicular to the basal plane (out-of-plane, z-direction) and not in plane (xdirection). The calculations might however have more uncertainty than in the case of in-plane Young’s modulus, because the calculated forces are significantly lower. Here, a four-point derivation method was used. The numerical differentiation was performed at 1.25%, 1.25%, 2.50%, and 3.75% strain. The results are summarized in Table 2 for the different DFT functionals. It has to be mentioned that the PBE functional gave too low forces, due to the lack of dispersion forces. If we analyze the results of Table 2 it is apparent that the ratio between out-of-plane Young’s moduli of graphite and double-layer graphene for a given method is dependent on the interlayer distance. The smaller the interlayer distance, the larger the Young’s modulus of double-layer graphene compared to the graphite. However, further analyzing the raw results for the Young’s modulus (Young’s modulus is 61.4 and 26.5 Gpa for compression and elongation, respectively, using a 2 point derivative
Table 2 Out-of-plane Young’s modulus (Ez) of double-layer graphene calculated under periodic boundary conditions. Method ⁄
LDA/6-31G PBE/6-31G⁄ HSE06H/6-31G⁄ M06L/6-31G⁄ Suggested value EGraphene ðExperimentalÞ z , EGraphene ðTheoreticalÞ z
Ez (GPa)
Ratioa
E0 z (GPa)b
50.0 >1c 14.1 70.9
0.830 1.63 1.56 0.441
41.5 – 22.0 31.3 25
a
Ratio ¼
b
E0 z = Ez * Ratio. Calculated forces are too low.
c
see Table 1.
method), it is apparent that the LDA functional gives a rather asymmetric potential energy surface with a too repulsive potential for compression. The latter effect is probably due to the too short interlayer distance. Based on the validation calculations in the case of graphite, we can estimate that the out-of-plane Young’s modulus of double layer graphene is about 2/3 of the out-of-plane Young’s modulus of graphite. The shear modulus of double-layer graphene has been calculated using Eq. (2) and with the Dz value equal to about 0.07 Å. The results for the shear modulus of double-layer graphene are summarized in Table 3, where two different values can be found. One is the shear modulus for AB stacking double-layer graphene (where half of the carbon atoms in the bottom layer have a middle of the ring on the top, the other half of the carbon in atoms in the bottom layer have another carbon atom on the top), the other one corresponds to the AA stacked double-layer graphene (each carbon atom in the bottom layer has another carbon atom on the top). The AB stacked double-layer graphene is more stable, and the corresponding shear is also higher in this alignment. Again, based on the validation calculations in the case of graphite (see Table 1), we can summarize that the (out-of-plane) shear modulus of double-layer graphene is about 1/2 of the (out-of-plane) shear modulus of graphite (4 GPa) [38]. This latter effect can be the consequence of the lack of support from layers below and under the double-layer graphene. The shear strain–stress curves are also calculated. In total 19 Dz values were calculated between the AB and AA stacking of doublelayer graphene. This gives the step size of about 0.07 Å for the Dz, depending on the relaxed geometry determined using different levels of theory. The shear strain–stress curves can be found in Figure 2. Based on the validation calculations in the case of graphite, it seems that the PBE and HSE06H functionals slightly underestimate the shear modulus and, therefore, the calculated stresses of graphitic systems and that M06L functional significantly overestimates the shear modulus and, therefore, the calculated stresses (by about factor of 3) of graphitic systems. It also apparent that the LDA functional gives possibly too large shear stresses due to the underestimated interlayer distances. In the case of the PBE strain–stress curve, the stress goes into a negative value indicating the reversal of the shear stress, and indeed, a very shallow local
Table 1 Different mechanical properties of graphite, compared to the experimental values. E is the in-plane Young’s modulus, Ez is the out-of-plane Young’s modulus, G is the out-of-plane shear modulus.
a
E (TPa) Ez (GPa) G (GPa) Graphite interlayer (Å) Double-layer graphene interlayer (Å) a
LDA/6-31G⁄
PBE/6-31G⁄
HSE06H/6-31G⁄
M06L/6-31G⁄
Experimental [38]
1.083(1.050) 44.0 4.84 3.24 3.26
0.8697(1.013) 22.4 5.31 3.91 3.92
0.9468(1.096) 23.4 4.26 3.88 3.89
1.062(1.087) 82.7 15.4 3.43 3.42
1.06 36.5 4 3.35 –
Values in parenthesis are calculated with the experimental interlayer distance (3.35 Å) [37].
40
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
Table 3 Out-of-plane shear modulus (G) of double-layer graphene calculated under periodic boundary conditions. Method ⁄
LDA/6-31G PBE/6-31G⁄ HSE06H/6-31G⁄ M06L/6-31G⁄ Suggested Value Graphene ðAB;ExperimentalÞ , GGraphene ðAB;TheoreticalÞ
a
Ratio ¼ G
b
G’ = G * Ratio.
G(AB) (GPa)
G(AA) (GPa)
Ratioa
G0 (AB) (GPa)b
G0 (AA) (GPa)b
2.63 1.41 1.93 5.87
2.03 0.229 0.268 2.36
0.83 0.75 0.94 0.26
2.2 1.1 1.8 1.5 2
1.7 0.17 0.25 0.61 0.7
see Table 1.
distances between graphite and double-layer graphene are practically matching, independently from the used level of theory. Acknowledgements The authors acknowledge the financial support of SOLVAY, S. A. The authors also wish to acknowledge support from the COST Materials, Physical and Nanosciences (MPNS) Action MP0901: ‘Designing Novel Materials for Nanodevices – from Theory to Practice (NanoTP)’. References
Figure 2. Shear strain–stress curve for double-layer graphene, using M06L/6-31G⁄, PBE/6-31G⁄, HSE06H/6-31G⁄, and LDA/6-31G⁄ levels of theory.
minimum can be found on the potential-energy surface around 9– 10% shear strain. This local minimum cannot be found on the potential-energy surface in the case of HSE06H and M06L functional, but a definitive deviation can be seen in the stress–strain curves in this region (around 5–10% strain) without a direction change in the shear stress. The origin of this deviation is not clear. However, the stress–strain curves show this effect whatever for each level of theory considered here. The interlayer distance increases about 6–7% when moving from AB to AA stacking; it changes from 3.26 to 3.54 Å, 3.42 to 3.64 Å, 3.89 to 4.12 Å, and 3.92 to 4.15 Å in the cases of LDA, M06L, HSE06H, and PBE functionals, respectively. The middle part (15–25% strain) of the LDA and M06L strain–stress curves are clearly not valid (marked with lighter colors), since these functionals cannot describe bond reordering, which seems to be happening in this regime. Signs of the latter are not present (or at least not in an apparent way) in the case of the other two functionals, again probably due to the inaccurate description of the weak bonding between the two layers (cf. layer distance).
4. Conclusion We have investigated the out-of-plane Young’s and the shear moduli of double-layer graphene. The out-of-plane Young’s modulus (E) of double-layer graphene is approximated to be 25 GPa, which is around 2/3 of graphite’s E value. The estimated value for the out-of-plane shear modulus (G) of double-layer graphene in AB stacking is 2 GPa, which is about half of graphite’s G value. The shear modulus is 0.7 GPa for the double-layer in AA stacking configuration. This latter value is less accurate than our other results, since the different functionals are in less accordance with each other. The use of the graphite experimental interlayer distance, in order to define plate thicknesses of single and few-layer graphenes, is an excellent choice, since the calculated interlayer
[1] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [2] R. Van Noorden, Nature 442 (2006) 228. [3] B. Hajgató, S. Güryel, Y. Dauphin, J.-M. Blairon, H.E. Miltner, G. Van Lier, F. De Proft, P. Geerlings, J. Phys. Chem. C 116 (2012) 22608. [4] S. Güryel, B. Hajgató, Y. Dauphin, J.M. Blairon, H.E. Miltner, F. De Proft, P. Geerlings, G. Van Lier, Phys. Chem. Chem. Phys. 15 (2013) 659. [5] X. Liu, T.H. Metcalf, J.T. Robinson, B.H. Houston, F. Scarpa, Nano Lett. 12 (2012) 1013. [6] X. Liu, T.H. Metcalf, J.T. Robinson, F.K. Perkins, B.H. Houston, Solid State Phenom. 184 (2012) 319. [7] X. Liu, J.T. Robinson, Z. Wei, P.E. Sheehan, B.H. Houston, E.S. Snow, Diamond Relat. Mater. 19 (2010) 875. [8] F. Scarpa, S. Adhikari, R. Chowdhury, Phys. Lett. A 374 (2010) 2053. [9] S.A. Hosseini Kordkheili, H. Moshrefzadeh-Sani, Comput. Mater. Sci. 69 (2013) 335. [10] I.W. Frank, D.M. Tanenbaum, A.M. van der Zande, P.L. McEuen, J. Vac. Sci. Technol., B 25 (2007) 2558. [11] J.C. Boettger, Phys. Rev. B 55 (1997) 11202. [12] M. Hasegawa, K. Nishidate, Phys. Rev. B 70 (2004) 205431. [13] N. Mounet, N. Marzari, Phys. Rev. B 71 (2005) 205214. [14] T. Tohei, A. Kuwabara, F. Oba, I. Tanaka, Phys. Rev. B 73 (2006) 064304. [15] M. Hasegawa, K. Nishidate, H. Iyetomi, Phys. Rev. B 76 (2007) 115424. [16] K.H. Michel, B. Verberck, Phys. Status Solidi B 245 (2008) 2177. [17] G.K.H. Madsen, L. Ferrighi, B. Hammer, J. Phys. Chem. Lett. 1 (2010) 515. [18] G. Savini, Y.J. Dappe, S. Öberg, J.C. Charlier, M.I. Katsnelson, A. Fasolino, Carbon 49 (2011) 62. [19] P.A.M. Dirac, Proc. R. Soc. London, Ser. A 123 (1929) 714. [20] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [21] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [22] J.C. Slater, The Self-Consistent Field for Molecules and Solids, McGraw-Hill, New York, 1974. [23] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200. [24] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [25] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 78 (1997) 1396. [26] T.M. Henderson, A.F. Izmaylov, G. Scalmani, G.E. Scuseria, J. Chem. Phys. 131 (2009) 044108. [27] J. Heyd, J.E. Peralta, G.E. Scuseria, R.L. Martin, J. Phys. Chem. 123 (2005) 174101. [28] J. Heyd, G.E. Scuseria, J. Chem. Phys. 121 (2004) 1187. [29] J. Heyd, G.E. Scuseria, J. Chem. Phys. 120 (2004) 7274. [30] J. Heyd, G.E. Scuseria, M. Ernzerhof, J. Chem. Phys. 124 (2006) 219906. [31] A.F. Izmaylov, G.E. Scuseria, M.J. Frisch, J. Chem. Phys. 125 (2006) 104103. [32] A.V. Krukau, O.A. Vydrov, A.F. Izmaylov, G.E. Scuseria, J. Chem. Phys. 125 (2006) 224106. [33] Y. Zhao, D.G. Truhlar, J. Chem. Phys. 125 (2006) 194101. [34] M.M. Francl, W.J. Pietro, W.J. Hehre, J.S. Binkley, M.S. Gordon, D.J. DeFrees, J.A. Pople, J. Chem. Phys. 77 (1982) 3654. [35] M.J. Frisch et al., GAUSSIAN 09, Gaussian, Inc., Wallingford CT, 2009. [36] G. Van Lier, C. Van Alsenoy, V. Van Doren, P. Geerlings, Chem. Phys. Lett. 326 (2000) 181. [37] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [38] O.L. Blakslee, D.G. Proctor, E.J. Seldin, G.B. Spence, T. Weng, J. Appl. Phys. 41 (1970) 3373.
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Out-of-plane shear and out-of plane Young’s modulus of double-layer graphene Balázs Hajgató a,⇑, Songül Güryel a, Yves Dauphin b, Jean-Marie Blairon b, Hans E. Miltner b, Gregory Van Lier a, Frank De Proft a, Paul Geerlings a a b
Free University of Brussels – Vrije Universiteit Brussel (VUB), Research Group General Chemistry (ALGC), Pleinlaan 2, B-1050 Brussels, Belgium SOLVAY S.A., Innovation Center, Ransbeekstraat 310, B-1120 Brussels, Belgium
a r t i c l e
i n f o
Article history: Received 27 December 2012 In final form 12 February 2013 Available online 19 February 2013
a b s t r a c t In this Letter we predict the out-of-plane Young’s modulus (perpendicular to the basal plane) and the out-of-plane shear modulus of double-layer graphene using density functional theory calculations with periodic boundary conditions using the GAUSSIAN 09 program package. These values are discussed in the context of the corresponding values of graphite. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction Until the 1985, only two ordered forms of carbon were known to scientists, namely diamond, with its perfect crystal structure, and graphite, also crystalline but black and flaky and not at all transparent. Besides those ordered forms, also coal, coke, soot, lampblack, and the many kinds of charcoal were known. The graphite structure determines its properties, since it is made up of sheets of carbon atoms arranged in a hexagonal lattice, like a honeycomb of fused benzene rings, and with weak bonding between adjacent sheets. This means that graphite easily forms flakes where the sheets can slide over each other, providing use of graphite as a lubricant. In 2007, researchers in Manchester found a way to mechanically peel single two-dimensional sheets from threedimensional graphite crystals [1]. Graphene is the name given to this flat monolayer of carbon atoms tightly packed into a twodimensional honeycomb lattice. Since the first experimental analysis [1], graphene has recently gained significant attention. In particular, its excellent mechanical properties are an important advantage for the practical applications of graphene [2]. These mechanical properties have been extensively investigated, and in particular, the most important elastic property, namely the (in-plane) Young’s modulus has been studied using a wide range of experimental and theoretical approaches (see references in [3]). Recently we also investigated the in-plane Young’s modulus and flexural moduli of double-layer graphene, including a short review about the currently available experimental and theoretical predictions of the (in-plane) Young’s modulus and flexural modulus [3,4]. If one runs through the available experimental an theoretical results of in-plane Young’s modulus, it is apparent that the results have a rather big deviation, ranging from 0.762 to 5.189 TPa (See references in [3]). The origin ⇑ Corresponding author. E-mail address: [email protected] (B. Hajgató). 0009-2614/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2013.02.018
of these big differences lies in the different interpretation of the thickness of atomic platelets, which is not well defined and greatly influences the results. The in-plane shear modulus of single and multiple-layer graphenes have been measured to be 280 and 53 GPa, respectively [5–7]. The out-of-plane shear modulus has been theoretically determined to be in the range of 0.482– 16.1 GPa using molecular structural mechanics and a mixed atomistic continuum finite element technique, respectively [8,9]. There are some experimental studies investigating the transverse elastic properties of suspended graphene sheets, but only the in-plane Young’s modulus was extracted from the measurements [10]. According to our knowledge, up to now, no theoretical or experimental predictions exist for the out-of-plane Young’s modulus of double layer graphene. The classical engineering definitions of these quantities are not unequivocally transferable to atomic scale materials, especially in the case of single-layer graphene, but it can be done in the case of multiple-layer graphene with no more complications than in the case of in-plane Young’s modulus. The available theoretical results for the out-of-plane Young’s modulus of graphite are in the range of 0.3–45.4 GPa [11–18], and the outof-plane shear modulus of graphite is in the range of 3.9–5.0 GPa [13,16,18]. In the case of the out-of-plane Young’s modulus the differences between the obtained results are rather big. If one analyses more deeply those results it turns out that the differences have a methodological origin, since many theoretical methods are not capable to describe correctly the interlayer interactions, and not a consequence of some fundamental difference as in the case of thickness definition for the in-plane Young’s modulus. For multi-layered structures it is a rather good assumption that the weak interlayer interactions do not influence the in-plane Young’s modulus. This assumption is not so evident in the cases of out-ofplane Young’s and shear moduli, since the reacting forces are significantly lower, and most probably a weak interaction can be of influence much more than in the case of the in-plane Young’s modulus.
38
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
The goal of this Letter is to predict the out-of-plane Young’s modulus and shear modulus of double-layer graphene, with special attention to the comparison with the corresponding values of graphite. The in-plane Young’s modulus of double-layer graphene can be approximated through the in-plane Young’s modulus of few-layer graphene and graphite; we investigate whether this approximation is especially valid in the case of the out-of-plane Young’s modulus and the out-of-plane shear modulus. 2. Computational methods and theory All quantum-chemical calculations have been performed using Density Functional Theory (DFT) using the pure LDA [19–23], and PBE [24,25], the short-range hybrid HSE06H [26–32] developed for solids (sometimes referred as HSEh1PBE), and the local M06L [33] functionals, in conjunction with 6-31G⁄ [34] split-valence basis sets using Periodic Boundary Conditions (PBC) implemented in the GAUSSIAN 09 [35] program package. We used the smallest possible unit cells (2 carbons in graphite and 4 carbons in double-layer graphene) in our PBC calculations. Based on our experience, the mechanical properties of single- and double-layer graphene can be calculated accurately enough using the minimum unit cells [3]. To minimize errors arising from numerical instabilities, the convergence in the self consistent field procedure was set to 1010 a.u. on the energy and 108 a.u. on the electron density, ultrafine integration grid (99 radial shells and 590 angular points per radial shells), and tight optimization criteria (15 la.u. on maximum force, 10 la.u. on RMS force, 60 la.u. on predicted maximum displacement, and 40 la.u. on RMS predicted displacement) were used. The Young’s modulus (E) was calculated according to [36]
E¼
tensile stress r AF FL0 ¼ ¼ ¼ ; tensile strain e DL L ADL
ð1Þ
0
where F is the force, A is the area where the force is applied, L0 is the relaxed length, and DL is the elongation/compression. The shear modulus (G) was calculated according to
G¼
tensile stress s AF Fl0 ¼ ¼ ¼ ; tensile strain c Dl z ADz
ð2Þ
0
where F is the force, A is the area where the force is applied, l0 is the relaxed interlayer distance, and Dz is the transverse displacement. (See Figure 1) We note that the methodology applied here is different from the one in Ref. [18], since we directly calculate forces from distorted structures, as we did in our previous studies [3,4,36]. We made use of the engineering approximation in the case of our Young’s and shear moduli calculations, so the geometrical changes arising from the elongation/compression or distortion have been taken into account (in other words, the area in the Eqs. (1) and (2) always corresponds to the actual elongation/compression or distortion geometry, respectively, not the zero strain situation). It is also important to stress that in the case of the in-plane Young’s modulus the thickness of the graphene sheet or the interlayer
Δz
A
F
Figure 1. Definition of the quantities of relevance for the determination of the shear modulus.
distance of graphite can greatly change the predicted value, meanwhile in the case of the out-of-plane Young’s modulus the thickness of a (double or more) layer graphene or the interlayer distance of graphite has almost no influence. This difference is arising from the fact that the area where the force is acting in the case of the out-of-plane Young’s modulus can always be pragmatically defined and can be determined with a relatively high precision (only the inplane atomic distances are needed), while this is not the case for the in-plane Young modulus, because the area cannot always be defined pragmatically and the uncertainty of the determination is also higher (one needs the out-of-plane atomic distances, and/or one should define a platelet thickness of an atomic scale material). We have to emphasize that for the zero-strain point the equilibrium interlayer distance should always be used during the calculation of the out-of-plane Young’s modulus from quantummechanical equations; otherwise the corresponding stress would not be zero. In all our calculations, the area in Eqs. (1) and (2) was calculated purely from the theoretical (optimized) geometries. In addition, in the case of the in-plane Young’s modulus of graphite we also used the experimental interlayer distance (3.35 Å) [37] for comparative purposes. Finally we have to mention that our methodology, due to its nature, is not suitable to determine the out-ofplane Young’s and shear moduli of single-layer graphene. During the out-of-plane Young’s modulus calculations we defined two hypothetical planes (for the two layers of the doublelayer graphene), and the carbon atoms were forced to remain in these planes. The out-of-plane Young’s modulus was calculated using the average force acting on the atoms at 1.25%, 2.50%, and 3.75% elongation and 1.25% compression with respect to the equilibrium interlayer distance. The only further constraint was the increased or decreased distance of the two hypothetical layers, and the atoms were able to freely move within the two defined planes during the geometry optimizations of the elongated/compressed structures. In the case of the out-of-plane shear modulus calculations, two hypothetical planes were defined representing the two layers of the double-layer graphene and one atom of each layer was forced to be in the corresponding plane. Geometry optimizations were performed with the only further constraint being the distance between one of these clamped atoms and the projection of the other clamped atom on the same hypothetical plane (the transverse displacement). Therefore, the transverse movement corresponds to the minimum energy path. The shear modulus was calculated using the average force acting on atoms in 0.07 Å steps in the transverse displacement. (In total 19 steps in between AA and AB stacking, using the minimum energy path). In the case of graphite the Young’s modulus was calculated at 0.5% compression and elongation with respect to the corresponding crystallographic parameter; the atoms in the unit cell and the other crystallographic parameters (lengths and angles) were fully optimized. The corresponding forces were calculated on the elongated or compressed crystallographic parameter. The shear modulus was calculated similarly, but one of the crystallographic angles was fixed, and all the other parameters, including the atomic positions in the unit cell, were optimized. The corresponding shear stress was 0.2% and 0.6% (0.5% and 1.5% change in the interlayer-shift).
3. Results and discussion In order to validate the quality of different functionals and the default parameters in our PBC calculations, we compare the calculated results of mechanical properties with accurate and widely accepted experimental results [38]. Unfortunately, mechanical properties determined by experiments are incomplete for graphene. As an alternative, we validated our methodology for graphite [38] The results are summarized in Table 1. The in-plane Young’s
39
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
modulus is described almost perfect with the LDA and M06L functionals, and calculations with the two other functionals considered yields less accurate results. The origin of this error is that the interlayer distance is overestimated by PBE and HSE06H functionals. The in-plane forces are still well reproduced, since we can see that the in-plane Young’s modulus is perfectly described with all the functionals under consideration if the experimental interlayer distance is used. If we run through the data for the out-of-plane Young’s modulus we see that none of the functionals describe well the forces in the out-of-plane direction. The PBE and HSE06H functionals seem to describe incorrectly the dispersion forces, therefore, it is not surprising that the interlayer distance is too large, and, consequently, the out-of-plane Young’s modulus is underestimated. The M06L functional overestimates significantly the out-of-plane Young’s modulus, but the predicted interlayer distance is almost perfect. The LDA functional is also overestimating the out-of-plane Young’s modulus, with the predicted interlayer distance being too short. In other words, the PBE and HSE06H functionals give too shallow potential energy surfaces with the minimum being too far, and the M06L functional describes the minimum correctly on the too deep potential energy surface. In the case of the LDA functional, the depth of the potential energy surface seems to be accurate, but the interlayer-distances are underestimated. Generally speaking, all the four functionals have some difficulties to describe the dispersion forces perfectly in this case. The shear modulus of graphite is surprisingly well described by the LDA, PBE, and HSE06H functionals, despite the non-correct description of dispersion forces. On the other hand, the M06L functional significantly overestimate the shear forces. Based on the comparison on interlayer distance in graphite and double-layer graphene, we can conclude that the use of graphite experimental interlayer distance in the case of double-layer graphene is a reasonable choice, since the calculated interlayer distances are practically matching. The Young’s modulus of double-layer graphene has been calculated following the same methodology described in our previous paper [3]. The only difference with this contribution is that the elongation and compression were now made perpendicular to the basal plane (out-of-plane, z-direction) and not in plane (xdirection). The calculations might however have more uncertainty than in the case of in-plane Young’s modulus, because the calculated forces are significantly lower. Here, a four-point derivation method was used. The numerical differentiation was performed at 1.25%, 1.25%, 2.50%, and 3.75% strain. The results are summarized in Table 2 for the different DFT functionals. It has to be mentioned that the PBE functional gave too low forces, due to the lack of dispersion forces. If we analyze the results of Table 2 it is apparent that the ratio between out-of-plane Young’s moduli of graphite and double-layer graphene for a given method is dependent on the interlayer distance. The smaller the interlayer distance, the larger the Young’s modulus of double-layer graphene compared to the graphite. However, further analyzing the raw results for the Young’s modulus (Young’s modulus is 61.4 and 26.5 Gpa for compression and elongation, respectively, using a 2 point derivative
Table 2 Out-of-plane Young’s modulus (Ez) of double-layer graphene calculated under periodic boundary conditions. Method ⁄
LDA/6-31G PBE/6-31G⁄ HSE06H/6-31G⁄ M06L/6-31G⁄ Suggested value EGraphene ðExperimentalÞ z , EGraphene ðTheoreticalÞ z
Ez (GPa)
Ratioa
E0 z (GPa)b
50.0 >1c 14.1 70.9
0.830 1.63 1.56 0.441
41.5 – 22.0 31.3 25
a
Ratio ¼
b
E0 z = Ez * Ratio. Calculated forces are too low.
c
see Table 1.
method), it is apparent that the LDA functional gives a rather asymmetric potential energy surface with a too repulsive potential for compression. The latter effect is probably due to the too short interlayer distance. Based on the validation calculations in the case of graphite, we can estimate that the out-of-plane Young’s modulus of double layer graphene is about 2/3 of the out-of-plane Young’s modulus of graphite. The shear modulus of double-layer graphene has been calculated using Eq. (2) and with the Dz value equal to about 0.07 Å. The results for the shear modulus of double-layer graphene are summarized in Table 3, where two different values can be found. One is the shear modulus for AB stacking double-layer graphene (where half of the carbon atoms in the bottom layer have a middle of the ring on the top, the other half of the carbon in atoms in the bottom layer have another carbon atom on the top), the other one corresponds to the AA stacked double-layer graphene (each carbon atom in the bottom layer has another carbon atom on the top). The AB stacked double-layer graphene is more stable, and the corresponding shear is also higher in this alignment. Again, based on the validation calculations in the case of graphite (see Table 1), we can summarize that the (out-of-plane) shear modulus of double-layer graphene is about 1/2 of the (out-of-plane) shear modulus of graphite (4 GPa) [38]. This latter effect can be the consequence of the lack of support from layers below and under the double-layer graphene. The shear strain–stress curves are also calculated. In total 19 Dz values were calculated between the AB and AA stacking of doublelayer graphene. This gives the step size of about 0.07 Å for the Dz, depending on the relaxed geometry determined using different levels of theory. The shear strain–stress curves can be found in Figure 2. Based on the validation calculations in the case of graphite, it seems that the PBE and HSE06H functionals slightly underestimate the shear modulus and, therefore, the calculated stresses of graphitic systems and that M06L functional significantly overestimates the shear modulus and, therefore, the calculated stresses (by about factor of 3) of graphitic systems. It also apparent that the LDA functional gives possibly too large shear stresses due to the underestimated interlayer distances. In the case of the PBE strain–stress curve, the stress goes into a negative value indicating the reversal of the shear stress, and indeed, a very shallow local
Table 1 Different mechanical properties of graphite, compared to the experimental values. E is the in-plane Young’s modulus, Ez is the out-of-plane Young’s modulus, G is the out-of-plane shear modulus.
a
E (TPa) Ez (GPa) G (GPa) Graphite interlayer (Å) Double-layer graphene interlayer (Å) a
LDA/6-31G⁄
PBE/6-31G⁄
HSE06H/6-31G⁄
M06L/6-31G⁄
Experimental [38]
1.083(1.050) 44.0 4.84 3.24 3.26
0.8697(1.013) 22.4 5.31 3.91 3.92
0.9468(1.096) 23.4 4.26 3.88 3.89
1.062(1.087) 82.7 15.4 3.43 3.42
1.06 36.5 4 3.35 –
Values in parenthesis are calculated with the experimental interlayer distance (3.35 Å) [37].
40
B. Hajgató et al. / Chemical Physics Letters 564 (2013) 37–40
Table 3 Out-of-plane shear modulus (G) of double-layer graphene calculated under periodic boundary conditions. Method ⁄
LDA/6-31G PBE/6-31G⁄ HSE06H/6-31G⁄ M06L/6-31G⁄ Suggested Value Graphene ðAB;ExperimentalÞ , GGraphene ðAB;TheoreticalÞ
a
Ratio ¼ G
b
G’ = G * Ratio.
G(AB) (GPa)
G(AA) (GPa)
Ratioa
G0 (AB) (GPa)b
G0 (AA) (GPa)b
2.63 1.41 1.93 5.87
2.03 0.229 0.268 2.36
0.83 0.75 0.94 0.26
2.2 1.1 1.8 1.5 2
1.7 0.17 0.25 0.61 0.7
see Table 1.
distances between graphite and double-layer graphene are practically matching, independently from the used level of theory. Acknowledgements The authors acknowledge the financial support of SOLVAY, S. A. The authors also wish to acknowledge support from the COST Materials, Physical and Nanosciences (MPNS) Action MP0901: ‘Designing Novel Materials for Nanodevices – from Theory to Practice (NanoTP)’. References
Figure 2. Shear strain–stress curve for double-layer graphene, using M06L/6-31G⁄, PBE/6-31G⁄, HSE06H/6-31G⁄, and LDA/6-31G⁄ levels of theory.
minimum can be found on the potential-energy surface around 9– 10% shear strain. This local minimum cannot be found on the potential-energy surface in the case of HSE06H and M06L functional, but a definitive deviation can be seen in the stress–strain curves in this region (around 5–10% strain) without a direction change in the shear stress. The origin of this deviation is not clear. However, the stress–strain curves show this effect whatever for each level of theory considered here. The interlayer distance increases about 6–7% when moving from AB to AA stacking; it changes from 3.26 to 3.54 Å, 3.42 to 3.64 Å, 3.89 to 4.12 Å, and 3.92 to 4.15 Å in the cases of LDA, M06L, HSE06H, and PBE functionals, respectively. The middle part (15–25% strain) of the LDA and M06L strain–stress curves are clearly not valid (marked with lighter colors), since these functionals cannot describe bond reordering, which seems to be happening in this regime. Signs of the latter are not present (or at least not in an apparent way) in the case of the other two functionals, again probably due to the inaccurate description of the weak bonding between the two layers (cf. layer distance).
4. Conclusion We have investigated the out-of-plane Young’s and the shear moduli of double-layer graphene. The out-of-plane Young’s modulus (E) of double-layer graphene is approximated to be 25 GPa, which is around 2/3 of graphite’s E value. The estimated value for the out-of-plane shear modulus (G) of double-layer graphene in AB stacking is 2 GPa, which is about half of graphite’s G value. The shear modulus is 0.7 GPa for the double-layer in AA stacking configuration. This latter value is less accurate than our other results, since the different functionals are in less accordance with each other. The use of the graphite experimental interlayer distance, in order to define plate thicknesses of single and few-layer graphenes, is an excellent choice, since the calculated interlayer
[1] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [2] R. Van Noorden, Nature 442 (2006) 228. [3] B. Hajgató, S. Güryel, Y. Dauphin, J.-M. Blairon, H.E. Miltner, G. Van Lier, F. De Proft, P. Geerlings, J. Phys. Chem. C 116 (2012) 22608. [4] S. Güryel, B. Hajgató, Y. Dauphin, J.M. Blairon, H.E. Miltner, F. De Proft, P. Geerlings, G. Van Lier, Phys. Chem. Chem. Phys. 15 (2013) 659. [5] X. Liu, T.H. Metcalf, J.T. Robinson, B.H. Houston, F. Scarpa, Nano Lett. 12 (2012) 1013. [6] X. Liu, T.H. Metcalf, J.T. Robinson, F.K. Perkins, B.H. Houston, Solid State Phenom. 184 (2012) 319. [7] X. Liu, J.T. Robinson, Z. Wei, P.E. Sheehan, B.H. Houston, E.S. Snow, Diamond Relat. Mater. 19 (2010) 875. [8] F. Scarpa, S. Adhikari, R. Chowdhury, Phys. Lett. A 374 (2010) 2053. [9] S.A. Hosseini Kordkheili, H. Moshrefzadeh-Sani, Comput. Mater. Sci. 69 (2013) 335. [10] I.W. Frank, D.M. Tanenbaum, A.M. van der Zande, P.L. McEuen, J. Vac. Sci. Technol., B 25 (2007) 2558. [11] J.C. Boettger, Phys. Rev. B 55 (1997) 11202. [12] M. Hasegawa, K. Nishidate, Phys. Rev. B 70 (2004) 205431. [13] N. Mounet, N. Marzari, Phys. Rev. B 71 (2005) 205214. [14] T. Tohei, A. Kuwabara, F. Oba, I. Tanaka, Phys. Rev. B 73 (2006) 064304. [15] M. Hasegawa, K. Nishidate, H. Iyetomi, Phys. Rev. B 76 (2007) 115424. [16] K.H. Michel, B. Verberck, Phys. Status Solidi B 245 (2008) 2177. [17] G.K.H. Madsen, L. Ferrighi, B. Hammer, J. Phys. Chem. Lett. 1 (2010) 515. [18] G. Savini, Y.J. Dappe, S. Öberg, J.C. Charlier, M.I. Katsnelson, A. Fasolino, Carbon 49 (2011) 62. [19] P.A.M. Dirac, Proc. R. Soc. London, Ser. A 123 (1929) 714. [20] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [21] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133. [22] J.C. Slater, The Self-Consistent Field for Molecules and Solids, McGraw-Hill, New York, 1974. [23] S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200. [24] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [25] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 78 (1997) 1396. [26] T.M. Henderson, A.F. Izmaylov, G. Scalmani, G.E. Scuseria, J. Chem. Phys. 131 (2009) 044108. [27] J. Heyd, J.E. Peralta, G.E. Scuseria, R.L. Martin, J. Phys. Chem. 123 (2005) 174101. [28] J. Heyd, G.E. Scuseria, J. Chem. Phys. 121 (2004) 1187. [29] J. Heyd, G.E. Scuseria, J. Chem. Phys. 120 (2004) 7274. [30] J. Heyd, G.E. Scuseria, M. Ernzerhof, J. Chem. Phys. 124 (2006) 219906. [31] A.F. Izmaylov, G.E. Scuseria, M.J. Frisch, J. Chem. Phys. 125 (2006) 104103. [32] A.V. Krukau, O.A. Vydrov, A.F. Izmaylov, G.E. Scuseria, J. Chem. Phys. 125 (2006) 224106. [33] Y. Zhao, D.G. Truhlar, J. Chem. Phys. 125 (2006) 194101. [34] M.M. Francl, W.J. Pietro, W.J. Hehre, J.S. Binkley, M.S. Gordon, D.J. DeFrees, J.A. Pople, J. Chem. Phys. 77 (1982) 3654. [35] M.J. Frisch et al., GAUSSIAN 09, Gaussian, Inc., Wallingford CT, 2009. [36] G. Van Lier, C. Van Alsenoy, V. Van Doren, P. Geerlings, Chem. Phys. Lett. 326 (2000) 181. [37] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [38] O.L. Blakslee, D.G. Proctor, E.J. Seldin, G.B. Spence, T. Weng, J. Appl. Phys. 41 (1970) 3373.