The Challenge of Inter-planar Bonding in Graphite

In Section6.1 and AppendixA, we are confronted with the challenge of determining the equilibrium crystal structure of graphite. The difficulty originates from the strongly anisotropic bonding in graphite. While the bond length of the strong covalent in-plane bonds are well described for a variaty of exchange-correlation functionals (see Tab. 6.1 and 6.2), the interaction between the layers is very weak and is generally attributed to vdWforces. Closely related to the question of the interlayer bond distance is the interlayer energy. Unfortunately, a computational investigation within DFTof the interlayer bonding in graphite is a challenge. The problem arises from the local and semi-local approximation to the exchange-correlation functional commonly used inDFTand the intrinsically long-range nature of vdW-bonding. For an accurate description of the interlayer bonding in graphite long range electron correlation effects have to be included in DFT calculations (Ch.2.4). The inclusion ofvdW effects is a very active scientific area of its own (e.g. see [271, 116, 182, 117, 188, 326, 325, 324, 5] and references therein). Björkman et al. [28] and Grazianoet al. [115]

gave a detailed discussion of incorporatingvdW-effects in the theoretical description of layered materials.

In this section, we focus on the comparison of two differentvdWcorrection schemes: the well established Tkatchenko-Scheffler (TS) scheme [326] and the latest progress incorporating many-body effects [325,324,5]. BothvdW schemes depend on the electronic density. Recapping briefly the summary of thevdWsection (Sec.2.4), the TS approach is based on a non-empirical method that includs a sum over pairwise interactions. The many body approach includes long-range many-body dispersion (MBD) effects employ-ing a range-separated (rs) self-consistent screenemploy-ing (SCS) of polarisabilities [324,5] (MBD@rsSCS). It captures non-additive contributions originating from simultaneous dipole fluctuations at different atomic sites, leading to an improved overall screening and accounts for long-range anisotropic effects.

However, the determination of the graphite interlayer binding energy is not just a challenge for theory. To our knowledge, four different experie-ments have been published reporting binding energies in the order of 30-55 meV/atom [106,17,350,201]. In the literature, three different ener-gies are used to describe the interlayer binding energy: exfoliation energy, binding energyandcleavage energy. The differences between these three terms are subtle. Björkmanet al.[29] gave an illustative description of this three

Figure 6.3.: The interlayer binding energy sketch for three different types of binding:a.) exfoliation energy, b.) binding energy, c.) cleavage energy. In (a) and (b), the layer of interest is marked by a box. The scissor and dashed line for the exfoliation (a) and cleavage (c) cases illustrates how the crystal is cut.

energies in the supplemental material of Ref. [29]. We here will follow their line of argumentation. For completeness we will briefly derive the three energies.

Figure 6.3 shows a sketch of the different terms for the binding energy.

The interaction energies between two layers are labelled ase₁ for nearest neighbours,e2for second nearest neighboursetc. The exfoliation energy (E_XF) is the energy to remove a single carbon layer (marked by a box in Fig.6.3(a)) from the graphite bulk structure. We assume that energy con-tributions originating from in-plane relaxation are neglectible, which is a good assumption forABstacked graphite¹. From Fig.6.3(a) we see, that EXF is a sum of the interaction energies

EXF =e₁+e2+. . . =

∑

en. (6.1)

We determine EXF by calculating the total energy of a graphite slab with n layers (En), of a graphite slab with n−1 layers (E_n−1) and an isolated graphene layerE₁. EXFthen is given as

EXF(n) = ¹ Natom

(En −(E_n−1+E₁)) (6.2) whereEstands for calculatedDFTtotal energies andNatomis the number of atoms in the graphite slab.

1The assumption does not necessarily hold for graphite layers rotated by an arbitray angle, because due to the rotation there might be regions close toAAand other close to ABstacking.AAandABstacked graphite have different interlayer distances [220]. This difference results in a layer corrugation and consequently in in-plane relaxation, which leads to sizeable energy contributions.

Another relevant energy is the energy needed to separate all layers from each other by increasing the interlayer distance – the binding energy (EB) (Fig.6.3(b)). When we stretch the crystal to separate the individual layers from each other we break two bonds, but since the bonds are shared the sum of interaction energies (en) is halved

EB = ¹ 2

∑

2·en. (6.3)

A comparison of Eq.6.1and Eq. 6.3shows that the exfoliation energy is equal to the binding energy per layer. We calculate the binding energy of an infinite graphite crystal separated into isolated graphene layers

EB = ¹

Natom E_graphite−N_layerE_graphene

(6.4) where Egraphite (Egraphene) is the total energy of graphite (graphene) and N_layeris the number of layers. In the case of an infinite large slab, Eq.6.2 gives the same energy as Eq. 6.4. In practice, Eq.6.4is used to calculate interlayer binding energies.

The cleavage energy is the energy required to split a graphite crystal in two halves (see Fig.6.3(c)) or, in other terms, it is the energy to create two units of surface. The cleavage energy expressed in terms ofen is

E_cleav =_e1+2·_e2+· · ·+n·_en =

∑

n·_en =2·E_{sur f} (6.5)

Comparing Eq.6.3and Eq.6.5shows, that Eq.6.5is greater than the binding energy and can only serve as an upper bound of the binding energy. The definitions of the exfoliation energy, binding energy and cleavage energy will be helpful when interpreting experimentally obtained graphene inter-planar binding energies.

Figure6.4.:Historyofthegraphiteinterlayerbindingenergy(meV/atom)sincethefirstexperimentin1954[106].Theexperimentalvalues [106,17,350,201,113]aregivenascrossesandthetheoreticalvaluesasfilledcircles.Thecolourofeachcircledependsontheexperimental workcitedinthepublication.GirifalcoandLad[106]gavethebindingenergyinunitsofergs/cm2,theirvalueof260ergs/cm2waserroneously convertedto23meV/atom(showninlightgreen)insteadofthecorrect42meV/atom(darkgreen).Afulllistofreferencesandbindingenergiesis giveninAppendixJ

Finding an accurate description of the graphite interlayer binding energy is an ongoing quest in experiment and theory alike. We first summarise the history of the search for the graphite interlayer binding energy (Fig. 6.4).

The four experiments [106,17,350,201] and the improved interpretation of the 2012 experiment [113] are shown as crosses. Girifalco and Lad [106]

gave their measured binding energy in units of ergs/cm², the value of 260 ergs/cm² was erroneously converted to 23 meV/atom (included in Fig.6.4as a dashed line) instead of the correct 42 meV/atom. In the literature both values were used reference. The different theoretical values are shown as filled circles. The colour of each circle depends on the experimental work cited. A full list of theoretical and experimental values of the graphite binding energy is given in AppendixJ.

Figure 6.5.: Binding energy E_B of graphite as a function of the interplanar distance d usingLDA,PBE,HSE06and two differentvdWcorrection schemes (TS [326] and MBD [325,324,5]). The symbols markE_Bat the equilibrium interlayer distance. The vertical line marks the experimental lattice spacing ofc₀ =6.694Å [226]. For comparisonE_B values from RPA@PBE calculations ( (1) Lebègue et al. [185], (2) Olsenet al. [239]) and from experiment ( (3) Girifalcoet al.[106], (4) Benedictet al.[17], (5) Zachariaet al.[350], (6) and (7) Liuet al.[201,113]) are included.

We first discuss the experimental data included in Fig.6.4. In 1954,Girifalco and Ladprovided the first experimental valueE_XF = (42.5±5)meV/atom for the graphite exfoliation energy obtained indirectly fromheat of wetting measurements [106]. In a unique experimentBenedictet al.determined the binding energy of twisted graphite based on transmission electron micro-scopy (TEM) measurements of the radial deformation of collapsed multiwall

carbon nanotubes. Carbon nanotubes with large radii collapse and lose their circular cross-section, because of an interplay between the intersheet attraction and the resistance to elastical deformation of the curvature of the carbon nanotube (mean curvature modulus). TheTEMdata was evaluated by continuum elasticity theory and a Lennard-Jones (LJ) potential for the description of the graphite sheet attraction. They found a binding energy ofEB = (35⁺₋¹⁵₁₀)meV/atom [17]. Zachariaet al.investigated the graphite interlayer bonding by thermal desorption of polyaromatic hydrocarbons, including benzene, naphthalene, coronene, and ovalene, from a basal plane of graphite. To determine the interlayer binding energy of graphiteZacharia et al.used the thermal desorption data. They assumed the interaction of carbon and hydrogen atoms to be pairwise additive to the total binding energy of the molecules. For the graphite cleavage energy they found Ecleav = (_61.65±₅) meV/atom and EB = (_52.65±₅) meV/atom for the binding energy [350]. Liuet al.performed the first direct measurement of the graphite cleavage energy based on displacement measurements facili-ating mechanical properties of graphite [353]. In the experiment a small graphite flake (a few micrometer) was sheared from a graphite island. The profile of the cantilever graphite flake was measured using atomic force microscopy (AFM). They derived a relationship between the cleaving force and the displacement from a parameterisedLJpotential. They estimated the binding energy to be 85% of the cleavage energy and obtained a value ofEB = (31±₂)meV/atom. In a recent workGould et al.suggested that the binding energy in the latter experiment might be substantially underes-timated, because the experimental data were analysed using aLJpotential.

At large separation theLJpotential gives qualitatively incorrect interlayer binding energies [75]. The other experimental results in the literature might also be affected by the same difficulty!

From a theoretical point of view the description of the interlayer bind-ing energy is equally challengbind-ing. Various methods were employed, ran-ging from semi-empirical approaches to advanced first-principles calcu-lations. The obtained values for the graphite interlayer binding energy scatter widely (see Fig.6.4). In Sec.2.4, we already pointed out that stand-ard local and semi-local approximations such as the LDA and the PBE approximation to the exchange-correlation functional cannot capture the inherently long range nature ofvdWinteractions. A common technique to include the long-range tail of dispersion interactions is to add ener-gies obtained using pairwise interatomic potential toDFT total energies [116,117,128,326, 118]. This approach led to binding energies between 43-82 meV/atom [119,128, 7, 114,200, 38, 42,45]. Another approach is the seamless inclusion ofvdWinteraction in the exchange-correlation func-tional [271,72,182,53,322,207,272] resulting in binding energies between

24-66 meV/atom [271, 167, 59, 207, 115, 272,323, 22]. Spanu et al. [299]

applied a quantum Monte Carlo (QMC) method (variational Monte Carlo and lattice regularized diffusion Monte Carlo [51]) to determine the layer energetics (E_B = (56±5)meV/atom) and equilibrium lattice parameter (d =3.426 Å) of graphite. The calculations were performed in very small simulation supercells (2×₂×2) and finite-size effects may limit the accur-acy of their results. Two groups appliedDFT(PBE) based random phase approximation (RPA) calculations of the correlation.RPA@PBE calculations provide a more sophisticated method for treatingvdWinteractions. Both groups found an equilibrium interplanar distance d = 3.34 Å. However their binding energies differ by 14 meV/atom (EB = 48 meV/atom [185]

andEB =62 meV/atom [239]).Lebègueet al.confirmed theEB(d)∼c3/d³ behavior of the correlation energy at large distances as predicted by ana-lytic theory for the attraction between planar π-conjugated systems [75].

Olsen and Thygesen [239] also calculated the binding energy with the same k-point sampling as Lebègue et al. [185] and found a binding energy of 47 meV/atom. The remaining difference might originate from a different setup of the projector augmented wave method. Using a very dense k-point sampling of (26×26×8) to obtain the Hartree-Fock wave function they found 62 meV/atom for the binding energy.

We calculated the graphite interlayer binding energy as a function of the interplanar distanced=c0/2 for different exchange-correlation function-als (LDA,PBE,HSE06) and vdW correction schemes (TS [326] and MBD [325,324,5]) shown in Fig.6.5. LDAis a purely local approximation to the exchange-correlation functional. Therefore it is not capable of capturing the long-range dispersion interactions. However, it gives a surprisingly good equilibrium lattice parameter of 6.643 Å (Tab.6.1), but it clearly un-derestimates the binding energy in comparison to available experimental and RPA@PBE values [275,54,318,127,185,239]. The PBEcurve shows a very weak bonding (c₀ =8.811 Å). More interesting is the performance of thevdWcorrected exchange-correlation functionals. We first compare the two different vdW correction schemes. For the pairwise TS scheme we find an equilibrium lattice parameter of c₀ = 6.669 Å for PBE+vdW and c0 = 6.641 Å for HSE06+vdW in very good agreement with the ex-perimental value of 6.694 Å [226]. The binding energies are overestimated with 84.0 meV/atom (87.1 meV/atom) for PBE+vdW (HSE06+vdW), re-spectively. Combining theMBD@rsSCSscheme with exchange-correlation functionals (PBEandHSE06) should improve the overall description of the graphite inter-planar bonding, because it additionally captures non-additive contributions and accounts for long-range anisotropic effects. Indeed, for the binding energy we find 47.6 meV/atom and 52.5 meV/atom forPBE in-cluding van-der-Waals effects with full many-body treatment (MBD@rsSCS)

[325,324,5] (PBE+MBD) andHSE06+MBD, respectively. In Fig.6.5in partic-ular thePBE+MBDbut also theHSE06+MBDcurve show a good agreement with theRPAdata of Lebègueet al.[185] close to the minimum.

As discussed above, the interpretation of the experimental data is usually based on theLJpotential, which might lead to an underestimation of the binding energy in graphite [113]. Here, high quality theoretical values of the graphite binding energy are a challenge because electronic structure calculations are difficult to converge [239]. In the case of graphite interlayer binding energies reliable benchmark data is still missing. However, a recent high qualityQMCstudy provides a benchmark value for bilayer graphene [220]. In the work of Mostaaniet al.[220], finite size effects were controlled by using large supercells (6×6). They obtained a binding energy forAB stacked bilayer graphene ofEB=17.7 meV/atom and an inter-planar lattice constant of 3.384 Å [220]. In comparison, we find binding energies of bilayer graphene ofE_B=22.0 meV/atom usingPBE+MBDand ofE_B=36.5 using PBE+vdW. For bilayer graphene, we find good agreement for the binding energy betweenPBE+MBDand high level benchmark data (QMC) and an overestimation forPBE+vdW.

TheMBD@rsSCSmethod is a very recent development. At the current state of the implementation it is not yet optimised well enough to treat large systems containg 1000s of atoms. Therefore, we will use the TS scheme for most of the results presented in this work. The TS scheme is an accurate non-empirical method to obtain theC₆coefficients from ground-state electron density and reference values for the free atoms allowing for an equally accurate treating of atoms in different chemical environments. On the other hand, most methods based on the pairwise summation of C₆/R⁶ coefficients rely on empirical or at least semiempirical determination of the C6coefficients [343,116,117]. In these methods the accuracy depends on the detailed determination of theC₆coefficients. However, all methods based on the pairwise summation ofC6/R⁶ coefficients cannot capture the non-additive contributions originating from simultaneous dipole fluctuations at different atomic sites.

In summary, a reliable benchmark value for the graphite interlayer binding energy is still missing both experimentally and theoretically. The main obstacle is the correct description of the binding energy at large separa-tions. In computational studies, a sum over pairwiseC6/R⁶ interactions is often used, leading to aEB(d)∼ _c5/d⁵behaviour, which is incorrect for graphite [75]. The same problem arises when interpreting experimental data [113]. The best result currently available is given byRPAcalculations, becauseRPAdoes not rely on assumptions of locality, additivity, norC6/R⁶

contributions [185,239]. However, a very dense sampling of the Brillouin zone (BZ) is needed for convergedRPAcalculations on graphite due to its semimetallic behavior at the Fermi level, which is computationally chal-lenging. In this section we demonstrated clearly the importance but at the same time problematic description of the anisotropic long-range dispersion effects in graphite. The best method currently available withFHI-aimsis theMBD@rsSCScorrection scheme which accounts for non-additive contri-butions. However, as this is very recent development most of the results presented in this work were calculated using the pairwise TS correction scheme.