Binding and interlayer force in the near-contact region of two graphite slabs: experiment and theory
Tim Gould, Ze Liu, Jefferson Zhe Liu, John F. Dobson, Quanshui Zheng, S. Lebègue
aa r X i v : . [ c ond - m a t . m t r l - s c i ] N ov Binding and interlayer force in the near-contact region of two graphite slabs:experiment and theory
Tim Gould, a) Ze Liu, Jefferson Zhe Liu, b) John F. Dobson, Quanshui Zheng,
2, 4 and S. Leb`egue Queensland Micro and Nano Technology Centre, Nathan campus, Griffith University, 170 Kessels Road, Nathan,QLD 4111, Australia Department of Engineering Mechanics and Center for Nano and Micro Mechanics, Tsinghua University,Beijing 100084, China Department of Mechanical and Aerospace Engineering, Monash University, Clayton, VIC 3800,Australia Institute of Advanced Study, Nanchang University, Nanchang, China Laboratoire de Cristallographie, R´esonance Magn´etique et Mod´elisations (CRM2, UMR CNRS 7036) InstitutJean Barriol, Universit´e de Lorraine BP 239, Boulevard des Aiguillettes 54506 Vandoeuvre-l`es-Nancy,France
Via a novel experiment, Liu et al. [Phys. Rev. B, , 205418 (2012)] estimated the graphite binding energy,specifically the cleavage energy, an important physical property of bulk graphite. We re-examine the dataanalysis and note that within the standard Lennard-Jones model employed, there are difficulties in achievinginternal consistency in the reproduction of the graphite elastic properties. By employing similar modelswhich guarantee consistency with the elastic constant, we find a wide range of model dependent bindingenergy values from the same experimental data. We attribute some of the difficulty in the determinationof the binding energy to: i) limited theoretical understanding of the van der Waals dispersion of graphitecleavage, ii) the mis-match between the strong bending stiffness of the graphite-SiO cantilever and the weakasymptotic inter-layer forces that are integrated over to produce the binding energy. We find, however, thatthe data does support determination of a maximum inter-layer force that is relatively model independent. Weconclude that the peak force per unit area is 1 . ± . . ± . I. INTRODUCTION
Graphene has attracted sustained interest in recentyears because of its unusual electronic, magnetic and me-chanical properties . Applications that depend on me-chanical properties include, for example, flexible touch-screens and graphene-coated oscillating sensor devices.These can be based on large-scale high-quality flexiblevapour deposited graphene sheets . To model such ap-plications one requires a reliable knowledge of the forceand binding energy between graphene layers, quantitiesthat have recently been controversial both at the theo-retical and experimental level. For example, results forthe graphene layer binding energy vary by at least a fac-tor of two between different experiments . An equallylarge spread of predictions is found amongst theoreticalanalyses . Any fresh experimental insight on this sys-tem is therefore important.A recent paper from some of the authors used thebulk mechanical properties of graphite to establish avalue for the inter-graphene-layer binding energy indi-rectly from displacement measurements. In the exper-iment, a flake of graphite with a SiO backing was al-lowed to cantilever into a stepped graphite substrate andatomic force microscopy (AFM) used to measure the pro-file of the cantilever. This profile was then fit to the a) Electronic mail: t.gould@griffith.edu.au b) Electronic mail: [email protected] prediction from a finite element analysis using a bind-ing force of cleavage vs. displacement relationship de-rived from a parameterised Lennard-Jones (LJ) poten-tial. One parameter (the effective binding energy) wasallowed to vary in the LJ potential and a best fit wasfound to the measured data predicting a binding energyof 0 . ± . .In this paper we shall first discuss certain problems inthe previous theoretical analysis of experiment . Themost significant of these is an internal inconsistency be-tween the effective elastic coefficient used in Hooke’slaw for the contact layer (20.0GPa), and other layers(36.5GPa), discussed in depth in Section II B. We willshow that this discrepancy comes from limitations of thepopular Lennard-Jones model when applied to layeredmaterials, and that layered materials present a particu-larly difficult case for indirect measurements of energy.This difficulty is due to a combination of poorly under-stood theory, and variable length and force scales for dif-ferent energetic contributions.We will then reanalyse the experimental data fromRef. 19 using improved binding energy models to bet-ter estimate the energetics of cleavage. These improvedmodels remove the discrepancies in the elastic coefficientby construction, and should thus be considered more re-liable in the near-contact regime. We will show that thepredicted binding energy varies greatly between the dif-ferent models and should thus be considered unreliable.However, the force vs distance curves for intermediateinterlayer distances show significantly less variation. Un-like the energy, the force is not affected by difficulties(discussed in Section II C) arising from the mismatch be-tween the large bending stiffness of the cantilever on onehand, and the weakness of the attractive interlayer vander Waals forces on the other hand.We will finally conclude that the previous theoreticalanalysis on experimental results is likely to substan-tially underestimate the cleavage energy value, but canbe used to make reliable predictions of the “peak force”(the minimum force required to cleave the layers). Thepeak force per unit area is found to be 1 . ± . ± II. THEORY AND ANALYSIS
As mentioned in the Introduction, we identified a num-ber of problems with the previous theoretical analysis ofexperimental results . These problems highlight the dif-ficulties of measuring the energetic properties of graphite,a material where direct experiments are difficult and the-ory is not always comprehensive or conclusive. We out-line the identified issues below, with the aim of account-ing for each in a reanalysis of the same data. A. Cleavage vs binding
We first note that Liu et al. analyzed the energy ofcleavage (the energy required to split a bulk along aplane) rather than the bulk-layer binding energy as nor-mally quoted (the energy required to split a bulk intowidely separated atomic planes). A discussion of thedifference between these energies is found in the supple-mentary material of Ref. 20, and from a more formal per-spective in Ref. 21. While these quantities are usually ofthe same order of magnitude, they will differ wheneverthe interaction between second-neighbor planes of atomsis significant, which it typically is in a van der Waalsbonded system such as graphene. We will subsequentlyuse E Clv to refer to the cleavage energy and E BLB to referto the bulk-layer binding energy.While high-level ab initio data is not available forcleavage of graphite, the binding energy of graphite wasfound by Spanu et al. with Quantum Monte-Carlo(QMC) method to be 56meV/Atom=0 . , and byLeb`egue et al. under the Random Phase Approxima-tion (RPA) to be 48meV/atom=0 . . Thrower etal. used lower-level theory benchmarked against exper-imental results for adsorption of polycyclic aromatic hy-drocarbons on graphite to predict a graphite binding en-ergy of 57 ± . ± . . The differ-ence between the cleavage and binding energies has alsobeen found on a similar system , and via lower-level the-ory adjusted to match known high-level theory .Podeszwa extrapolated from anthracene, pyrene andcoronene to estimate the exfoliation, binding and cleav-age energies of graphite. While his model neglects long- ranged (in the plane) plasmon interaction effects due tothe finite size of the coronenes, it is expected to yield agood estimate of the relative strengths near contact. Hefinds a cleavage energy 14% greater than the binding en-ergy of graphene based on a four sheet model. This givesan absolute cleavage energy of E Clv = 48 . . . As Podeszwa notes, this extrapolation ne-glects certain bulk binding properties of graphite, andmay underestimate the binding and cleavage energies.In a recent publication , Bj¨orkman et al. investi-gated bulk properties in various two dimensional ma-terials using various theories. For graphite they founda 2 . / Ang = 0 . difference between thegraphite exfoliation and graphene binding energies. Un-published results suggest that there is a further 10%-15% increase in energy for cleavage compared to exfolia-tion (an absolute value of 0 . − . based on RPAand QMC results). This is of similar relative magnitudeto the energy difference found by Podeszwa , and theenergy difference predicted by LJ pair summation mod-els using experimental layer spacing.These various theories thus predict a cleavage energy E Clv between 0 . and 0 . , with a bindingenergy E BLB between 0 . and 0 . with thelower bound likely underestimated. The cleavage energy E Clv = 0 . ± . estimated via a theoretical anal-ysis of the cantilever experiment would correspond to abinding energy of E BLB = 0 . ± . , around halfthat of the higher-level theories for bulks, and at least30% less than the lowest reasonable theoretical predic-tion. B. Near-contact force
Putting aside the distinction between binding andcleavage energies, and taking the model in the previouswork prima facie we note another difficulty that war-rants further analysis. The potential model employed isderived directly from the two-parameter LJ potential ofRef. 24. It gives the energy Φ LJ as a function of the dis-tance parameter ν describing the stretching of a singleinterplanar distance during the cleavage process, whileleaving other interplanar distances unchanged. For thispaper we rewrite the energy in terms of x ≡ D − D ≡ ν − ν , the difference of the distance D between the surfaces aftercleavage, and the graphite equilibrium interlayer spacing D such that ddx Φ LJ (0) = 0 ( ν is defined as the value of ν that minimises Φ LJ ).Equation A6 of Ref. 19’s supplementary material gives,after some algebra,Φ LJ ( x ) := ∞ X n =1 n Φ ( x + ν + nσ ) , (1)where Φ ( δ ) = αE Clv σ (cid:20)
25 ( σ/δ ) − ( σ/δ ) (cid:21) (2)is the interlayer potential between two layers and issummed over all possible layer interactions to form Φ LJ .Here E Clv is the cleavage energy (called the “bindingenergy” in Ref. 19) while α and σ are constants, thelatter being related to the interlayer spacing and tak-ing the value σ = 0 . x , the energy takes its minimum valueat x = 0 which gives ν = − . D =0 . α = 10 . -2 is chosen to make Φ LJ (0) = − E Clv andcan be derived from A2-A6 of the supplementary mate-rial for Ref. 19.The potential Φ LJ can be considered a non-linearextension to Hooke’s law between planes, reducing toHooke’s law for small displacements. When two graphitesurfaces are parallel and near equilibrium lattice spacing,the resulting force should be proportional to the displace-ment, and should correctly reproduce the C coefficientof bulk graphite measured in previous experiments .In the model (1), small displacements occur when x ≡ D − D ≈ LJ ( x ) ≈ − E Clv + ˜ C D x , F ( x ) ≈ ˜ C xD (3)where F ( x ) is the cleavage force near ‘contact’. Here˜ C = D d Φ LJ ( x ) dx (cid:12)(cid:12) x =0 (4)is, in fact the elastic coefficient via cleavage, and cannottrivially be compared with existing experiments on bulkgraphite. As discussed in Appendix A we can correct (4)to obtain the true C coefficient for stretching of bulkgraphite via C = ˜ C + ∆ C , (5)where ∆ C is a small correction.It is thus clear that we can use Φ LJ ( x ) to evaluate theinterlayer elastic coefficient C via (4) and (5). Using E Clv = 0 . and D = 0 . C =19 . C = 20 . which givevalues between 36.5GPa and 40.7GPa.This causes a discrepancy in the overall physical andmechanical model in Ref. 19. For all layers except theopen surface, an interplane elastic coefficient of C =36 . C is proportional to E Clv for thepotential Φ LJ ( x ) employed, and gives an elastic constant C = 20 . LJ ( x ) defined via (1)]. The experimental elastic con-stant can be matched by setting the cleavage energy to0 . ± . , however this is inconsistent with thevalue 0 . deduced by other means in Ref. 19. FIG. 1. AFM measurements of one of the cantilevers showingthe effect of in-plane stiffness on the profile. Dotted lines areprovided as a visual guide. The experimental setup is shownbelow, showing the scanning direction in (e). (e) ∆ ∆ =53.7±0.9nm ∗∆ =53.8±2.7nm5 µ mMM3A tip(a)(b)(c)(d) (f) h=155nm h C. Intermediate and long distance force
In the experimental setup reported in Ref. 19, graphitewas arranged with its layers lying parallel to the base ofthe cantilevered segment. The bending stiffness of thegraphite cantilever is proportional to the in-plane elas-tic coefficient of graphite (which is one of the largestknown at ∼ ∼ intermediate part of the binding force willbe emphasised in any analysis of the experimental data,while the long distance binding will be harder to deter-mine. This makes evaluation of the binding energy diffi-cult as it depends on the integral of the force from infiniteseparation to finite separation. This also means that thebest fit [using (1)] to a given experiment must matchboth the difficult to measure longer-ranged forces, andthe near and intermediate forces. With only one free pa-rameter matching all three ranges is very difficult, andthe intermediate region will likely dominate the fit.The experimental difficulties are further amplified bythe fact that theoretical understanding of the longer-ranged van der Waals potential of graphite cleavage islimited. Here high-level theory results are valid onlyfor the extreme outer limit where corrections for finiteDirac cones and other non-asymptotic physics are notrequired.This lack of theoretical insight comes from difficul-ties in understanding the contribution of the π z and sp orbitals on the dispersion processes involved in cleav-age. While these dispersion forces contribute a signficantamount to the intermediate region energy profile ( ∼ ), they contribute comparitively less to theintermediate force. Thus a failure to properly account fordispersion is less problematic for interpretation of peakforce measurements compared to the cleavage energy.Furthermore, the dispersion contribution is itself splitinto two components, which dominate at different lengthscales. Here the sp orbitals dominate the dispersion inthe intermediate region, while π z orbitals dominate in theasymptotic region. The different force and length scalesinvolved in dispersion thus make extrapolation from thenear-contact region to the asymptotic region very diffi-cult. This makes mathematical fitting of an all-regionforce unreliable without firm theoretical backing, withconsequences for the energy. It does not, however, af-fect measurement of near-contact forces, which dependminimally on the unknown asymptotics. III. ALTERNATE MODELS
As mentioned in Section II B, the Lennard-Jones model(1) with cleavage energy determined by best fit to ex-periment failed to reproduce the known near-contactforce of graphite. Ideally any model potential shouldreproduce the correct near-contact force by construc-tion, via the experimentally measured interlayer dis-tance D = 0 . C = 36 . as the van der Waals power law − C D − well away from contact. A model potential can then beformed that allows the binding energy to vary as a param-eter, while matching these known properties of graphite.The LJ potential, with just two parameters, is not ableto simultaneously fit D , C and E Clv , and leads to C ∝ E Clv . We thus propose two new models of thebinding energy that decouple the elastic coefficient fromthe cleavage energy, and thus allow better reproduction of known properties of graphite.Our first model is designed to reproduce the RPA en-ergy curve in the near-contact region when the input en-ergy takes the RPA value E Clv = 0 . . Here, inaddition to matching the binding distance and elasticcoefficient, we match the non-linear coefficient C = − x/D the force obeys F ( x ) ≈ C ( x/D ) + C ( x/D ) . Onesuch “near-contact model” (NCM) isΦ NCM ( x ) = − E Clv c x + c x e − kx (6)where x = D − D is the deviation from the experi-mental lattice spacing. Here c = C / (2 D E Clv ) and c = C / (6 D E Clv ) ensure that the first three deriva-tives of Φ
NCM are equal to their RPA values at contact D = D for arbitrary E Clv . Free parameter k >
NCM ( x → ∞ ) ∝ x − . We choose k = 8 toensure that the force dΦ NCM d x is single peaked for x > E Clv = 0 . is used.Alternatively, one may wish to keep the LJ-like form ofthe model potential, while reproducing the known physi-cal properties D and C as well as the correct − C /D asymptotic form (if not its coefficient C ). This suggestsa model potential of general formΦ LJ(p) ( x ) = E Clv p − (cid:20) x/D ) p − p (1 + x/D ) (cid:21) (7)where the exponent p = C D / (2 E Clv ) is chosen to en-sure the experimental C elastic coefficient can be re-produced. Such an approach is perhaps physically lessjustified than the above in the near contact region, butdoes not require an RPA third-order elastic coefficient asadditional input, and may be more accurate away fromcontact. IV. RESULTS
Following the previous work , we use the three differ-ent model potentials Φ LJ , Φ NCM and Φ
LJ(p) [from re-spectively equations (1), (6) and (7)] in FEA modelsto determine a model-dependent cleavage energy E Clv .This was defined to be the value (for each model) thatcaused the elastically deformed curves determined viaFEA to best match the experimentally measured pro-files. The variation of the predicted binding energy wasvery great even within a given example geometry, rang-ing from E Clv = 0 . for the near-contact modelΦ NCM , to E Clv = 0 . for the LJ model Φ and E Clv = 0 . for the modified LJ model Φ LJ(p) . Themodified LJ model best matches previous theoretical pre-dictions of the cleavage energy, but the error bar acrossthe three models is too signficant for any result to be con-sidered reliable. This is unsurprising as the cleavage en-ergy depends partly on the very weak dispersive forces at
FIG. 2. Cleavage force per unit area under the three differentmodels plotted against the inter-layer distance D (where theequilibrium interlayer spacing is D = 0 . bulk layer binding force is included for illustrative purposes,although it is not in general the same as the cleavage force. F ( G P a ) D (nm)Forces with best-fit energies LJLJ(p)ModelRPA
TABLE I. Best fit cleavage energy, and corresponding peakforce and position for different models. RPA results are pro-vided for illustration, not comparison.Model E Clv (J/m ) F m (GPa) D m (nm)LJ Φ 0.19 0.96 0.388NCM Φ NCM p ) Φ LJ(p) a a Using data from Ref. 16 and additional data points. large distance, which are difficult to determine accuratelyvia the current experimental setup . Additionally, theasymptotic van der Waals behaviour for cleavage is un-known, and none of the models can be guaranteed to holdtrue in the more distant limit.From the potential it is trivial to determine the force F = dΦd x , and in the near-contact region this should be lessmodel-dependent and given more reliably by experiment.Indeed by construction both Φ NCM and Φ
LJ(p) will pre-dict identical linearly varying forces near the inter-layerbinding distance. Using the best fit E Clv parameters foreach model yields distance-dependent forces displayed inFigure 2. It is immediately apparent that while the forcein the tail varies greatly between the models, the posi-tion and magnitude of the “peak” force (the maximum)varies much less signficantly. The “peak force” is a physi-cally important property of graphite, giving the force perunit area require to “cleave” graphite in two at a singlesurface. This agrees with the arguments given in Sec-tion II C that the analysis of experiment of Ref. 19 willbe most sensitive to the force in the intermediate range.Using the best fit E Clv for each of the three models,we semi-analytically determine the peak force F m and the inter-layer distance D m at which it occurs. We tabu-late the results in Table I, including the peak bulk layerbinding (as opposed to cleavage) force predicted by theRPA for illustrative comparison. The peak force variesby up to 27% over all three models, which is very goodagreement for such different models. For the two modelsΦ NCM and Φ
LJ(p) with the experimental elastic coefficientof bulk graphite included as a parameter the variation iseven smaller still, with only a 5pm difference to the po-sition and an 8% difference in the peak force. This smallvariation is a very positive feature, as the cleavage energyparameter varies by over 100%.The peak force predicted by the RPA is significantlylarger than those found via the analysis of the experi-ment with any of the models. While at first glance thisis worrying, we note that the RPA data is for the differentproblem of bulk layer binding (under uniaxial stretchingperpendicular to the planes) and cannot be directly com-pared. Indeed one expects the force laws for cleavage,exfoliation and bulk layer binding to show greater varia-tion than the ‘binding’ energies as each obeys a differentasymptotic power law , with bulk Coulomb screeningreducing the dispersion force in cleavage and exfoliation.This screening effect suggests that the peak force of cleav-age should be less than that of bulk layer binding, andthus one should expect an apparent discrepancy betweeni) the values of peak forces deduced from the experimentsusing the three models, and ii) theoretical predictionsfrom the RPA. Full resolution of the differences betweenpeak force in bulk-layer binding and cleavage would re-quire further experiments and/or further RPA calcula-tions (or other high-level theory) for the cleavage geom-etry. These are beyond current capabilities, however. V. CONCLUSION
Overall we conclude that the binding energy evaluatedthrough the theoretical analysis of experimental resultsin the previous work (Ref. 19) is likely to be substantiallyunderestimated. Firstly, the model employed in Ref. 19yields a near-contact force law with an elastic coefficientof C = 20 . with C = 38 . ± . C ∝ E Clv suggesting that E Clv ≈ . would be required to match the ex-periment elastic parameter. Secondly, the experimentalsetup and model potential [see our (1) or (A6) of the sup-plementary material of Ref. 19] are designed for the mea-surement of cleavage energies E Clv , not the bulk-layerbinding energies E BLB with which their results were com-pared. It is believed that E BLB ≈ . E Clv in graphiticsystems, leading to a value E BLB ≈ . from theexperimental analysis, approximately half the value esti-mated by high-level theory .By attempting new fits to the same data, using modelpotentials which include the experimental C coefficientas an input, we find that the predicted cleavage en-ergy E Clv varies dramatically (by 100%) depending onthe model potential used. However, the size and posi-tion of the peak force deduced from the experiment (theminimum force per unit area required to fully “cleave”graphite at a surface) is much less dependent on themodel, showing only a 27% variation in magnitude acrossall three models, reduced to 8% between the two mod-els that reproduce the experimental C coefficient. Wethus conclude that the peak force per unit area can beaccurately predicted by this experiment, and that it is1 . ± . . ± . π z dispersion.We thus recommend against the use of the Lennard-Jonesmodels for layered systems, and recommend instead theuse of something like the LJ( p ) model (7) when mod-elling interlayer forces in systems with known elastic co-efficients. ACKNOWLEDGMENTS
We would like to thank Torbj¨orn Bj¨orkman for helpfuldiscussion on cleavage, exfoliation and binding in lay-ered solids. T.G. and J.F.D. were supported by ARCDiscovery Grant DP1096240. S.L. acknowledges finan-cial support from the Universit´e de Lorraine throughthe program “Soutien `a la dimension internationale dela recherche”. Q.S.Z. acknowledges the financial sup-port from NSFC through Grants No. 10832005 andNo. 10972113, the 973 Program through Grants No.2007CB936803 and No. 2013CB934200. J.Z.L. acknowl-edges Seed Grant 2013 from the engineering faculty ofMonash University.
Appendix A: Elastic coefficients
In general the C elastic coefficient is the ratio of inter-layer force to interlayer displacement induced by stretch-ing all layers evenly (with appropriate scaling). By con-trast the model used in Ref. 19 involves displacementbetween two particular layers (cleavage) with other layerspacing kept essentially fixed, a process associated withelastic constant ˜ C . Employing the interlayer potentialΦ [equation (2)] we can find the coefficient in the twodifferent models via [using x = ν − ν in equation (4) for˜ C and a similar derivation for C ]˜ C = D d dν ∞ X n =1 n Φ ( nσ + ν ) (cid:12)(cid:12) ν = ν , (A1) C = D d dν ∞ X n =1 Φ ( nσ + nν ) (cid:12)(cid:12) ν = ν , (A2)where (A2) should reproduce the true C coefficient ofgraphite for a good interlayer potential model. The dif-ference between these two coefficients is∆ C ≡ C − ˜ C = D ∞ X n =2 d dν [Φ ( nσ + nν ) − n Φ ( nσ + ν )] ν = ν (A3)since the n = 1 terms cancel. We thus expect ∆ C tobe small compared to C as Φ (2 σ ) ≪ Φ ( σ ) and | ν | ≪ σ . Inserting the experimental parameters ν /σ = − . D = 0 . E BLB = 0 . gives ∆ C =0 . C = 19 . . K. Novoselov, A. Geim, S. Morozov, D. Jiang, Y. Zhang,S. Dubonos, I. Grigorieva, and A. Firsov, Science , 666(2004). A. K. Geim and K. S. Novoselov, Nature Materials , 183 (2007). M. I. Katsnelson, Materials Today , 20 (2006). A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov,and A. K. Geim, Rev. Mod. Phys. , 109 (2009). M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, NaturePhysics , 620 (2006). L. Hu and G. Gruner, U.S. Patent No. 8390589 B2 (5 Mar 2013). K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim,J. H. Ahn, P. Kim, J. Choi, and B. H. Hong, Nature , 706(2009). L. A. Girifalco and R. A. Lad, The Journal of Chemical Physics , 693 (1956). L. X. Benedict, N. G. Chopra, M. L. Cohen, A. Zettl, S. G. Louie,and V. H. Crespi, Chem. Phys. Letters , 490 (1998). R. Zacharia, H. Ulbricht, and T. Hertel, Phys. Rev. B , 155406(2004). M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, and B. I.Lundqvist, Phys. Rev. Lett. , 246401 (2004). S. D. Chakarova-K¨ack, E. Schr¨oder, B. I. Lundqvist, and D. C.Langreth, Phys. Rev. Lett. , 146107 (2006). E. Ziambaras, J. Kleis, E. Schr¨oder, and P. Hyldgaard, Phys.Rev. B , 155425 (2007). M. Hasegawa and K. Nishidate, Phys. Rev. B , 205431 (2004). L. Spanu, S. Sorella, and G. Galli, Phys. Rev. Lett. , 196401(2009). S. Leb`egue, J. Harl, T. Gould, J. G. ´Angy´an, G. Kresse, andJ. F. Dobson, Phys. Rev. Lett. , 196401 (2010). J. D. Thrower, E. E. Friis, A. L. Skov, L. Nilsson, M. Andersen,L. Ferrighi, B. Jrgensen, S. Baouche, R. Balog, B. Hammer, andL. Hornekr, The Journal of Physical Chemistry C , 13520(2013). T. c. v. Buˇcko, S. Leb`egue, J. Hafner, and J. G. ´Angy´an, Phys.Rev. B , 064110 (2013). Z. Liu, J. Z. Liu, Y. Cheng, Z. Li, L. Wang, and Q. Zheng, Phys.Rev. B , 205418 (2012). T. Bj¨orkman, A. Gulans, A. V. Krasheninnikov, and R. M.Nieminen, Phys. Rev. Lett. , 235502 (2012). T. Gould, K. Simpkins, and J. F. Dobson, Phys. Rev. B ,165134 (2008). R. Podeszwa, The Journal of Chemical Physics , 044704(2010). T. Bj¨orkman and A. Gulans, (private communication). L. A. Girifalco and M. Hodak, Phys. Rev. B , 125404 (2002). O. L. Blakslee, D. G. Proctor, E. J. Seldin, G. B. Spence, andT. Weng, Journal of Applied Physics , 3373 (1970). W. B. Gauster and I. J. Fritz, Journal of Applied Physics ,3309 (1974). N. Wada, R. Clarke, and S. Solin, Solid State Communications , 675 (1980). A. Bosak, M. Krisch, M. Mohr, J. Maultzsch, and C. Thomsen,Phys. Rev. B , 153408 (2007). T. Gould, E. Gray, and J. F. Dobson, Phys. Rev. B , 113402(2009). T. Gould, J. F. Dobson, and S. Leb`egue, Phys. Rev. B ,165422 (2013). T. Gould, S. Leb`egue, and J. F. Dobson, Journal of Physics:Condensed Matter25