From Graphite to Graphene: The Confinement of the Fermi Surface to the Line KH

Abstract

The graphite structure is self-consistently calculated by use of the all electron Modified Augmented Plane Wave (MAPW) scheme with lattice constants 𝑐 considerably enlarged above the experimental value of graphite. Overall, the band structures of the series are found to be quite similar: the energy levels of the highly symmetric states K and H almost coincide, essentially fixing the Fermi level of the semimetallic solid. The dispersion along lines parallel to the atomic planes, already small in graphite at the experimental value of 𝑐, continues to flatten with increasing value of 𝑐. The structure with an interlayer distance enlarged by the factor 3 over the experimental value provides a good approximation of the behaviour of a monoatomic sheet. In this context, the unusual behaviour of graphene appears in a new light.

1. Introduction

Because of its layered structure with relatively large separation between the layers, graphite has been modelled as a two-dimensional solid in the first theoretical investigations based on various methods [1]. In the following decade, the full three-dimensional crystal structure has been investigated by all modern band structure schemes: pseudo-potential [2–5], full-potential linearised augmented plane wave [6, 7], full-potential linear muffin-tin-orbital [8], full-potential linear combination of Gaussian-type orbitals fitting-function technique (LCGTO-FF) [9] which all, due to inherent systematic restrictions, produce different degrees of accuracy. To understand the extraordinary behaviour of graphene we use the reverse strategy and study how the electronic structure of the graphite structure changes with increasing interlayer distance. We find that the semimetallic behaviour is drastically reduced as the Fermi surface shrinks to the surrounding of the KH line accompanied with a loss of the dispersion in the relevant bands. Only highly accurate schemes can properly describe this phenomenon. By a proper choice of the intrinsic parameters and with some specific care the MAPW method [10, 11] satisfies this criterion. The present investigations raise doubts on previous investigations of graphite and graphene based on the tight binding scheme with the nearest and second nearest neighbours couplings [12–14].

The outline of the paper is as follows: in Section 2 some issues of the crystallography of graphite and graphene are described. Then we sketch some features and the proper choice of intrinsic parameters of MAPW as well as further details of the calculation. Section 3 reports the results of self-consistent calculations done with increasing interlayer distance. Especially, we discuss some common features in the band structure. Then we describe in Section 4 the strategies to deal with the semimetallic behaviour. The construction of the Fermi surface (FS) clearly illustrates the confinement in 𝑘 space, which has far reaching consequences. Finally, Section 5 gives a short outlook on other electronic features of graphite and graphene. Unfortunately, a detailed comparison with experiments is not possible since their results are mostly discussed in context of the Slonczewski-Weiss-McClure [13, 14] model.

2. Electronic Structure

2.1. Crystallography of Graphite and Graphene

The Bravais lattice of graphite is primitive hexagonal with the basis 𝑎₁=(0,−1,0)𝑎, 𝑎₂=(√3/2,1/2,0)𝑎, 𝑎₃=(0,0,1)𝑐 [15]. Its crystalline structure consists of equidistant planes of carbon atoms, each forming a hexagonal net, stacked in a manner such that half of the carbon atoms (type A) are located directly above each other in adjacent planes, while the other half (type B) are located above the centres of the hexagons in the adjacent plane [1]. There are four atoms per unit cell. The space group of graphite is 𝑃6₃/mmc=𝐷⁴_6ℎ, No. 194 [15]. The corresponding point group has 24 elements and is isomorphic to 𝐷_6ℎ. Half of the operations of the point group have nonprimitive translations 𝑣=𝑎₃/2 associated with them, carrying from an atom on one plane to an atom of the same type on the adjacent plane.

The numerical work is considerably reduced if we choose the origin to be located halfway between adjacent layers on a 𝑐-axis through atoms of type A. It is an inversion centre of the point group with the advantage that the Fourier transforms of the charge density, the crystal potential, and other quantities invariant under the operations of the space group are real-valued. In this choice of the elementary cell the atoms of the basis are located at Aatoms:𝑅¹=𝑎₃4,Batoms:𝑅²=𝑎₁+2𝑎₂3+𝑎₃4 (1) and two further atoms in inversion positions.

In graphene we have only the basis vectors 𝑎₁ and 𝑎₂. The pattern of carbon atoms looks the same as that of one plane in graphite. But there exist no crystallographic reasons to distinguish closely neighboured atoms on a hexagon ring. Again its 2-dimensional space group is composed of the point group 𝐷_6ℎ where the nonprimitive translation 𝑣=(𝑎₁+2𝑎₂)/3 is associated with half of its operations.

2.2. Some Details of the Band Structure Calculations

The electronic structure was self-consistently determined by use of the modified augmented plane-wave (MAPW) formalism [10]. Common to most modern schemes, the constituents within the atomic polyhedron are surrounded by nonoverlapping spheres (APW spheres). Outside these spheres the Bloch functions ⟨𝑟⃗∣𝑛𝑘⟩ are approximated by superposition of plane waves, whereas within each APW sphere the plane waves are suitably augmented to properly describe the atomic character of the wave function. Specifically, in the MAPW scheme, the angular decomposition of ⟨𝑟⃗∣𝑛𝑘⟩ is a sum of products of spherical harmonics 𝑌_𝑙𝑚 and a set of radial functions 𝑅_𝑠𝑙 where 𝑠 is associated with various orbital energies. The 𝑅_𝑠𝑙 are generated in a spherical averaged potential, mostly by requiring that their logarithmic derivative 𝑑𝑅_𝑠𝑙/𝑑𝑟 is either 1 or −1. To avoid any truncation error, the augmentation of the plane waves is only performed for the leading angular momenta 𝑙≤𝑙_max, whereas for 𝑙>𝑙_max the spherical Bessel functions 𝑗_𝑙 are kept. In contrast to all other versions of the APW scheme [6], the full-wave function and its derivative are made exactly continuous on the surface of the APW spheres by use of additional constraints in the extremal problem, with the consequence that the corresponding eigenvalue problem is a generalised one, but is real symmetric.

Test investigations in the full atomic configuration 1s², 2s², 2p² of C had the results that the two lowest bands are well separated from the other bands by more than 18 Ry and do not significantly influence the valence electrons. Thus the 1s electrons are considered to constitute the core which, however, is not kept fixed in the self-consistent cycles. Thanks to the strong localisation of the core wave functions the valence states are almost fully orthogonal on the core states. Due to the strong s-p character of graphite we choose 𝑙_max=1 and a set of radial functions 𝑅_𝑠𝑙 consisting of 4 elements for each 𝑙=0,1 with alternating sign of the logarithmic derivative.

In the present study which focuses on the influence of increasing interlayer distance the choice of the plane waves used to describe a Bloch state 𝑘  is rather subtle. We have found that the inequality 𝑘_‖+𝐾_‖²+𝑘_⟂+𝐾_⟂²≤QMXSQ, (2) which limits for a given value of QMXSQ the kinetic energy of the plane waves and guarantees that the correct symmetry behaviour in 𝑘 space is most appropriate. ‖ denotes the components of 𝑘 and 𝐾 parallel to the graphite sheets and ⟂ their perpendicular components. This requirement has the effect that the wave functions in the interplanar space are properly described as the number of reciprocal vectors 𝐾 increasing due to their perpendicular components when the interlayer distance grows. For any value of the wave vector 𝑘 the actual value of QMXSQ fixes the rank of the MAPW eigenvalue problem and thus controls the accuracy of the Bloch eigenvalues 𝐸_𝑛𝑘. Cursory investigations have been done with the value QMXSQ=20.16 a.u. corresponding to a maximal kinetic energy of the same value in Rydberg.

The density of the electrons 𝜌(𝑟)⃗ obtained from the occupied states is, as a consequence of the MAPW ansatz, a symmetric combination of plane waves everywhere, superimposed by an angular-dependent contributions that are treated in the accuracy wanted by use of a suitably chosen fine mesh of 𝑟⃗ points [16].

The Brillouin zone (BZ) integrations over all occupied states yielding the electron density 𝜌(𝑟)⃗ and the total energy 𝐸_tot are approximated by a sum over 𝑘 points properly chosen in the irreducible wedge with the concept of Monkhorst and Pack [17] applied to the hcp structure. They are part of a hexagonal Bravais lattice spanned by the basis vectors (𝜋/MG₁𝑎,𝜋/MG₂𝑎,𝜋/MG₃𝑐) which avoids locating the points H and K in the centre of the parallelepipeds. Here, the MG_𝑖,𝑖=1,2,3 are specially chosen integers. Thus the influence of possible energy degeneracies is kept small. Within the self-consistent cycles the same value of QMXSQ is used which has been chosen to keep the fluctuation of the rank of the eigenvalue problems small.

The characteristic ground-state properties are calculated within the local-density approximation (LDA) by use of the exchange-correlation functionals [18, 19]. As we found no significant difference we restrict ourselves in the following to results obtained with the functional [19].

3. Results: Some Common Features of the Band Structure

Self-consistent calculations were performed for the lattice constant 𝑎=2.456Å=4.641167 Bohr [1] and the lattice constant 𝑐 varies from 6.696Å=12.653605 Bohr [1] up to 36.0 Bohr.

Figure 1 shows the band structure along different special directions of the Brillouin zone characteristic for the experimental value of 𝑐. Up to 𝑐=36 Bohr the band structures look quite similar to previous investigations apart from a rigid shift caused by the Ewald potential and a global compression by 5.9% against the LCGTO-FF calculations [9] in the energy scale.

Figure 1 

Band structure along different lines in the Brillouin zone for the value 𝑐=12.6536 Bohr. The notation of Bradley-Cracknell [15] is used. The broken line shows the Fermi level. In the outermost panel the twofold degenerated bands are marked red. All energies are in Ry.

(1)The low lying bands are well separated from the excited states. Both groups are only connected at the special points H and K which have almost coinciding energies. Table 1 lists the values of these energies for different values of the lattice parameter 𝑐. With increasing value of 𝑐 the energies 𝐸_K and 𝐸_H next to the Fermi energy get closer without any order in the sign. The low lying bands completely accommodate the 2s and 2p electrons of graphite for all values of 𝑐 considered. The well-separated excited states are not occupied. The Fermi level is close to the eigenvalues of the H and K states quoted in Table 1. Its relative location to these values determines the semimetallic behaviour.(2)Parallel to the carbon planes, for example, along the lines Δ, U, and P, the dispersion of the bands weakens with increasing lattice constant 𝑐 and is almost negligible at 𝑐≥24.0 Bohr. To make this more clear, the four bands closest to the Fermi level 𝜖_𝐹 are displayed in Figure 2 for different values of 𝑐 on an enlarged scale. Both outermost curves mark nondegenerate states, whereas the curve inbetween mark’s two-fold degenerate states. In previous investigations these bands were denoted by 𝐸₁,𝐸₂, and 𝐸₃, respectively. The twofold degeneracy is lifted as we move away from the zone edge; one of these bands will host the Fermi level. All energies are given as offsets to the energy values of the twofold states at the point K. Their scales have been chosen in such a way that their shapes become similar. It is worth mentioning that the scales vary by more than two order of magnitude when the interlayer distance is doubled. For the values of 𝑐 considered the twofold generated band is cut by the Fermi level and thus separates occupied and non occupied states. Further details are discussed in the next section.(3)The metallicity becomes distinctly apparent in the density-of-states-curve, displayed in Figure 4, by a fine dip near the energies of the twofold generated band (red curve) which sharpens with increasing value of the lattice constant 𝑐. Within the accuracy of the plotting our band structure fairly well agrees with that obtained with a special version of the Gaussian-type orbitals technique (LCGTO-FF) [9]. Especially, both calculations show a splitting at the bottom of the occupied bands at the point Γ in contrast to the FLMTO investigations [8].

Table 1 

The energy eigenvalues of the point K and the eigenvalues of the twofold degenerated state K in units of Rydberg with increasing interlayer distance. For each value of 𝑐, by proper choice of QMXSQ, the set of the plane waves is generated which minimises the fluctuations of the number of plane waves NPW along the line KH.

(a)

(b)

Figure 2 

Band structure along the line KH for the experimental lattice constant 𝑐=12.6536 Bohr, and the enlarged values 𝑐=16.0 Bohr, 𝑐=20.0 Bohr and 𝑐=24.0 Bohr. The middle, almost flat, curve shows the twofold degenerate band, whereas the other mark the nondegenerate bands. Horizontal axis in units of 2𝜋/𝑐. The notation of the relevant bands is in accord with previous work. Note the close agreement in the case of 𝑐=12.6536 Bohr with the result displayed in Figure 7(a) [20] or Figure 2 [5] obtained by LDA and GW calculations.

4. Strategies to Deal with the Semimetallic Behaviour

With increasing value of the lattice constant 𝑐, the surface dividing occupied and nonoccupied states (Fermi surface) almost completely contracts to the line KH. In the energy region of interest the H point hosts two twofold degenerate eigenstates with energies differing by 0.034 mRy, whereas in the K point four eigenvalues cover the range of 0.160 mRy with a twofold degenerate state inbetween. Both values refer to 𝑐=24.0 Bohr with the tendency to shrink in the limit 𝑐→∞.

Because of the confinement in 𝑘 space all schemes to locate the Fermi surface, originally designed for ordinary metals, fail even if the number of sampling points of the wave vector is increased. On the other hand, it allows to use the 𝑘⋅𝑝 approximation, which has already been applied within the initial investigations of graphite [12, 13], based on the four almost degenerate eigenstates on the KH line. For fixed value of 𝑘_⟂, 0≤𝑘_⟂≤2𝜋/𝑐, we get an eigenvalue problem of rank 4 Σ_𝑗𝐻_𝑖𝑗𝑥_𝑗=𝜀𝑥_𝑖, (3) where the leading contribution to the matrix elements 𝐻_𝑖𝑗 originates from the components of the momentum operator 𝑝. Apart from 𝑘_⟂ the elements 𝐻_𝑖𝑗 depend on the two-dimensional vector 𝜅 perpendicular to the vector characterising the KH line, for example, (2𝜋/𝑎√3,2𝜋/3𝑎). As already mentioned, it is useful to define the energy eigenvalues 𝜀 relative to the energy of the degenerate eigenstate on the line 𝜅=0.

Nontrivial solutions of the equation are defined by the characteristic equation 𝜀⁴−𝜀³𝐴₁+𝜀²𝐴₂−𝜀¹3𝐴₃+𝐴₄=0, (4) which defines four eigenvalues 𝜀_𝑖. The scalar magnitudes 𝐴_𝑖 are approximated by polynomials in the two-dimensional vector 𝜅. Symmetry operations of the small group with respect to the line KH (see Bradley-Cracknell [15, page 373]) leave their values invariant. Thus the 𝑘⋅𝑝 approximation only allows even powers of the norm of the 2d vector 𝜅 with the important consequence that the Fermi surface remains unchanged under any rotation around the line KH. The evaluation of the matrix elements 𝐻_𝑖𝑗 in the MAPW scheme is straightforward but quite tedious. We can avoid it by the fact that the characteristic equation (4) may be rewritten by use of the eigenvalues 𝜀_𝑖𝜀−𝜀₁𝜀−𝜀₂𝜀−𝜀₃𝜀−𝜀₄=0, (5) which conversely allows to express the scalars 𝐴_𝑖 by the four roots 𝜀_𝑖. Well known are the examples 𝐴₁=𝜀₁+𝜀₂+𝜀₃+𝜀₄ and 𝐴₄=𝜀₁𝜀₂𝜀₃𝜀₄. Analogous relations for 𝐴₂ and 𝐴₃ are easily derived.

Thus we can proceed in the following way.(1)For a fixed value of 𝑘_⟂ and a set of values |𝜅| determine the roots 𝜀_𝑖 by solving the MAPW problem and evaluate the scalars 𝐴_𝑖 by use of (4). (2)In accord with the rotation invariance use the power expansions: 𝐴₁=𝛼₁+𝛽₁|𝜅|²,𝐴₂≈𝛼₂+𝛽₂|𝜅|²,𝐴₃=𝛽₃|𝜅|⁴,𝐴₄≈𝛽₄|𝜅|⁴. (6) for any value of 𝑘_⟂ and determine the coefficients 𝛼_𝑖 and 𝛽_𝑖 by a least squares fit. Since the energy eigenvalues have been chosen relative to the energy value of the twofold degenerate state on the line KH, the power expansion of 𝐴₃ and 𝐴₄ starts with |𝜅|⁴. The omitted terms in 𝐴₂ and 𝐴₄ are by a factor |𝜅|² smaller than the included terms and are negligible for the modest values of 𝑐.

With this background, the strategy to find in 𝑘 space the points of equal energy is established. The characteristic equation (4) combined with the power expansions equations (6) is an implicit equation which determines the radius |𝜅| of the circle connecting points of equal energy. As a consequence of the nonzero values 𝛽₃ and 𝛽₄ it is a quadratic equation defining either two different sheets or no sheet at all. Depending on the sign of the derivative 𝑑|𝜅|²/𝑑𝜀, which can also be obtained from the characteristic equation, these circles run on a electron or hole sheet. Finally, the number of states 𝑍₊(𝜀) per unit volume with energies lower than 𝜀 is obtained by integrating over the area of the circles with positive sign of 𝑑|𝜅|²/𝑑𝜀: 𝑍₊(𝜀)2=±(2𝜋)^302𝜋/𝑐𝜋|𝜅|²𝑑𝑘_⟂. (7)

This 1d integral is mostly restricted to subdomains where the integrand is defined. The analogous expression 𝑍₋(𝜀) counts the hole states with energies greater than 𝜀.

The semimetallic behaviour of graphite requires that at least one electron and one hole-like sheet exist. In the case that the energy 𝐸_K is larger than 𝐸_H, the hole-like surface originates from the point K, whereas the electron-like surface develops from the point H. The exact compensation of the electron- and hole-like states defines the correct Fermi energy by the implicit equation 𝑍₊𝜀_𝐹+𝑍₋𝜀_𝐹=0, (8) which is suitably solved near an estimated value 𝜀₀𝜀_𝐹=𝜀₀−𝑍₊𝜖₀+𝑍₋𝜖₀𝑁₊𝜖₀+𝑁₋𝜖₀. (9)

Here 𝑁_±=𝑑𝑍_±/𝑑𝜀 is nothing else than the density of electron and hole-like states per unit volume. This strategy produces the value the Fermi energy with the same accuracy as the one-particle energies 𝐸_𝑘.

At this point the extraordinary small energy scale used in Figure 2 needs a special consideration since we cannot expect that the MAPW scheme guarantees this high accuracy even when more than 600 plane waves are chosen according to (2). Thus it may happen that with slightly different values 𝑘_⟂, distinctly different sets of plane waves are used with the consequence that the corresponding radii |𝜅(𝜀,𝑘_⟂)| show distinctly visible jumps on the relevant length scale in 𝑘 space. Such a jump may be seen in the case of 𝑐=20.0 Bohr near 𝑘_⟂=0.4𝜋/𝑐. It allows to estimate the accuracy of the Bloch energies to be less than 0.25 mRy. Consequently the bands in the case 𝑐=24.0 Bohr show distinct fluctuations. To avoid them the bands displayed have been found by use of the same set of plane waves along the whole line.

The result of this refined analysis for different values of the interlayer distance 𝑐 is compiled in Table 2. It is remarkable that the number of electron or hole-like states 𝑍_± is decay by three orders of magnitude when the interlayer distance is doubled. The density of states at the Fermi level 𝑁_± shows a weaker decay, but the decrease of the metallicity is still pronounced.

Table 2 

Fermi energy, volume of elementary cell 𝑉_𝑐, number of the states with energies below 𝑍₋ and above the Fermi energy 𝑍₊, and density of electron and hole states per unit volume at the Fermi level 𝑁_±(𝜖_𝐹). All magnitudes are expressed in corresponding atomic units.

Only for 𝑐=12.6536 Bohr a comparison with measurements is possible. Even our sophisticated treatment of the metallicity overestimates the number of conduction electrons and holes per C atom, which is defined by 𝑉_𝑐|𝑍_±(𝜖_𝐹)|/4, by a factor of 4. McClure [13] found that his calculated carrier densities 2.4×10⁻⁵  per atom for electrons and 1.8×10⁻⁵  per atom for holes are in rough agreement with estimates derived from galvanomagnetic data.

Measurements of the heat capacity [22, 23] give the value 𝑁(𝜖_𝐹)=0.145 states per Ry and cell. This discrepancy can certainly be ascribed to many body effects beyond the LDA approximation which are completely neglected in the present work. According to [5], the electron-electron correlation produces the main corrections. It is suspected that in graphite the phonon enhancement is negligible.

5. Further Results

5.1. The Fermi Surface

According to the preceding section the Fermi surfaces are rotationally symmetric around the line KH even in the case of the experimental value of 𝑐. This is in contradiction to a series of previous investigations [1, 2, 5]. With our choice of the parameter QMXSQ the electron-like surfaces are found to be located at the point K, whereas the compensating hole-like states are located around the point H in the case 𝑐=12.6536, 16.0 and 24.0 a.u. This is in accordance with previous investigations [24–26]. In the case 𝑐=20.0 a.u. the arrangement is reversed.

Figure 3 shows sections of a plane through the line KH with surfaces of equal energies magnified in radial direction that yield similarly shaped curves. The black lines describe electron-like states, and hole-like states are marked red. The Fermi surfaces around the plane KΓM have approximately the shape of extremely thin spindles. As a consequence of the existence of a mirror plane, the operation 𝜎 of the point group 𝐷_6ℎ, those around the plane HLA, are double spindles with minimal calibres near the mirror plane. Values of the extremal cross sections are listed in Table 3. Again they are decayed by 3 orders of magnitude with the doubling of the interlayer distance. They clearly demonstrate how the Fermi surfaces shrink to the line HK with increasing value of 𝑐, with the consequence of extremely small de-Haas-van-Alphen frequencies. The cyclotron masses 𝑚_𝑐 as well as the inverse transport masses 𝑚⁻¹_trns also listed in the table show a distinctly weaker decrease. The latter is defined by the sum in 𝑘 space 𝑚⁻¹_trns=13𝑁Spur_𝑛𝑘𝜕𝜖_𝑛𝑘𝜕𝑘_𝜇𝜕𝜖_𝑛𝑘𝜕𝑘_𝜈𝜕𝑓𝜖𝜕𝜖^𝑛𝑘, (10) where 𝑁 is the number of electron per cell. A highly simplified treatment of the dc conductivity based on the assumption of a relaxation time 𝜏 describing the scattering of electrons by lattice defects or phonons gives for an element of the dc tensor 𝜎_𝜇𝜈=𝑒²_𝑛𝑘𝜏_𝑛𝑘𝜕𝜖_𝑛𝑘𝜕𝑘_𝜇𝜕𝜖_𝑛𝑘𝜕𝑘_𝜈𝜕𝑓𝜖𝜕𝜖^𝑛𝑘, (11) which constitutes, as consequence of (10), nothing else than a generalised Drude formula. Thus the transport mass allows to estimate the magnitude of the dc conductivity. It also determines the plasma frequency 𝜔²_pl=4𝜋𝑒²𝑁𝑚⁻¹_trns/ℏ.

Table 3 

Extremal areas, cyclotron masses, transport masses of electron (𝑒) and hole  (ℎ) Fermi surfaces. The areas arein units of (2𝜋/𝑎)², and the masses are multiples of the free electron mass.

(a)

(b)

Figure 3 

Intersections of planes through line KH with the surface of equal energies. Black line: electron states, red line: hole states. Dimensions in units of 2𝜋/𝑎, 2𝜋/𝑐; respectively; however, the radii are enlarged by a factor of 10000. Rotation around the line KH generates the corresponding Fermi surface.

Figure 4 

The density of states (DOS) for graphite. The full line shows the calculated DOS, while the broken line is obtained by convoluting with a Gaussian of  FWHM=0.05 Ry−2. The dotted line shows the experimental X-ray photo emission spectra (XPS) [21].

Measurements of the oscillatory magneto resistance (Shubnikov-de Haas Effect) of a high-quality graphite single crystal yield effective mass values in the basal plane of 0.039 free electron masses for electrons and 0.057 free electron masses for holes [27] which are significantly larger than theoretical values quoted in Table 3. The de-Haas-van-Alphen periods are 1.55 and 2.1×10⁻⁵ G⁻¹, respectively, electrons at K and holes at H [28].

5.2. Density of States (DOS)

The electronic density of states as function of the energy 𝑁(𝐸) has been calculated using the Gilat-Raubenheimer scheme [29] using 150𝑘 points in the Brillouin zone. For details see [30]. The typical step function in Figure 4 is due to the lack of dispersion along the line M–L. The pronounced peak at 𝐸≈0.85 Ry can be associated with the highest occupied state at the point L. In accordance with the previous discussion we find a deep notch at the Fermi level. In a crude manner we took into account the correlation enhancement by folding 𝑁(𝐸) with a Gaussian of width 0.1 Ry. This broadening results in a curve is compared favourably with the spectra obtained by X-ray photo emission spectroscopy [21].

5.3. Charge Density

Figure 5 shows the charge-density contours of the valence electrons in the fundamental triangle consisting of the two basis atoms in the lower corners and the centre of the hexagonal ring in the upper corner. The charge-density 𝜌_val assumes pronounced values near the two basis atoms; that is, 𝜌_val=4.9568 and 4.3869 electrons Bohr⁻³ at atoms of type A and type B, respectively, decreases rapidly and reaches minimal values near the line through the centre. Outside the circle which is the cut of the APW sphere with the plane: of the hexagonal ring, solely the plane wave part of the wave functions contributes to the charge density. Between the two basis atoms we find a secondary maximum with the value 0.28 electrons Bohr⁻³ which explains the covalent bonds between the different atoms of the hexagonal ring. As consequence of the symmetry 𝐷_6ℎ the contours around the centre of the hexagonal ring are almost spherically symmetric. These results qualitatively agree with previous investigations [8]. The charge density in the plane spanned by the vectors (𝑎₁−𝑎₂)/3 and 𝑎₃/2 is displayed in Figure 6. The quasi two-dimensional behaviour of graphite is obvious from the charge-density distribution in the plane containing the hexagonal ring, since most of the charge density is found in the planes, with very little overlap between the different planes. Again the secondary maximum of the density distribution on the line connecting the atoms of type A and type B is quite pronounced. By comparing Figure 7 which also shows the charge-density contours in the fundamental quadrangle for the enlarged 𝑐 value 36.0 Bohr with Figure 6 we see that the overlap and interaction between the different planes decrease substantially with growing values of 𝑐. Thus the results derived for the biggest value of 𝑐 well describe the charge-density distribution of graphene.

Figure 5 

Charge-density contours for hexagonal graphite at the experimental value of the lattice constants, 𝑎=4.6412 Bohr, 𝑐=12.6536 Bohr. The vector (𝑎1+2𝑎2)/3 points from the A-type atom, right-hand side corner to the B-type atom, left-hand side corner. The vector (𝑎1−𝑎2)/3 points from the A-type atom to the centre of the hexagonal ring. The charge-density contours 𝜌val=0.025, 0.050, 0.075, 0.100, 0.125, 0.25, and 0.28 electrons Bohr−3 are indicated by the colours red, bright red, yellow, green, bright blue, blue, and black, respectively.

Figure 6 

Charge-density contours for hexagonal graphite at the experimental value of the lattice constants in the fundamental quadrangle spanned by the vectors (𝑎1+2𝑎2)/3 pointing in the horizontal direction and the vector 𝑎3/2 pointing in the vertical direction. In both directions the same scale, in Bohr, is used. The median of the quadrangle is the cut with the hexagonal ring. The charge-density contours 𝜌val=0.005, 0,025, 0.05, 0.1, 0.2, 0.3, 0.4, and 0.8 electrons Bohr−3 are indicated by the colours black, red, light-red, yellow, green, bright-blue, black, respectively. Note that the vertical midline is not a cut with the reflection plane.

Figure 7 

Charge-density contours for hexagonal graphite at the lattice constants 𝑎=4.6411 Bohr and 𝑐=36.0 Bohr in the fundamental quadrangle. To show the weak dependence on 𝑐 in the limit of large values the horizontal scale has been reduced by the factor 12.6526/36.0. For further details see Figure 6.

6. Conclusions and Outlook

With the careful precautions described perviously it is shown that the Fermi surface (FS) of crystalline graphite, even in the case of the measured interlayer distance, is essentially located around the line KH in the reciprocal space. It evolves from a twofold degenerate eigenstate on this line. In the light of the 𝑘⋅𝑝 approximation and the small extension of the FS in the plane perpendicular to the line KH it is evident that the angular decomposition of the Bloch functions defined on the Fermi surface consists mainly of a term depending on the 𝑧 coordinate which is especially located at the atoms of B type. Since approximately only half of the charge density is found to be within the APW spheres we have great doubts whether investigations based on the tight binding scheme or on the concept of localised Wannier functions give reliable results. The expansion of the interlayer distance 𝑐 favours the semimetallic property of graphite accompanied with a further shrinking of the Fermi surface, for example, for 𝑐=36.0𝑎 Bohr the FS dimensions are reduced by an order of magnitude. This insight will considerably simplify the evaluation of other ground state properties based on the band structure near the Fermi surface. For example in (11) we can assume the relaxation time 𝜏_𝑛𝑘 on the Fermi surface to depend only on 𝑘_⟂. As the total charge density strongly is decayed outside the atomic layers the results derived for the biggest value of 𝑐 provide a good description for the case of graphene which is planned to be investigated in the near future.

Acknowledgments

Many thanks are to Dr. Herbert Stöhr and Dr. Reinhold Bader for helpful comments and the critical reading of the paper. The generous hospitality of Professor Jan von Delft and Professor U. Schollwöck is gratefully acknowledged.