Graphene: A partially ordered non-periodic solid Dongshan Wei and Feng Wang Citation: The Journal of Chemical Physics 141, 144701 (2014); doi: 10.1063/1.4897255 View online: http://dx.doi.org/10.1063/1.4897255 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/14?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Diffusion, adsorption, and desorption of molecular hydrogen on graphene and in graphite J. Chem. Phys. 139, 044706 (2013); 10.1063/1.4813919 Vibrational stability of graphene AIP Advances 3, 052101 (2013); 10.1063/1.4804244 Reaction between graphene and hydrogen under oblique injection J. Appl. Phys. 110, 084320 (2011); 10.1063/1.3651394 Diffusion and drift of graphene flake on graphite surface J. Chem. Phys. 134, 104505 (2011); 10.1063/1.3557819 Study of the Critical Behavior of Nonequilibrium Systems: Application to Driven Diffusive Lattice Gases AIP Conf. Proc. 661, 90 (2003); 10.1063/1.1571296
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
THE JOURNAL OF CHEMICAL PHYSICS 141, 144701 (2014)
Graphene: A partially ordered non-periodic solid Dongshan Wei1 and Feng Wang2,a) 1
Key Laboratory of Multi-Scale Manufacturing Technology, Chongqing Institute of Green and Intelligent Technology, Chinese Academy of Sciences, Chongqing 400714, China 2 Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, USA
(Received 2 June 2014; accepted 24 September 2014; published online 9 October 2014) Molecular dynamics simulations were performed to study the structural features of graphene over a wide range of temperatures from 50 to 4000 K using the PPBE-G potential [D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134, 184704 (2011)]. This potential was developed by force matching the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional and has been validated previously to provide accurate potential energy surface for graphene at temperatures as high as 3000 K. Simulations with the PPBE–G potential are the best available approximation to a direct Car-Parrinello Molecular Dynamics study of graphene. One advantage of the PBE-G potential is to allow large simulation boxes to be modeled efficiently so that properties showing strong finite size effects can be studied. Our simulation box contains more than 600 000 C atoms and is one of the largest graphene boxes ever modeled. With the PPBE-G potential, the thermal-expansion coefficient is negative up to 4000 K. With a large box and an accurate potential, the critical exponent for the scaling properties associated with the normal-normal and height-height correlation functions was confirmed to be 0.85. This exponent remains constant up to 4000 K suggesting graphene to be in the deeply cooled regime even close to the experimental melting temperature. The reduced peak heights in the radial distribution function of graphene show an inverse power law dependence to distance, which indicates that a macroscopic graphene sheet will lose long-range crystalline order as predicted by the Mermin-Wagner instability. Although graphene loses long-range translational order, it retains long range orientational order as indicated by its orientational correlation function; graphene is thus partially ordered but not periodic. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/1.4897255] I. INTRODUCTION
Graphene is a two-dimensional (2D) material that has the structure of a graphite monolayer.1 Perfectly flat graphene will have a zero band gap. This zero band gap leads to many superior electronic properties, making graphene a promising material for building ground breaking electronics. However, at any finite temperature, a graphene sheet cannot be perfectly flat. Spontaneous rippling of the sheet has been observed experimentally.2 Such rippling will cause the bandgap to deviate from zero3–6 and produce electron-ripple scattering.7, 8 In order to fully understand the potential of graphene as a new material, it is important to investigate the effect of rippling for large graphene sheets. Due to thermo-fluctuations, no 2D material can be perfectly flat at any finite temperature. Ignoring the surface tension term, the free energy of a 2D membrane can be written as 1 dxdy(∇ 2 h)2 F = κ 2 1 (1) + dxdy λ u2xx + u2yy + 2μu2xy , 2
a) [email protected]
0021-9606/2014/141(14)/144701/8
where h is the normal displacement, κ is the bending rigidity, λ and μ are Lamé’s first and second parameters, respectively. In Eq. (1), the strain matrix elements uαβ are defined as ∂uβ ∂h ∂h 1 ∂uα + + uαβ ≈ , (2) 2 ∂β ∂α ∂α ∂β where uα and uβ are in-plane displacements. Equation (1) is applicable to both solid and liquid membranes. The experimental observation of large graphene sheets surprised the scientific community initially. This was due to the argument that no 2D material can exist over a large length scale. It is worth noting that there are two similar arguments regarding the stability of 2D materials, one by Peierls9 and Landau10 and another one by Peliti and Leibler.11 The PeierlsLandau argument states that there could be no long range crystalline order in 1D and 2D materials. This is generally referred to as the Mermin-Wagner instability.12 This statement allows the existence of 1D or 2D materials but requires such materials not to be periodic at macroscopic scale. In other words, such reduced dimensional material can exist but lacks long range translational order. The argument by Peliti and Leibler is stronger. They argues that no 2D material can exisit at macrosopic scale. This statement was made by assuming the second integral in Eq. (1) to vanish, which is true when the coupling between
141, 144701-1
© Author(s) 2014
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-2
D. Wei and F. Wang
J. Chem. Phys. 141, 144701 (2014)
normal and in-plane displacements is negligible. Under this condition, one can show that the renormalized bending rigidity depends on membrane size and to the first order can be described by the equation, κ(L) = κ0 − (3kB T /4π ) ln(L/a),
(3)
where κ 0 is the bare bending rigidity and a is a certain microscopic length. In this scenario, the bending rigidity of the 2D material diverges to negative infinity when its size is asymptotically large; the material spontaneously wraps around to form 3D structures, such as vesicles. The assumption made by Peliti and Leibler holds for liquid membranes. A solid membrane can potentially satisfy Eq. (3) if the coupling between in-plane and out-of-plane displacement is weak. The renormalized bending rigidity of a material that satisfies Eq. (3) will have a negative temperature dependence. Whether graphene satisfies Eq. (3) is a subject of debate. Simulations using the REBO-II potential13 shows negative temperature dependence14 of κ(L) consistent with Eq. (3) while the LCBOPII potential15 produces a positive dependence.16 For 2D materials that do satisfy Eq. (3), the length at which κ(L) vanishes is frequently referred to as the persistence length ξ .11 ξ = a · e4πκ0 /3kB T .
(4)
Such a material will undergo the crumpling phase transition when its size is larger than ξ . Another popular definition for the persistence length is through the normalnormal correlation function G(R) = n(0) · n(R). For weakly coupled membranes, the normal-normal correlation function is expected to decay exponentially with respect to the separation R, n(0) · n(R) = e−R/ξn .
(5)
The decay rate can be fit by a single exponential to obtain the persistence length ξ n of the normal. For a liquid membrane, such as an oil-water interface saturated with surfac2πκ 0
tants, De Gennes and Taupin17 argue that ξn = a · e 3kB T thus 2πκ 0
ξ = ξn · e 3kB T or ξ = ξn2 /a. In the following discussion, we will refer to ξ as the crumpling persistence length and ξ n as the normal persistence length. Whether graphene is susceptible to a crumpling phase transition depends on the model potential being used. For example, the LCBOPII potential16, 18 shows a complete breakdown of Eq. (3), but simulation with the REBO-II potential14 shows a behavior that is compatible to Eq. (3). Although one could argue that the experimental observation of large suspended graphene will rule out the crumpling instability predicted in Eq. (3), this is not true. The bare bending rigidity of graphene has been estimated to be in the range of 1 eV or higher.16, 19, 20 Even if the thermodynamics of graphene were governed by Eq. (3), at 300 K, with an a of 1.4 Å, Eq. (4) predicts a ξ of 3.3 × 1060 m and a ξ n of 2.2 × 1025 m, both are astronomically large. At 3000 K, ξ would still be microscopically large at 1.5 mm and ξ n would be 0.46 μm. We note referring to graphene as a solid membrane or a membrane with local crystalline order is probably correct
but not very precise. A lipid membrane can be in a 2D crystalline state at low temperatures. Although there is local crystalline order in such a state, there are no covalent bonds between lipid molecules. Graphene is more precisely referred to as a tethered surface or polymerized membrane.21 For a tethered membrane, the bonding connectivity between neighboring molecules cannot change. This is true for graphene under a wide range of temperatures. Experimentally, graphite will not melt until around 4800 K.22 The Stone-Wales topological defect is expected to have a formation energy of 4.7 eV.23 Even at 4000 K, the thermo-population of such defects is around 10−6 . The difference between a tethered membrane with fixed connectivity and a free diffusing membrane with weak intermolecular interactions can be illustrated by a vivid example given by Kantor.24 Although a fluid membrane can crumple and form a sphere, such as a vesicle, a tethered membrane, such as a piece of paper, can never form a sphere. The overall “shape” of a tethered surface cannot change. The Lamé’s parameters for a tethered surface can be large since bending to the third dimension cannot be accomplished without causing strain in the in-plane direction by stretching and compressing covalent bonds. The thermodynamics of a tethered surface is far from well-understood. It has been argued that a tethered membrane will crumple when κ is less than a critical value κ c . When κ is sufficiently large, a tethered membrane will remain macroscopically flat. We note having the tethered membrane to remain macroscopically flat does not contradict with the Mermin–Wagner instability.12 As mentioned previously, the Mermin-Wagner instability only rules out long range translational order but allows the existence of 2D materials. Another interesting property that can be derived for nonself-avoiding (phantom) flat membrane is the existence of a universal scaling factor as a function of size for the power spectra of normal-normal correlation function and the heightheight correlation function. The existence of universal scaling constant is frequently observed in 3D materials near a critical point. A universality class behavior in 2D away from a critical point is very intriguing.18 Several theoretical25, 26 and numerical27–29 investigations of the critical scaling exponent, η, has been derived for 2D membrane in 3D space. The most accurate estimate available is probably η = 0.849 from renormalization group theory,25 although other estimates such as 0.821,26 0.750,29 and 0.8127 have also been reported. In this work, we investigate the thermodynamics of a large graphene sheet using the recently developed PPBEG potential.30 The PPBE-G potential was created by forcematching30–33 the Perdew-Burke-Ernzerhof (PBE) forces calculated with the Projector Augmented Wave (PAW) treatment in the temperature range from 300 to 3000 K.34, 35 In the temperature range of parameterization, the PPBE-G potential produces an atomic force root-mean-square-error (RMSE) of 4.4 kcal/(mol Å) when compared to PAW-PBE references.30 Since other graphene force fields surveyed produce a force RMSE of at least 14 kcal/(mol Å) when compared to PAWPBE, a PPBE-G simulation is the best available approximation to a direct Car-Parrinello molecular dynamics36 (CPMD) study for graphene.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-3
D. Wei and F. Wang
One clear advantage of the PPBE-G potential is its high computational efficiency. This efficiency allows graphene sheets with millions of atoms to be studied with molecular dynamics (MD) or Monte Carlo. This efficiency enables long wave-length rippling to be captured and minimizes finite size effects that are important for the investigation of graphene properties. In this work, we used a graphene sheet of around 0.13 μm in each dimension. To the best of our knowledge, this size is at least a factor of 10 larger than the largest graphene ever simulated for finite size scaling properties.37 Graphene of similar size has recently been studied for vibrational spectral densities.38 The renormalization group derivation of the critical scaling exponent η assumes the membrane to be in the deep flat regime.25 The scaling behavior may change when graphene is getting closer to the crumpling temperature. Most existing numerical investigations of η are limited to lower temperatures at much smaller system sizes. In this work, we will provide accurate determination of η over a wide range of temperatures with a large graphene sheet. We also calculate the renormalized bending rigidities and normal-normal correlation functions of graphene to investigate the possible crumpling and normal persistence lengths that may exist for graphene. In addition, we investigate the thermal-expansion coefficient (TEC), radial distribution function (RDF), and translational and orientational correlation functions of graphene. The RDF and the translational and orientational correlation functions provide direct evidence regarding to the validity of the Mermin-Wagner instability. In this paper, the computational details will be described in Sec. II. Results and discussions are reported in Sec. III. We provide a summary in Sec. IV.
II. COMPUTATIONAL DETAILS
MD simulations were performed at 50, 100, 150, 200, 300, 500, 800, 1000, 1500, 2000, 2500, 3000, 3500, and 4000 K. At low temperature, quantum nuclear effect may be important and a classical simulation may become inappropriate. It is our assumption that graphene at 50 K can still be described with classical mechanics. This is likely a good approximation considering the heavy mass of a C atom. It is our assumption that defect density in graphene is low at 4000 K and the PPBE-G potential still remains sufficiently accurate. Although the PPBE-G potential allows breaking of C–C bond, its accuracy is questionable when graphene starts to disintegrate. Fortunately, the defect formation energy in graphene is relatively high at around 4.7 eV,23 which is much higher than the thermo-energy at 4000 K. A graphene sheet with 606 208 12 C atoms was used in our simulations. The initial configuration was perfectly flat with a surface area of 1282.3× 1284.0 Å2 . After equilibrium, the size of graphene decreased slightly due to rippling. The simulation was performed with the LAMMPS program39 with a time step of 1.0 fs below 1500 K and 0.5 fs at 1500 K or higher. A 3D periodic box was used with the graphene sheet placed in the x-y plane. The z dimension of the box was set to be 200.0 Å.
J. Chem. Phys. 141, 144701 (2014) TABLE I. Equilibrium surface areas of graphene at different temperatures. T (K)
Area (106 Å2 )
50 100 150 200 300 500 800 1000 1500 2000 2500 3000 3500 4000
1.6489 1.6483 1.6479 1.6474 1.6466 1.6452 1.6432 1.6424 1.6408 1.6389 1.6375 1.6364 1.6359 1.6351
In order to determine the equilibrium box size, each graphene sheet was simulated in the NSxy T ensemble for 9 ns with no strain in the x-y plane and the z dimension held as a constant. During the NSxy T simulation, the Nosé-Hoover thermostat and barostat were used with a temperature relaxation time of 0.5 ps and a pressure relaxation time of 1.0 ps. The equilibrium area of graphene was calculated by averaging the last 3 ns of the NSxy T simulation and the results were summarized in Table I. The graphene sheet contracts upon heating in the entire temperature range. A fit to the surface area in Table I yields a TEC of −2.7 × 10−6 K−1 at 300 K in good agreement with our previous estimate of −3.0 × 10−6 K−1 with this model. However, this value is about a factor of 2 smaller than current experimental estimates.6, 40 At 4000 K, the size of the graphene sheet reduces to approximately 1277.6 × 1279.8 Å2 . Since the sheet would be approximately 1282.3 × 1284.0 Å2 staying flat, the corrugation is much less than the z-dimension of the simulation box. After the NSxy T simulation, an additional 13 ns of MD was performed in the NVT ensemble for graphene sheets below 1000 K. 7 ns of NVT simulation was performed for graphene at or above 1000 K. The first 1 ns of NVT trajectories were discarded and the remaining trajectories were saved once every 10 ps for data analysis. In this work, the renormalized bending rigidity, κ, was calculated with the usual method26 through the power spectrum of the height correlation function: 2 ˜ = |h(q)|
kB T , κq 4
(6)
˜ where h(q) is the Fourier transform of the height h, kB is the Boltzmann constant; q in Eq. (6) is the norm of the wave vector q, which was calculated as q = (qx , qy ) = 2π (nx /Lx , ny /Ly ),
(7)
where Lx , Ly is the size of graphene and nx , ny are the grid indices of the discrete Fourier transform along the x- and ydirections, respectively. In this work, a 512 × 512 grid was used. Each grid point has on average 2.3 C atoms. The height h at each grid point was estimated as the average height of
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-4
D. Wei and F. Wang
FIG. 1. Thermal-expansion coefficient of graphene as a function of temperature. The solid line is a guide for the eyes.
the carbon atoms in the grid. To apply Eq. (6), the power spectrum was normalized with the area of the graphene.
III. RESULTS AND DISCUSSIONS
A recently survey of the temperature dependent thermalexpansion coefficient of graphene showed conflicting trends from atomistic scale simulations.41 Fig. 1 reports the TEC of graphene for the PPBE-G potential. The TEC of graphene is negative in the entire temperature range, consistent with a recent study by Pozzo et al. with PAW-PBE in a much smaller 512 C simulation box at up to 2000 K.42 Due to the existence of large scale ripples that contribute to the negative TEC,2, 16 large simulation boxes are required to properly converge TEC. This finite size effect was made clear in the Pozzo study, where a different trend was observed when smaller boxes were used. A negative graphene TEC persisting up to 2300 K was also predicted with the quasi-Harmonic approximation with the PBE functional.43 With a 8640 C atoms periodic box, previous simulations using the LCBOPII potential predicts the TEC of graphene to change from negative to positive when the temperature is higher than approximately 900 K.44 Due to contributions from long wave-length out-of-plane vibrations, a larger box is expected to allow more rippling thus leading to a more negative TEC. Our results are more consistent with predictions based on PBE. At the same time, the PPBE-G potential allows the use of box sizes three orders of magnitude larger than what is possible with PBE and may be able to better incorporate the contribution from low frequency vibrational modes to TEC. 2 ˜ ) as a function of log q at difFig. 2 reports log(|h(q)| ferent temperatures. A slope of approximately −4.0 can be observed at the range of −1.75 ≤ log q ≤ −0.25. Figure 3 reports κ as a function of temperature. κ increases monotonically from 1.22 eV at 50 K to 1.89 eV at 4000 K. In the inset of Fig. 3, κ from our model was compared to that from Fasolino et al.16 using the LCBOPII potential and by Liu et al. using the REBO-II potential.14 The PPBE-G potential showed a trend similar to that predicted by LCBOPII but with weaker temperature dependence. Only REBO-II predicted a behavior
J. Chem. Phys. 141, 144701 (2014)
2 ) as a function of log q for graphene at different tem˜ FIG. 2. Log(|h(q)| peratures. A linear region with a slope of −4.0 is observed at −1.75 ≤ log q ≤ −0.25.
consistent with Eq. (3). Thus the larger scale PPBE-G simulation also supports a strong anharmonic coupling in graphene consistent with previous predictions by Fasolino et al.16 Although Fig. 2 does not rule out the possibility that a very large sheet may still crumple, the estimation of the crumpling persistence length according to Eq. (4) is not appropriate for graphene. In order to understand long-range crystallinity, the C-C RDF of graphene is reported in Figure 4. A crystal will have a RDF that is always oscillating. The RDF of a liquid generally decays to one after a few solvation shells. As a 2D membrane, the RDF of graphene does remain oscillatory at separations as large as 600 Å. The oscillatory behavior at large separation is observed even for the highest temperature investigated. In order to further investigate the decay of structural order in graphene, we report the reduced peak height, which is defined as the absolute difference between the peak height
FIG. 3. Renormalized bending rigidity, κ, of graphene as a function of temperature. The solid line is just a guide for the eyes. The inset compares the κ calculated by the PPE-G potential (this work) to that obtained with the LCBOPII potential16, 47 and the REBOII potential.14
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-5
D. Wei and F. Wang
J. Chem. Phys. 141, 144701 (2014)
FIG. 4. Radial distribution function (RDF) of graphene at different temperatures. The RDF from 620 to 640 Å is enlarged and shown in panel (b).
in RDF and one. For a liquid, the reduced peak height will be zero as the liquid loses long-range order. It is clear from Figure 5 that the highest peaks for the reduced peak height scale as (1/R)δ , with δ being 0.5 at 100 K and 0.8 at 3000 K. To the best of our knowledge, such a scaling behavior for graphene has not been reported previously. Due to the inverse power law scaling, the crystallinity of graphene will be lost in all practical sense at large separations consistent with the prediction of the Mermin–Wagner instability.
FIG. 6. Log-log plots of the translational correlation function GT (R) (a) and the orientational correlation function G6 (R) (b) of graphene at 300, 1000, 3000, and 4000 K.
We note such a power law scaling is similar to that observed for 2D crystals modeled with a square shoulder square well (SSSW) potential, where a δ of 0.7 was observed with reduced temperature of 0.3 and 0.4.45 One key difference between SSSW and graphene is that graphene is a 2D crystal in 3D space and the previous work with the SSSW potential is strictly a 2D system. Figure 6 reports the translational correlation function GT (R) and orientational correlation functions G6 (R)45, 46 of graphene at temperatures 300, 1000, 3000, and 4000 K. We followed the definition: GT (R) = |exp(iG0 · Rj )|,
(8)
where G0 is a vector of the reciprocal lattice and Rj with magnitude R is the position vector of atom j relative to an origin taken to be one of the atom positions,45 and G6 (R) = |Q6 (R)Q∗6 (0)|,
(9)
with N
FIG. 5. Log-log plot of the reduced peak height of the C-C RDF as a function of distance at 100 K (a), 300 K (b), 1000 K (c), and 3000 K (d).
j 1 Q6 (Rj ) = exp(6iθj k ), Nj k
(10)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-6
D. Wei and F. Wang
FIG. 7. (a) Normal-normal correlation function G(R) at different temperatures calculated with the 0.13 μm × 0.13 μm box. (b) G(R) calculated for three different box sizes: A (225.4 × 234.3 Å2 ), B (641.2 × 642.0 Å2 ), and C (1282.3 × 1284.0 Å2 ) at 300, 1000, and 3000 K show almost identical crossover from fast decay to slow decay.
where Nj is the number of neighbors around atom j, which is 3 for graphene; θ jk is the angle of the bond between neighboring atoms j and k measured with respect to an arbitrary reference axis. In the definitions of GT (R) and G6 (R), the angle brackets indicate averages and the correlation functions were defined based on the amplitude of complex averages. For a perfect crystal, both GT (R) and G6 (R) will reach one. Figure 6 shows that GT (R) undergoes an inverse powerlaw decay, which is consistent with graphene losing long range translational order as a result of the Mermin-Wagner instability. Similar to the reduced peak height, the exponent for the GT (R) decay also depends on temperature, ranging from 0.49 at 300 K, 0.57 at 1000 K, to 0.62 at 3000 K and 4000 K. Contrary to GT (R), no appreciably decay of G6 (R) can be observed for graphene beyond 30 Å, indicating that graphene retains long range orientational order. Fig. 7(a) reports the normal-normal correlation function G(R) = n(0) · n(R) at different temperatures. The normal vector at each C atom was determined as the average of the normals of three planes formed by the central C atom and its three nearest neighbors.47 From Fig. 7(a), the G(R) shows a fast initial decay followed by a very slow decay. The slow de-
J. Chem. Phys. 141, 144701 (2014)
cay does not behave exponentially. At lower T, ξ n is expected to be large. It is indeed possible that a slow exponential decay may not look exponential over a short distance. However, the tail of G(R) does not look exponential even at 4000 K, where ξ n is expected to be comparable to the size of our simulation box. Thus, we conclude that G(R) of graphene does not exhibit an exponential decay. Consequently, the definition of normal persistence length is not appropriate for graphene. Since the onset of the slow decay is less than a tenth of the size of our simulation box, it is unlikely the behavior is caused by finite size limitations. Nonetheless, simulations were repeated at 3 temperatures with two smaller box sizes of 222.3 × 231.0 and 632.3 × 633.1 Å2 , respectively, when flat. As shown in Fig. 7(b), the G(R) for the smaller boxes switches to the slow decay at approximately the same position. We are thus certain the behavior is not a result of finite size effects. In order to better understand the slow decay, the power spectrum of G(R) was reported in Fig. 8 and is designated ˜ as G(q). For tethered membrane in a deep flat state, it is an˜ ticipated that G(q) should satisfy a power law scaling with −(2−η) ˜ at small q (large separations).27, 37 The critG(q) ∝ q ical exponent, η, is expected to be 0.849 at low enough temperature according to renormalization group theory, although other estimates ranging from 0.75 to 0.82 has been proposed.26, 27, 29 Clearly, a linear scaling can be observed in Fig. 8. It is interesting to note that the critical exponent η is approximately 0.85 for all the temperatures investigated up to 4000 K, whereas the renormalization group estimate of 0.849 is expected to hold only for membrane far from the crumpling transition. The η of 0.85 is in good agreement with earlier numerical estimate obtained by Los37 with the LCBOPII potential. We note the LCBOPII potential showed different temperature dependence of bending rigidity and TEC. The insensitivity of η with respect to the interaction potential is anticipated for a universal scaling behavior. Fig. 9(a) reports the height-height correlation function H(R) = [h(R) − h(0)]2 at different temperatures. The increase of H(R) slows down toward 620 Å, which is approximately half box size. Fig. 9(b) shows a log-log plot of H(R)
FIG. 8. Log-log plot of the power spectrum of the normal-normal correlation ˜ ˜ function, G(q). At small q, G(q) ∝ q −(2−η) , with a universal scaling constant η = 0.85.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-7
D. Wei and F. Wang
FIG. 9. (a) Height-height correlation function H(R) as a function of distance. (b) Log-log plots of H(R) at the distances from 20.0 to 400.0 Å. A dashed line with a slope of 1.15 is shown for reference.
as a function of distance. A linear region with a slope of 1.15 can be seen when R is no more than approximately 150 Å. The deviation from linear scaling at large R is most likely caused by finite size effects. A slope of 1.15 indicates that ˜ this universal H(R) ∝ R2−η , with η = 0.85. As for G(q), scaling constant holds for graphene up to 4000 K. The exponent is in good agreement with the finding of Los et al.37 obtained with a smaller box at 300 K. A more detailed view is available from Fig. 7(b) since a denser distance grid was used in this work when compared with the Los study.
IV. SUMMARY
We performed MD simulations on large graphene sheets with the realistic and efficient PPBE-G potential recently developed by force matching PAW-PBE at finite temperatures. Simulation with the PPBE-G potential is the best approximation to a direct CPMD simulation with PAW-PBE. Large simulation boxes allow better modeling of long wave-length oscillations that are important for converging TEC and the critical scaling exponent η at low q. To the best of our knowledge, our simulation box is significantly larger than any prior simulations on either TEC or the critical scaling properties
J. Chem. Phys. 141, 144701 (2014)
of graphene. The temperature dependence of η has not been investigated prior to our work. The TEC of graphene remains negative for the entire temperature range investigated consistent with previous studies with the PBE functional at temperatures up to 2000 K. Although the RDF of graphene is still oscillatory at 620 Å, the reduced peak heights of RDF decay as inverse power law to the distance. Such a decay is consistent with the requirement of the Mermin-Wagner instability. Although graphene loses long range translational order as indicated by an inverse power law decay of GT (R), graphene does possess long range orientational order. The PPBE-G potential predicts the renormalized bending rigidity to increase with temperature consistent with earlier simulations with the LCBOPII potential,16, 18 which was contradictory to the REBO-II prediction.14 The normal-normal correlation of graphene G(R) does not show an exponential decay. This behavior makes the standard definition of normal persistence length lose its significance for graphene. Consistent with previous studies, even with much larger simulation boxes, the power spectrum of G(R) decays according inverse power law with a critical exponent of –(2 − η), with η = 0.85. Assuming our potential remains reliable up to the highest temperature simulated, the critical exponent η of 0.85 is a constant for graphene up to 4000 K, indicating graphene to be in the deep flat regime even at less than 1000 K below the experimental melting temperature. It is thus intriguing whether graphene has a crumpling transition before it melts. Graphene is a tethered membrane that loses translational order at large R. By losing its long range translational order, graphene is no longer periodic and does not satisfy the general accepted definition of a crystal in 3D. Although being referred to as a crystal in 2D, graphene does not possess conventional crystalline long-range order. At the same time, graphene shows strong anharmonic coupling and does not follow thermodynamics of a liquid membrane. We note that one limitation of a force-field based MD simulation is the inability to treat electronic excitations that may be important for large graphene sheets. This is a limitation shared also by DFT based ab initio MD simulations, such as CPMD. With AFM based force fields, we can map DFT potential energy surfaces to less computational intensive formulas; however, the accuracy of an AFM force field will be bound by the reference method, in this case Kohn-Sham DFT. On the other hand, the scaling properties of graphene shows universality class behavior, we do not anticipate minor inaccuracy of the PBE functional or explicit treatment of electronic excitation will change the overall scaling. Nonetheless, developing more accurate electronic structure models for graphene and explicit inclusion of phonon-exciton coupling are important future directions of research. ACKNOWLEDGMENTS
This work is supported by key scientific and technological projects of Chongqing, China (cstc2012ggC50002) and NSF CAREER Award No. CHE0748628. The computer resource for the study was provided by the High Performance Center of University of Arkansas and the Super-Computer
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-8
D. Wei and F. Wang
Center of Chongqing Institution of Green and Intelligent Technology, Chinese Academy of Sciences. F.W. acknowledges support from the Arkansas Biosciences Institute. 1 K.
S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306(5696), 666 (2004). 2 J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature (London) 446(7131), 60 (2007). 3 F. Guinea, B. Horovitz, and P. Le Doussal, Phys. Rev. B 77(20), 205421 (2008). 4 R. Miranda and A. L. Vazquez de Parga, Nat. Nano 4(9), 549 (2009). 5 A. L. Vázquez de Parga, F. Calleja, B. Borca, M. C. G. Passeggi, J. J. Hinarejos, F. Guinea, and R. Miranda, Phys. Rev. Lett. 100(5), 056807 (2008). 6 W. Z. Bao, F. Miao, Z. Chen, H. Zhang, W. Y. Jang, C. Dames, and C. N. Lau, Nat. Nanotechnol. 4(9), 562 (2009). 7 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81(1), 109 (2009). 8 M. I. Katsnelson and A. K. Geim, Philos. Trans. Royal Soc., A 366(1863), 195 (2008). 9 R. Peierls, Ann. Inst. Henri Poincare 5, 177 (1935). 10 L. Landau, Phys. Z. Sowjetunion 11, 26 (1937). 11 L. Peliti and S. Leibler, Phys. Rev. Lett. 54(15), 1690 (1985). 12 N. D. Mermin, Phys. Rev. 176(1), 250 (1968). 13 D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J. Phys.: Condens. Matter 14(4), 783 (2002). 14 P. Liu and Y. W. Zhang, Appl. Phys. Lett. 94(23), 231912 (2009). 15 J. H. Los, L. M. Ghiringhelli, E. J. Meijer, and A. Fasolino, Phys. Rev. B 72(21), 214102 (2005). 16 A. Fasolino, J. H. Los, and M. I. Katsnelson, Nat. Mater. 6(11), 858 (2007). 17 P. G. De Gennes and C. Taupin, J. Phys. Chem. 86(13), 2294 (1982). 18 M. I. Katsnelson and A. Fasolino, Acc. Chem. Res. 46(1), 97 (2013). 19 K. N. Kudin, G. E. Scuseria, and B. I. Yakobson, Phys. Rev. B 64(23), 235406 (2001). 20 R. C. Thompson-Flagg, M. J. B. Moura, and M. Marder, Europhys. Lett. 85(4), 46002 (2009).
J. Chem. Phys. 141, 144701 (2014) 21 F.
L. Braghin and N. Hasselmann, Phys. Rev. B 82(3), 035407 (2010). I. Savvatimskiy, Carbon 43(6), 1115 (2005). 23 J. Ma, D. Alfè, A. Michaelides, and E. Wang, Phys. Rev. B 80(3), 033407 (2009). 24 Y. Kantor, M. Kardar, and D. R. Nelson, Phys. Rev. Lett. 57(7), 791 (1986). 25 J. P. Kownacki and D. Mouhanna, Phys. Rev. E 79(4), 040101 (2009). 26 P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett. 69(8), 1209 (1992). 27 Z. Zhang, H. T. Davis, and D. M. Kroll, Phys. Rev. E 48(2), R651 (1993). 28 M. J. Bowick and A. Travesset, Phys. Rep. 344(4–6), 255 (2001). 29 M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleifsson, and K. N. Anagnostopoulos, J. Phys. I (France) 6(10), 1321 (1996). 30 D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134(18), 184704 (2011). 31 O. Akin-Ojo, Y. Song, and F. Wang, J. Chem. Phys. 129(6), 064108 (2008). 32 O. Akin-Ojo and F. Wang, J. Comput. Chem. 32, 453 (2011). 33 O. Akin-Ojo and F. Wang, J. Phys. Chem. B 113(5), 1237 (2009). 34 P. E. Blöchl, Phys. Rev. B 50(24), 17953 (1994). 35 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77(18), 3865 (1996). 36 K. V. Zakharchenko, R. Roldán, A. Fasolino, and M. I. Katsnelson, Phys. Rev. B 82(12), 125435 (2010). 37 J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Zakharchenko, and A. Fasolino, Phys. Rev. B 80(12), 121405 (2009). 38 D. Liesegang and C. Oligschleger, J. Mod. Phys. 5, 149 (2014). 39 S. Plimpton, J. Comput. Phys. 117(1), 1 (1995). 40 D. Yoon, Y.-W. Son, and H. Cheong, Nano Lett. 11(8), 3227 (2011). 41 Y. Magnin, G. D. Förster, F. Rabilloud, F. Calvo, A. Zappelli, and C. Bichara, J. Phys.: Condens. Matter 26(18), 185401 (2014). 42 M. Pozzo, D. Alfè, P. Lacovig, P. Hofmann, S. Lizzit, and A. Baraldi, Phys. Rev. Lett. 106(13), 135501 (2011). 43 N. Mounet and N. Marzari, Phys. Rev. B 71(20), 205214 (2005). 44 K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino, Phys. Rev. Lett. 102(4), 046808 (2009). 45 A. M. Almudallal, S. V. Buldyrev, and I. Saika-Voivod, J. Chem. Phys. 140(14), 144505 (2014). 46 U. Gasser, J. Phys.: Condens. Matter 21(20), 203101 (2009). 47 K. V. Zakharchenko, J. H. Los, M. I. Katsnelson, and A. Fasolino, Phys. Rev. B 81(23), 235439 (2010). 22 A.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
THE JOURNAL OF CHEMICAL PHYSICS 141, 144701 (2014)
Graphene: A partially ordered non-periodic solid Dongshan Wei1 and Feng Wang2,a) 1
Key Laboratory of Multi-Scale Manufacturing Technology, Chongqing Institute of Green and Intelligent Technology, Chinese Academy of Sciences, Chongqing 400714, China 2 Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, USA
(Received 2 June 2014; accepted 24 September 2014; published online 9 October 2014) Molecular dynamics simulations were performed to study the structural features of graphene over a wide range of temperatures from 50 to 4000 K using the PPBE-G potential [D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134, 184704 (2011)]. This potential was developed by force matching the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional and has been validated previously to provide accurate potential energy surface for graphene at temperatures as high as 3000 K. Simulations with the PPBE–G potential are the best available approximation to a direct Car-Parrinello Molecular Dynamics study of graphene. One advantage of the PBE-G potential is to allow large simulation boxes to be modeled efficiently so that properties showing strong finite size effects can be studied. Our simulation box contains more than 600 000 C atoms and is one of the largest graphene boxes ever modeled. With the PPBE-G potential, the thermal-expansion coefficient is negative up to 4000 K. With a large box and an accurate potential, the critical exponent for the scaling properties associated with the normal-normal and height-height correlation functions was confirmed to be 0.85. This exponent remains constant up to 4000 K suggesting graphene to be in the deeply cooled regime even close to the experimental melting temperature. The reduced peak heights in the radial distribution function of graphene show an inverse power law dependence to distance, which indicates that a macroscopic graphene sheet will lose long-range crystalline order as predicted by the Mermin-Wagner instability. Although graphene loses long-range translational order, it retains long range orientational order as indicated by its orientational correlation function; graphene is thus partially ordered but not periodic. © 2014 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http://dx.doi.org/10.1063/1.4897255] I. INTRODUCTION
Graphene is a two-dimensional (2D) material that has the structure of a graphite monolayer.1 Perfectly flat graphene will have a zero band gap. This zero band gap leads to many superior electronic properties, making graphene a promising material for building ground breaking electronics. However, at any finite temperature, a graphene sheet cannot be perfectly flat. Spontaneous rippling of the sheet has been observed experimentally.2 Such rippling will cause the bandgap to deviate from zero3–6 and produce electron-ripple scattering.7, 8 In order to fully understand the potential of graphene as a new material, it is important to investigate the effect of rippling for large graphene sheets. Due to thermo-fluctuations, no 2D material can be perfectly flat at any finite temperature. Ignoring the surface tension term, the free energy of a 2D membrane can be written as 1 dxdy(∇ 2 h)2 F = κ 2 1 (1) + dxdy λ u2xx + u2yy + 2μu2xy , 2
a) [email protected]
0021-9606/2014/141(14)/144701/8
where h is the normal displacement, κ is the bending rigidity, λ and μ are Lamé’s first and second parameters, respectively. In Eq. (1), the strain matrix elements uαβ are defined as ∂uβ ∂h ∂h 1 ∂uα + + uαβ ≈ , (2) 2 ∂β ∂α ∂α ∂β where uα and uβ are in-plane displacements. Equation (1) is applicable to both solid and liquid membranes. The experimental observation of large graphene sheets surprised the scientific community initially. This was due to the argument that no 2D material can exist over a large length scale. It is worth noting that there are two similar arguments regarding the stability of 2D materials, one by Peierls9 and Landau10 and another one by Peliti and Leibler.11 The PeierlsLandau argument states that there could be no long range crystalline order in 1D and 2D materials. This is generally referred to as the Mermin-Wagner instability.12 This statement allows the existence of 1D or 2D materials but requires such materials not to be periodic at macroscopic scale. In other words, such reduced dimensional material can exist but lacks long range translational order. The argument by Peliti and Leibler is stronger. They argues that no 2D material can exisit at macrosopic scale. This statement was made by assuming the second integral in Eq. (1) to vanish, which is true when the coupling between
141, 144701-1
© Author(s) 2014
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-2
D. Wei and F. Wang
J. Chem. Phys. 141, 144701 (2014)
normal and in-plane displacements is negligible. Under this condition, one can show that the renormalized bending rigidity depends on membrane size and to the first order can be described by the equation, κ(L) = κ0 − (3kB T /4π ) ln(L/a),
(3)
where κ 0 is the bare bending rigidity and a is a certain microscopic length. In this scenario, the bending rigidity of the 2D material diverges to negative infinity when its size is asymptotically large; the material spontaneously wraps around to form 3D structures, such as vesicles. The assumption made by Peliti and Leibler holds for liquid membranes. A solid membrane can potentially satisfy Eq. (3) if the coupling between in-plane and out-of-plane displacement is weak. The renormalized bending rigidity of a material that satisfies Eq. (3) will have a negative temperature dependence. Whether graphene satisfies Eq. (3) is a subject of debate. Simulations using the REBO-II potential13 shows negative temperature dependence14 of κ(L) consistent with Eq. (3) while the LCBOPII potential15 produces a positive dependence.16 For 2D materials that do satisfy Eq. (3), the length at which κ(L) vanishes is frequently referred to as the persistence length ξ .11 ξ = a · e4πκ0 /3kB T .
(4)
Such a material will undergo the crumpling phase transition when its size is larger than ξ . Another popular definition for the persistence length is through the normalnormal correlation function G(R) = n(0) · n(R). For weakly coupled membranes, the normal-normal correlation function is expected to decay exponentially with respect to the separation R, n(0) · n(R) = e−R/ξn .
(5)
The decay rate can be fit by a single exponential to obtain the persistence length ξ n of the normal. For a liquid membrane, such as an oil-water interface saturated with surfac2πκ 0
tants, De Gennes and Taupin17 argue that ξn = a · e 3kB T thus 2πκ 0
ξ = ξn · e 3kB T or ξ = ξn2 /a. In the following discussion, we will refer to ξ as the crumpling persistence length and ξ n as the normal persistence length. Whether graphene is susceptible to a crumpling phase transition depends on the model potential being used. For example, the LCBOPII potential16, 18 shows a complete breakdown of Eq. (3), but simulation with the REBO-II potential14 shows a behavior that is compatible to Eq. (3). Although one could argue that the experimental observation of large suspended graphene will rule out the crumpling instability predicted in Eq. (3), this is not true. The bare bending rigidity of graphene has been estimated to be in the range of 1 eV or higher.16, 19, 20 Even if the thermodynamics of graphene were governed by Eq. (3), at 300 K, with an a of 1.4 Å, Eq. (4) predicts a ξ of 3.3 × 1060 m and a ξ n of 2.2 × 1025 m, both are astronomically large. At 3000 K, ξ would still be microscopically large at 1.5 mm and ξ n would be 0.46 μm. We note referring to graphene as a solid membrane or a membrane with local crystalline order is probably correct
but not very precise. A lipid membrane can be in a 2D crystalline state at low temperatures. Although there is local crystalline order in such a state, there are no covalent bonds between lipid molecules. Graphene is more precisely referred to as a tethered surface or polymerized membrane.21 For a tethered membrane, the bonding connectivity between neighboring molecules cannot change. This is true for graphene under a wide range of temperatures. Experimentally, graphite will not melt until around 4800 K.22 The Stone-Wales topological defect is expected to have a formation energy of 4.7 eV.23 Even at 4000 K, the thermo-population of such defects is around 10−6 . The difference between a tethered membrane with fixed connectivity and a free diffusing membrane with weak intermolecular interactions can be illustrated by a vivid example given by Kantor.24 Although a fluid membrane can crumple and form a sphere, such as a vesicle, a tethered membrane, such as a piece of paper, can never form a sphere. The overall “shape” of a tethered surface cannot change. The Lamé’s parameters for a tethered surface can be large since bending to the third dimension cannot be accomplished without causing strain in the in-plane direction by stretching and compressing covalent bonds. The thermodynamics of a tethered surface is far from well-understood. It has been argued that a tethered membrane will crumple when κ is less than a critical value κ c . When κ is sufficiently large, a tethered membrane will remain macroscopically flat. We note having the tethered membrane to remain macroscopically flat does not contradict with the Mermin–Wagner instability.12 As mentioned previously, the Mermin-Wagner instability only rules out long range translational order but allows the existence of 2D materials. Another interesting property that can be derived for nonself-avoiding (phantom) flat membrane is the existence of a universal scaling factor as a function of size for the power spectra of normal-normal correlation function and the heightheight correlation function. The existence of universal scaling constant is frequently observed in 3D materials near a critical point. A universality class behavior in 2D away from a critical point is very intriguing.18 Several theoretical25, 26 and numerical27–29 investigations of the critical scaling exponent, η, has been derived for 2D membrane in 3D space. The most accurate estimate available is probably η = 0.849 from renormalization group theory,25 although other estimates such as 0.821,26 0.750,29 and 0.8127 have also been reported. In this work, we investigate the thermodynamics of a large graphene sheet using the recently developed PPBEG potential.30 The PPBE-G potential was created by forcematching30–33 the Perdew-Burke-Ernzerhof (PBE) forces calculated with the Projector Augmented Wave (PAW) treatment in the temperature range from 300 to 3000 K.34, 35 In the temperature range of parameterization, the PPBE-G potential produces an atomic force root-mean-square-error (RMSE) of 4.4 kcal/(mol Å) when compared to PAW-PBE references.30 Since other graphene force fields surveyed produce a force RMSE of at least 14 kcal/(mol Å) when compared to PAWPBE, a PPBE-G simulation is the best available approximation to a direct Car-Parrinello molecular dynamics36 (CPMD) study for graphene.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-3
D. Wei and F. Wang
One clear advantage of the PPBE-G potential is its high computational efficiency. This efficiency allows graphene sheets with millions of atoms to be studied with molecular dynamics (MD) or Monte Carlo. This efficiency enables long wave-length rippling to be captured and minimizes finite size effects that are important for the investigation of graphene properties. In this work, we used a graphene sheet of around 0.13 μm in each dimension. To the best of our knowledge, this size is at least a factor of 10 larger than the largest graphene ever simulated for finite size scaling properties.37 Graphene of similar size has recently been studied for vibrational spectral densities.38 The renormalization group derivation of the critical scaling exponent η assumes the membrane to be in the deep flat regime.25 The scaling behavior may change when graphene is getting closer to the crumpling temperature. Most existing numerical investigations of η are limited to lower temperatures at much smaller system sizes. In this work, we will provide accurate determination of η over a wide range of temperatures with a large graphene sheet. We also calculate the renormalized bending rigidities and normal-normal correlation functions of graphene to investigate the possible crumpling and normal persistence lengths that may exist for graphene. In addition, we investigate the thermal-expansion coefficient (TEC), radial distribution function (RDF), and translational and orientational correlation functions of graphene. The RDF and the translational and orientational correlation functions provide direct evidence regarding to the validity of the Mermin-Wagner instability. In this paper, the computational details will be described in Sec. II. Results and discussions are reported in Sec. III. We provide a summary in Sec. IV.
II. COMPUTATIONAL DETAILS
MD simulations were performed at 50, 100, 150, 200, 300, 500, 800, 1000, 1500, 2000, 2500, 3000, 3500, and 4000 K. At low temperature, quantum nuclear effect may be important and a classical simulation may become inappropriate. It is our assumption that graphene at 50 K can still be described with classical mechanics. This is likely a good approximation considering the heavy mass of a C atom. It is our assumption that defect density in graphene is low at 4000 K and the PPBE-G potential still remains sufficiently accurate. Although the PPBE-G potential allows breaking of C–C bond, its accuracy is questionable when graphene starts to disintegrate. Fortunately, the defect formation energy in graphene is relatively high at around 4.7 eV,23 which is much higher than the thermo-energy at 4000 K. A graphene sheet with 606 208 12 C atoms was used in our simulations. The initial configuration was perfectly flat with a surface area of 1282.3× 1284.0 Å2 . After equilibrium, the size of graphene decreased slightly due to rippling. The simulation was performed with the LAMMPS program39 with a time step of 1.0 fs below 1500 K and 0.5 fs at 1500 K or higher. A 3D periodic box was used with the graphene sheet placed in the x-y plane. The z dimension of the box was set to be 200.0 Å.
J. Chem. Phys. 141, 144701 (2014) TABLE I. Equilibrium surface areas of graphene at different temperatures. T (K)
Area (106 Å2 )
50 100 150 200 300 500 800 1000 1500 2000 2500 3000 3500 4000
1.6489 1.6483 1.6479 1.6474 1.6466 1.6452 1.6432 1.6424 1.6408 1.6389 1.6375 1.6364 1.6359 1.6351
In order to determine the equilibrium box size, each graphene sheet was simulated in the NSxy T ensemble for 9 ns with no strain in the x-y plane and the z dimension held as a constant. During the NSxy T simulation, the Nosé-Hoover thermostat and barostat were used with a temperature relaxation time of 0.5 ps and a pressure relaxation time of 1.0 ps. The equilibrium area of graphene was calculated by averaging the last 3 ns of the NSxy T simulation and the results were summarized in Table I. The graphene sheet contracts upon heating in the entire temperature range. A fit to the surface area in Table I yields a TEC of −2.7 × 10−6 K−1 at 300 K in good agreement with our previous estimate of −3.0 × 10−6 K−1 with this model. However, this value is about a factor of 2 smaller than current experimental estimates.6, 40 At 4000 K, the size of the graphene sheet reduces to approximately 1277.6 × 1279.8 Å2 . Since the sheet would be approximately 1282.3 × 1284.0 Å2 staying flat, the corrugation is much less than the z-dimension of the simulation box. After the NSxy T simulation, an additional 13 ns of MD was performed in the NVT ensemble for graphene sheets below 1000 K. 7 ns of NVT simulation was performed for graphene at or above 1000 K. The first 1 ns of NVT trajectories were discarded and the remaining trajectories were saved once every 10 ps for data analysis. In this work, the renormalized bending rigidity, κ, was calculated with the usual method26 through the power spectrum of the height correlation function: 2 ˜ = |h(q)|
kB T , κq 4
(6)
˜ where h(q) is the Fourier transform of the height h, kB is the Boltzmann constant; q in Eq. (6) is the norm of the wave vector q, which was calculated as q = (qx , qy ) = 2π (nx /Lx , ny /Ly ),
(7)
where Lx , Ly is the size of graphene and nx , ny are the grid indices of the discrete Fourier transform along the x- and ydirections, respectively. In this work, a 512 × 512 grid was used. Each grid point has on average 2.3 C atoms. The height h at each grid point was estimated as the average height of
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-4
D. Wei and F. Wang
FIG. 1. Thermal-expansion coefficient of graphene as a function of temperature. The solid line is a guide for the eyes.
the carbon atoms in the grid. To apply Eq. (6), the power spectrum was normalized with the area of the graphene.
III. RESULTS AND DISCUSSIONS
A recently survey of the temperature dependent thermalexpansion coefficient of graphene showed conflicting trends from atomistic scale simulations.41 Fig. 1 reports the TEC of graphene for the PPBE-G potential. The TEC of graphene is negative in the entire temperature range, consistent with a recent study by Pozzo et al. with PAW-PBE in a much smaller 512 C simulation box at up to 2000 K.42 Due to the existence of large scale ripples that contribute to the negative TEC,2, 16 large simulation boxes are required to properly converge TEC. This finite size effect was made clear in the Pozzo study, where a different trend was observed when smaller boxes were used. A negative graphene TEC persisting up to 2300 K was also predicted with the quasi-Harmonic approximation with the PBE functional.43 With a 8640 C atoms periodic box, previous simulations using the LCBOPII potential predicts the TEC of graphene to change from negative to positive when the temperature is higher than approximately 900 K.44 Due to contributions from long wave-length out-of-plane vibrations, a larger box is expected to allow more rippling thus leading to a more negative TEC. Our results are more consistent with predictions based on PBE. At the same time, the PPBE-G potential allows the use of box sizes three orders of magnitude larger than what is possible with PBE and may be able to better incorporate the contribution from low frequency vibrational modes to TEC. 2 ˜ ) as a function of log q at difFig. 2 reports log(|h(q)| ferent temperatures. A slope of approximately −4.0 can be observed at the range of −1.75 ≤ log q ≤ −0.25. Figure 3 reports κ as a function of temperature. κ increases monotonically from 1.22 eV at 50 K to 1.89 eV at 4000 K. In the inset of Fig. 3, κ from our model was compared to that from Fasolino et al.16 using the LCBOPII potential and by Liu et al. using the REBO-II potential.14 The PPBE-G potential showed a trend similar to that predicted by LCBOPII but with weaker temperature dependence. Only REBO-II predicted a behavior
J. Chem. Phys. 141, 144701 (2014)
2 ) as a function of log q for graphene at different tem˜ FIG. 2. Log(|h(q)| peratures. A linear region with a slope of −4.0 is observed at −1.75 ≤ log q ≤ −0.25.
consistent with Eq. (3). Thus the larger scale PPBE-G simulation also supports a strong anharmonic coupling in graphene consistent with previous predictions by Fasolino et al.16 Although Fig. 2 does not rule out the possibility that a very large sheet may still crumple, the estimation of the crumpling persistence length according to Eq. (4) is not appropriate for graphene. In order to understand long-range crystallinity, the C-C RDF of graphene is reported in Figure 4. A crystal will have a RDF that is always oscillating. The RDF of a liquid generally decays to one after a few solvation shells. As a 2D membrane, the RDF of graphene does remain oscillatory at separations as large as 600 Å. The oscillatory behavior at large separation is observed even for the highest temperature investigated. In order to further investigate the decay of structural order in graphene, we report the reduced peak height, which is defined as the absolute difference between the peak height
FIG. 3. Renormalized bending rigidity, κ, of graphene as a function of temperature. The solid line is just a guide for the eyes. The inset compares the κ calculated by the PPE-G potential (this work) to that obtained with the LCBOPII potential16, 47 and the REBOII potential.14
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-5
D. Wei and F. Wang
J. Chem. Phys. 141, 144701 (2014)
FIG. 4. Radial distribution function (RDF) of graphene at different temperatures. The RDF from 620 to 640 Å is enlarged and shown in panel (b).
in RDF and one. For a liquid, the reduced peak height will be zero as the liquid loses long-range order. It is clear from Figure 5 that the highest peaks for the reduced peak height scale as (1/R)δ , with δ being 0.5 at 100 K and 0.8 at 3000 K. To the best of our knowledge, such a scaling behavior for graphene has not been reported previously. Due to the inverse power law scaling, the crystallinity of graphene will be lost in all practical sense at large separations consistent with the prediction of the Mermin–Wagner instability.
FIG. 6. Log-log plots of the translational correlation function GT (R) (a) and the orientational correlation function G6 (R) (b) of graphene at 300, 1000, 3000, and 4000 K.
We note such a power law scaling is similar to that observed for 2D crystals modeled with a square shoulder square well (SSSW) potential, where a δ of 0.7 was observed with reduced temperature of 0.3 and 0.4.45 One key difference between SSSW and graphene is that graphene is a 2D crystal in 3D space and the previous work with the SSSW potential is strictly a 2D system. Figure 6 reports the translational correlation function GT (R) and orientational correlation functions G6 (R)45, 46 of graphene at temperatures 300, 1000, 3000, and 4000 K. We followed the definition: GT (R) = |exp(iG0 · Rj )|,
(8)
where G0 is a vector of the reciprocal lattice and Rj with magnitude R is the position vector of atom j relative to an origin taken to be one of the atom positions,45 and G6 (R) = |Q6 (R)Q∗6 (0)|,
(9)
with N
FIG. 5. Log-log plot of the reduced peak height of the C-C RDF as a function of distance at 100 K (a), 300 K (b), 1000 K (c), and 3000 K (d).
j 1 Q6 (Rj ) = exp(6iθj k ), Nj k
(10)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-6
D. Wei and F. Wang
FIG. 7. (a) Normal-normal correlation function G(R) at different temperatures calculated with the 0.13 μm × 0.13 μm box. (b) G(R) calculated for three different box sizes: A (225.4 × 234.3 Å2 ), B (641.2 × 642.0 Å2 ), and C (1282.3 × 1284.0 Å2 ) at 300, 1000, and 3000 K show almost identical crossover from fast decay to slow decay.
where Nj is the number of neighbors around atom j, which is 3 for graphene; θ jk is the angle of the bond between neighboring atoms j and k measured with respect to an arbitrary reference axis. In the definitions of GT (R) and G6 (R), the angle brackets indicate averages and the correlation functions were defined based on the amplitude of complex averages. For a perfect crystal, both GT (R) and G6 (R) will reach one. Figure 6 shows that GT (R) undergoes an inverse powerlaw decay, which is consistent with graphene losing long range translational order as a result of the Mermin-Wagner instability. Similar to the reduced peak height, the exponent for the GT (R) decay also depends on temperature, ranging from 0.49 at 300 K, 0.57 at 1000 K, to 0.62 at 3000 K and 4000 K. Contrary to GT (R), no appreciably decay of G6 (R) can be observed for graphene beyond 30 Å, indicating that graphene retains long range orientational order. Fig. 7(a) reports the normal-normal correlation function G(R) = n(0) · n(R) at different temperatures. The normal vector at each C atom was determined as the average of the normals of three planes formed by the central C atom and its three nearest neighbors.47 From Fig. 7(a), the G(R) shows a fast initial decay followed by a very slow decay. The slow de-
J. Chem. Phys. 141, 144701 (2014)
cay does not behave exponentially. At lower T, ξ n is expected to be large. It is indeed possible that a slow exponential decay may not look exponential over a short distance. However, the tail of G(R) does not look exponential even at 4000 K, where ξ n is expected to be comparable to the size of our simulation box. Thus, we conclude that G(R) of graphene does not exhibit an exponential decay. Consequently, the definition of normal persistence length is not appropriate for graphene. Since the onset of the slow decay is less than a tenth of the size of our simulation box, it is unlikely the behavior is caused by finite size limitations. Nonetheless, simulations were repeated at 3 temperatures with two smaller box sizes of 222.3 × 231.0 and 632.3 × 633.1 Å2 , respectively, when flat. As shown in Fig. 7(b), the G(R) for the smaller boxes switches to the slow decay at approximately the same position. We are thus certain the behavior is not a result of finite size effects. In order to better understand the slow decay, the power spectrum of G(R) was reported in Fig. 8 and is designated ˜ as G(q). For tethered membrane in a deep flat state, it is an˜ ticipated that G(q) should satisfy a power law scaling with −(2−η) ˜ at small q (large separations).27, 37 The critG(q) ∝ q ical exponent, η, is expected to be 0.849 at low enough temperature according to renormalization group theory, although other estimates ranging from 0.75 to 0.82 has been proposed.26, 27, 29 Clearly, a linear scaling can be observed in Fig. 8. It is interesting to note that the critical exponent η is approximately 0.85 for all the temperatures investigated up to 4000 K, whereas the renormalization group estimate of 0.849 is expected to hold only for membrane far from the crumpling transition. The η of 0.85 is in good agreement with earlier numerical estimate obtained by Los37 with the LCBOPII potential. We note the LCBOPII potential showed different temperature dependence of bending rigidity and TEC. The insensitivity of η with respect to the interaction potential is anticipated for a universal scaling behavior. Fig. 9(a) reports the height-height correlation function H(R) = [h(R) − h(0)]2 at different temperatures. The increase of H(R) slows down toward 620 Å, which is approximately half box size. Fig. 9(b) shows a log-log plot of H(R)
FIG. 8. Log-log plot of the power spectrum of the normal-normal correlation ˜ ˜ function, G(q). At small q, G(q) ∝ q −(2−η) , with a universal scaling constant η = 0.85.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-7
D. Wei and F. Wang
FIG. 9. (a) Height-height correlation function H(R) as a function of distance. (b) Log-log plots of H(R) at the distances from 20.0 to 400.0 Å. A dashed line with a slope of 1.15 is shown for reference.
as a function of distance. A linear region with a slope of 1.15 can be seen when R is no more than approximately 150 Å. The deviation from linear scaling at large R is most likely caused by finite size effects. A slope of 1.15 indicates that ˜ this universal H(R) ∝ R2−η , with η = 0.85. As for G(q), scaling constant holds for graphene up to 4000 K. The exponent is in good agreement with the finding of Los et al.37 obtained with a smaller box at 300 K. A more detailed view is available from Fig. 7(b) since a denser distance grid was used in this work when compared with the Los study.
IV. SUMMARY
We performed MD simulations on large graphene sheets with the realistic and efficient PPBE-G potential recently developed by force matching PAW-PBE at finite temperatures. Simulation with the PPBE-G potential is the best approximation to a direct CPMD simulation with PAW-PBE. Large simulation boxes allow better modeling of long wave-length oscillations that are important for converging TEC and the critical scaling exponent η at low q. To the best of our knowledge, our simulation box is significantly larger than any prior simulations on either TEC or the critical scaling properties
J. Chem. Phys. 141, 144701 (2014)
of graphene. The temperature dependence of η has not been investigated prior to our work. The TEC of graphene remains negative for the entire temperature range investigated consistent with previous studies with the PBE functional at temperatures up to 2000 K. Although the RDF of graphene is still oscillatory at 620 Å, the reduced peak heights of RDF decay as inverse power law to the distance. Such a decay is consistent with the requirement of the Mermin-Wagner instability. Although graphene loses long range translational order as indicated by an inverse power law decay of GT (R), graphene does possess long range orientational order. The PPBE-G potential predicts the renormalized bending rigidity to increase with temperature consistent with earlier simulations with the LCBOPII potential,16, 18 which was contradictory to the REBO-II prediction.14 The normal-normal correlation of graphene G(R) does not show an exponential decay. This behavior makes the standard definition of normal persistence length lose its significance for graphene. Consistent with previous studies, even with much larger simulation boxes, the power spectrum of G(R) decays according inverse power law with a critical exponent of –(2 − η), with η = 0.85. Assuming our potential remains reliable up to the highest temperature simulated, the critical exponent η of 0.85 is a constant for graphene up to 4000 K, indicating graphene to be in the deep flat regime even at less than 1000 K below the experimental melting temperature. It is thus intriguing whether graphene has a crumpling transition before it melts. Graphene is a tethered membrane that loses translational order at large R. By losing its long range translational order, graphene is no longer periodic and does not satisfy the general accepted definition of a crystal in 3D. Although being referred to as a crystal in 2D, graphene does not possess conventional crystalline long-range order. At the same time, graphene shows strong anharmonic coupling and does not follow thermodynamics of a liquid membrane. We note that one limitation of a force-field based MD simulation is the inability to treat electronic excitations that may be important for large graphene sheets. This is a limitation shared also by DFT based ab initio MD simulations, such as CPMD. With AFM based force fields, we can map DFT potential energy surfaces to less computational intensive formulas; however, the accuracy of an AFM force field will be bound by the reference method, in this case Kohn-Sham DFT. On the other hand, the scaling properties of graphene shows universality class behavior, we do not anticipate minor inaccuracy of the PBE functional or explicit treatment of electronic excitation will change the overall scaling. Nonetheless, developing more accurate electronic structure models for graphene and explicit inclusion of phonon-exciton coupling are important future directions of research. ACKNOWLEDGMENTS
This work is supported by key scientific and technological projects of Chongqing, China (cstc2012ggC50002) and NSF CAREER Award No. CHE0748628. The computer resource for the study was provided by the High Performance Center of University of Arkansas and the Super-Computer
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
144701-8
D. Wei and F. Wang
Center of Chongqing Institution of Green and Intelligent Technology, Chinese Academy of Sciences. F.W. acknowledges support from the Arkansas Biosciences Institute. 1 K.
S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306(5696), 666 (2004). 2 J. C. Meyer, A. K. Geim, M. I. Katsnelson, K. S. Novoselov, T. J. Booth, and S. Roth, Nature (London) 446(7131), 60 (2007). 3 F. Guinea, B. Horovitz, and P. Le Doussal, Phys. Rev. B 77(20), 205421 (2008). 4 R. Miranda and A. L. Vazquez de Parga, Nat. Nano 4(9), 549 (2009). 5 A. L. Vázquez de Parga, F. Calleja, B. Borca, M. C. G. Passeggi, J. J. Hinarejos, F. Guinea, and R. Miranda, Phys. Rev. Lett. 100(5), 056807 (2008). 6 W. Z. Bao, F. Miao, Z. Chen, H. Zhang, W. Y. Jang, C. Dames, and C. N. Lau, Nat. Nanotechnol. 4(9), 562 (2009). 7 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81(1), 109 (2009). 8 M. I. Katsnelson and A. K. Geim, Philos. Trans. Royal Soc., A 366(1863), 195 (2008). 9 R. Peierls, Ann. Inst. Henri Poincare 5, 177 (1935). 10 L. Landau, Phys. Z. Sowjetunion 11, 26 (1937). 11 L. Peliti and S. Leibler, Phys. Rev. Lett. 54(15), 1690 (1985). 12 N. D. Mermin, Phys. Rev. 176(1), 250 (1968). 13 D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni, and S. B. Sinnott, J. Phys.: Condens. Matter 14(4), 783 (2002). 14 P. Liu and Y. W. Zhang, Appl. Phys. Lett. 94(23), 231912 (2009). 15 J. H. Los, L. M. Ghiringhelli, E. J. Meijer, and A. Fasolino, Phys. Rev. B 72(21), 214102 (2005). 16 A. Fasolino, J. H. Los, and M. I. Katsnelson, Nat. Mater. 6(11), 858 (2007). 17 P. G. De Gennes and C. Taupin, J. Phys. Chem. 86(13), 2294 (1982). 18 M. I. Katsnelson and A. Fasolino, Acc. Chem. Res. 46(1), 97 (2013). 19 K. N. Kudin, G. E. Scuseria, and B. I. Yakobson, Phys. Rev. B 64(23), 235406 (2001). 20 R. C. Thompson-Flagg, M. J. B. Moura, and M. Marder, Europhys. Lett. 85(4), 46002 (2009).
J. Chem. Phys. 141, 144701 (2014) 21 F.
L. Braghin and N. Hasselmann, Phys. Rev. B 82(3), 035407 (2010). I. Savvatimskiy, Carbon 43(6), 1115 (2005). 23 J. Ma, D. Alfè, A. Michaelides, and E. Wang, Phys. Rev. B 80(3), 033407 (2009). 24 Y. Kantor, M. Kardar, and D. R. Nelson, Phys. Rev. Lett. 57(7), 791 (1986). 25 J. P. Kownacki and D. Mouhanna, Phys. Rev. E 79(4), 040101 (2009). 26 P. Le Doussal and L. Radzihovsky, Phys. Rev. Lett. 69(8), 1209 (1992). 27 Z. Zhang, H. T. Davis, and D. M. Kroll, Phys. Rev. E 48(2), R651 (1993). 28 M. J. Bowick and A. Travesset, Phys. Rep. 344(4–6), 255 (2001). 29 M. J. Bowick, S. M. Catterall, M. Falcioni, G. Thorleifsson, and K. N. Anagnostopoulos, J. Phys. I (France) 6(10), 1321 (1996). 30 D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134(18), 184704 (2011). 31 O. Akin-Ojo, Y. Song, and F. Wang, J. Chem. Phys. 129(6), 064108 (2008). 32 O. Akin-Ojo and F. Wang, J. Comput. Chem. 32, 453 (2011). 33 O. Akin-Ojo and F. Wang, J. Phys. Chem. B 113(5), 1237 (2009). 34 P. E. Blöchl, Phys. Rev. B 50(24), 17953 (1994). 35 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77(18), 3865 (1996). 36 K. V. Zakharchenko, R. Roldán, A. Fasolino, and M. I. Katsnelson, Phys. Rev. B 82(12), 125435 (2010). 37 J. H. Los, M. I. Katsnelson, O. V. Yazyev, K. V. Zakharchenko, and A. Fasolino, Phys. Rev. B 80(12), 121405 (2009). 38 D. Liesegang and C. Oligschleger, J. Mod. Phys. 5, 149 (2014). 39 S. Plimpton, J. Comput. Phys. 117(1), 1 (1995). 40 D. Yoon, Y.-W. Son, and H. Cheong, Nano Lett. 11(8), 3227 (2011). 41 Y. Magnin, G. D. Förster, F. Rabilloud, F. Calvo, A. Zappelli, and C. Bichara, J. Phys.: Condens. Matter 26(18), 185401 (2014). 42 M. Pozzo, D. Alfè, P. Lacovig, P. Hofmann, S. Lizzit, and A. Baraldi, Phys. Rev. Lett. 106(13), 135501 (2011). 43 N. Mounet and N. Marzari, Phys. Rev. B 71(20), 205214 (2005). 44 K. V. Zakharchenko, M. I. Katsnelson, and A. Fasolino, Phys. Rev. Lett. 102(4), 046808 (2009). 45 A. M. Almudallal, S. V. Buldyrev, and I. Saika-Voivod, J. Chem. Phys. 140(14), 144505 (2014). 46 U. Gasser, J. Phys.: Condens. Matter 21(20), 203101 (2009). 47 K. V. Zakharchenko, J. H. Los, M. I. Katsnelson, and A. Fasolino, Phys. Rev. B 81(23), 235439 (2010). 22 A.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 155.33.16.124 On: Sat, 22 Nov 2014 07:19:38
