Interplanar binding in graphite studied with the Englert-Schwinger equation

A model of a graphite crystal is used which consists of a set of parallel slabs of positive charge immersed in an electron sea. The density of electrons in the region between slabs is calculated from the Englert-Schwinger equation. That equation improves Thomas-Fermi theory by including exchange and inhomogeneity corrections to the kinetic energy. The results are in semiquantitative agreement with empirical data and are slightly better than previous calculations of the interplanar binding of graphite.


I. INTRODUCTION
In a series of papers, Schwinger and co-workers, '   have reconsidered the Thomas-Fermi theory, obtaining new equations that they have tested in atoms with a reasonable degree of success.This suggests the conveni- ence of applying the equations to more complex electron systems, like molecules or solids, where other forms of Thomas-Fermi theory have been found useful.'   A system with convenient features for a test of Thomas-Fermi theory is graphite, in particular the study of the interplanar binding.Indeed, a number of papers have been devoted to the subject.'" Graphite has the ad- vantage, with respect to other systems, that it has two en- tirely different types of interactions.In the basal plane, the carbon atoms are held in a two-dimensional hexagonal lattice by strong covalent bonds, while the planes are held by much weaker energies which are neither covalent nor directional.The electrons binding the planes are most adequate for the Thomas-Fermi theory.Graphite is also convenient because many reliable experimental data ex- ist." We shall use a model of graphite that has been described elsewhere." It consists of a set of parallel "ion- ic slabs" of positive charge separated by regions with a slowly varying electron density.The "slabs" contain homogeneously distributed positive charges that simulate the carbon nuclei and a part of the electrons.In order to test the model, three different choices have been made (see Ref. 11 for details).
(1) Ionic slabs of width 2.36 a.u.(atomic units such that fi=e =rn = 1 are used throughout) containing the positive charges corresponding to the nuclei plus five electrons per atom (two electrons 1s and three electrons in the sigma bonds).The electric field at the boundary of each slab is i Eo i =0.6706 a.u.We shall use Thomas-Fermi theory for the interslab regions that contain one electron per car- bon atom.
(2} Homogeneous planar distributions of positive charge (i.e. , slabs of zero width) corresponding to the charges of nuclei.The electric field at the layer is i Eo i =4.024 a.u. and Thomas-Fermi theory is applied to all electrons (i.e., six per carbon atom).
(3} Intermediate choice, having slabs of width 0.35 a.u., with charges corresponding to the nuclei plus two elec- trons 1s per atom.The electric field at the slab surface is i Eo i =2.6824 a.u., and four electrons per carbon atom are put in the interslab regions.
This model of graphite is a kind of laboratory for the test of different approximations to Thomas-Fermi theory.
Its main advantage is that the (partial nonlinear) differen- tial equations of the theory can be solved analytically with the consequence that the various approximations can be tested in a clean way.with V= V+ ( -2V)'i', 6m.
(2) where V(r) is the potential energy of an electron at r minus the Thomas-Fermi energy.A new, mare.rigorous, derivation of Eq. ( 1) has been given in Ref. 5(b), when the present calculation was already in progress.The function V involved in that derivation is shghtly different from (2), so that the new equation would give different results for the interplanar binding of graphite.The modification, however, is likely to be not great, as discussed in the Appendix, but certainly this point should be tested by an actual calculation.
It can be realized that the ES equation (1) coincides with the Thomas-Fermi one in the limit i V i »1.On the other hand, the ES equation has the same form as the Thomas-Fermi-Dirac (TFD) equation, but it has different numerical constants (e.g. , -" , instead of 1).The important point is that Eq. (1}contains corrections for inhomogenei-

II. THE EQUATION OF ENGI.ERT AND SCHWINGER
Our starting point has been the equation derived by Englert and Schwinger in 1982, which they applied to the calculation of the energy and the mean-square radii of atoms.%e shall not be concerned with the equations de- rived later by the same authors, ' which are more accurate but lose the simplicity of the Thomas-Fermi theory.
The Englert-Schwinger (ES) equation is (in atomic units) '  1986 The American Physical Society ty to the kinetic energy, which are not contained in the TFD equation.In fact, it has the same degree of accuracy as the Thomas-Fermi-Dirac-Kirzhnits equation (although it is simpler) as is shown in the following.' We may start with the energy functional where n(r} is the number of electrons per unit volume and -V(r} is the electrostatic potential.{Actually, we defined V to be the potential energy of an electron- which is minus the potential in atomic unitsminus the Fermi energy.) The first term is the total Coulomb energy (it is infinite if the nuclei are considered pointlike, but this fact does not give rise to any problem, as we are interested only in energy differences), the second term is the usual Thomas-Fermi kinetic energy, the third term is the ex- change energy in the form of Dirac, and the last one the inhomogeneity correction to the kinetic energy in the form of Kirzhnits.[Actually, Kirzhnits obtained an ad- ditional term, of the form Vzn, which gives no contribution to the integral (3}.] The functional (3} is assumed stationary with respect to variations of n and V related through the Poisson equa- f P5/3 p 4/3+ (P V)2 d 3 where the exchange contribution has been approximated by (3m n) / =p /, i.e. , we neglect corrections to the ex- change correction.
The relations (2) and ( 5) are so chosen that the term containing the gradient of the density disappears from (6}, and also the Poisson equation ( 4) is transformed into where nz(r) is the nuclear charge density.The Euler- I agrange equation corresponding to the minimization of (3), with the constraint that the total number of electrons is fixed"gives an integrodifferential equation.Putting n from (4) in (3) leads to the Thomas-Fermi-Dirac- Kirzhnits (TFDK) equation.This equation has been used previously" for the calculation of the interlayer force in graphite, but a difficulty appeared due to the ambiguity in fixing the (four) boundary conditions needed for the solu- tion of the differential equation.One of the purposes of the calculation reported in this paper is to test whether the ES equation is an adequate substitute for the TFDK equa- tion.
In the derivation of the ES equation a pseudopotential Vis introduced, see (2), and a pseudodensity p given by P =3+n --, '(3H)'"V'(n'") (5) so that the functional (3) can be rewritten in terms of p and V.This gives V V= -p. 3m- (For the sake of simplicity we write only the expression valid outside the nuclei.} The functional ( 6) is stationary for variations of p and V related through (7), which leads to (1).An equivalent equation, easily obtained from (1) and {7)is 2 V&(&r}-'=0, which gives a constraint stronger than (9), namely -1/3 p A more detailed discussion about the alternative con- straints ( 9) and (11) (in relation with atoms), can be smn in the original paper by Englert and Schwinger.We shall use the weakest form (9). In the version of the ES equation of Ref. 5(b) the same equation ( 8) is obtained so both constraints (9} and (11) remain as they stand.The changes in V will modify the energy, Eq. ( 6), and the boundary condition for V, see Eq. ( 15) below, but these changes will likely not be very relevant for our problem, as discussed in the Appendix.

III. CALCULATION OF THE ENERGY IN TERMS OF THE INTBRLAYER DISTANCE
As said in the Introduction, we consider a model of a graphite crystal consisting of a set of positively charged ionic slabs having width 2Zo, separated by regions of width R where the electron density is calculated by ES theory.We assume that the total energy is the sum of ad- ditive contributions associated with these two regions and that the energy of a slab is a constant when R changes.
Then, we put the slab energy equal to zero.The calculation follows closely that of Ref. 11 for the three choices of slab width mentioned in the Introduction.An essential point is that it is not necessary to work in full the solution of the ES, Eq. ( 1), in order to calculate the energy.
Instead of using Eq. ( 1), we start with Eqs.(17} and {18).In our problem, V and p are functions of a single coordinate x, this being perpendicular to the graphite layers.We choose x=O at the surface of an ionic slab and find the solution from x=O up to x =Ri2, i.e. , to the plane halfway between the layers, the potential being symmetrical with respect to that plane.The boundary conditions for the potential are Equations ( 7) and ( 8) and the functional (6) will be the basis for our calculation of the interlayer force of gra- phite.
From Eq. ( 8) we obtain that the condition p 1/3) 11 9m is necessary for p, being a single valued function of V.
Actually this is a consistency condition once Eqs. ( 7) and ( 8) are given, but our derivation of ( 8) implies [see ( 2 Eo being the electric field at the surface of the slab (see the Introduction for the value of Eo).Equation (2) allows finding the boundary conditions for V from (12).We start by differentiating Eq. ( 8) with respect to x, then re- placing dp/dx by (4p-/3m )dp/d V' [see ( 7)], and finally performing an integration of the resulting differential equation (with p and V' as variables) to get The density at x=0 can be obtained from the second Eq.
(12) and Eq. ( 2), which gives V'(P(0 where Co is a constant.Putting the first Eq. ( 12 into ( 13) we obtain the constant Co in terms of the density j at R/2: Now V(0) and V'(0) can be related to p(0) and Co through Eqs. ( 8) and (13), respectively.In this way we are able to obtain p and V at x=0 and at x =R/2 for each value of the constant Co.
After that, we are in a position to calculate the interslab distance R, for each value of Cc" from ( 7) and ( 13).We V (P Co} 4P (0) The possible values of Cp are constrained by those of p(R /2) [see ( 14)], and these by the conditions ( 9) and ( 11) The constraint (9) gives rise to a maximum possible value for R, which we label R~.If the slabs are separat- ed by a distance greater than R~, the density becomes discontinuous, a fact similar to the one appearing in the Thomas-Fermi-Dirac theory of the atom.The values of R obtained for the three different models described in the Introduction are given in the first column of Table I.
They are greater than the empirical equilibrium value 6.3 a.u., which is a necessary condition if we want a theoreti- cal prediction for this equilibrium distance.The calcula- tion of the energy is straightforward.Starting from (6}we obtain the energy per carbon atom associated to the elec- tron in the interslab region to be 'dp   (17)   where 3=18. 735a.u., is the area corresponding to a car- bon atom (see Ref. 11 for details).In practice, we have used a computer program such that, for each value of Co, we get the interslab distance R and the energy per atom The results for the energy E(R) with model II are sho~n in Fig. 1.The theory predicts binding with the equilibrium distance at R =7.23 a.u., a value that is far away from the experimental one and also greater than the value obtained, with the same model, in Ref. 9, where the , inhomogeneity correction to the kinetic energy was in- cluded as a perturbation but it was not in the minirniza- tion of the functional.
Taking the zero of energy at the maximum interlayer separation we obtain, as a first estimate for the exfoliation energy, the value of the energy at the minimum, 0.4X10 a.u., to be compared with the experimental value 0.84&10 a.u.In consequence, we arrive at the conclusion that the quantum corrections do not produce enough binding between the layers.(However, the method for the calculation of the binding energy is somewhat am- biguous.) In Table I the corresponding values for the other two models are also shown.
To make easier the comparison of our values with the empirical values and with other theoretical work, we have fitted to our results for the energy the following Morse TABLE I. Calculated values of (1) the maximum interslab distance R~, (2) two equilibrium dis- tances [for E(R) and for E(R)+E,(R)] R, and (3) the energy without correlation energy included [E(R)], and with it [E(R)+E,(R)] at this distance R (all in atomic units).The experimental equili- brium distance is 6.33 a.u.This curve for the total energy as a function of the inter- layer distance is of the scaled form discussed by several authors in connection with the existence of a universal binding-energy-distance relation for metals, molecules, and other systems.'3 The parameters of the ftt are given in the first row of Table II, as well as the parameters cor- responding to the experimental values and to other calcu- lations.
For the sake of comparison with previous works " we have also fitted an unscaled curve to the en- The values of the parameters for the best fit are given in the first row of Table III.
We see that, with the extrapolation procedure used, the results from our calculation of the energy do not give a clear improvement of the values obtained in Ref. 10,   where other terms were included in the energy functional, the most important one being the correlation energy.In fact, we are not able to reproduce the right position of the equilibrium distance, and we obtain a too high value for the compressibility.Anyway as the cutoff for the model (maximum interlayer distance) is very near the binding re- gion, we think that the extrapolation of the data by the use of any prescribed function is highly doubtful.As we will see later in Sec.IV the corrections due to the correla- tion energy greatly improve these results.
From the Poisson equation (7) and using (13) it is possi- ble to obtain a relation between the density p(x) and the distance x from the slab.This relation in integral form is =3~- 1. Calculation of the energy for graphite model 2. E is the Engelert-Schwinger energy given by ( 17) and referred to its value at R~.We include the correlation energy, given by ( 22), in E+E,.The value at R, E(R~)+E,(R~), is taken as zero of energy for the E+E, curve.The third curve, E+E,+E"", ',   includes gradient corrections to the exchange and correlation energies, given by (28); it is also referred to its value at R.
for each interslab distance R. In order to obtain n (x) it is necessary to invert Eq. ( 5), and we obtain The second derivative of p '~i s readily obtained from (7) TABLE II.Morse-potential parameters and structural constants for graphite model 2, with three different prescriptions for the calculation of the energy.The last two rows correspond to the calculated and experimental values given in Ref. 10.The force con- stant, E=[d E/dR ja, and the compressibility a;=A(2R E) ', are obtained from the fitting.This fit is performed between R = 5.4 and R~for E and between R =4.3 and R~in the other two cases.and (13) in terms of p(x) as will be discussed in Sec.V.

E E+E
The results for n (x) are given in Fig. 2 for two different values of the interslab distance.

IV. INFLUENCE OF THE CORRELATION ENERGY
As is discussed in Sec.II and in Ref. 12, further terms of the energy density functional (1) are lacking in the ES equation.The most important one is the correlation ener- gy.We have used the following functional proposed by where use has been made of the boundary condition V'(P(R /2), Cc) =0, the first of Eqs. ( 12).In Fig. 1 the results are given for E(R)+E,(R).The values for the parameters of the fitted curves (18}and ( 19) are given in the second row of Tables 11 and III, respec- tively.In Fig. 3  i.e. , E, only about 6% of E.
The agreement with thi: experimental values for the curve given by Eq. ( 18) is quite good: 3% deviation for the equilibrium distance and 36% and 42% for the bind- ing energy and the compressibility, respectively, everything by excess.The parameter az, which gives the decay length of the total energy towards zero, is reproduced better by the calculation of DiVincenzo et al. , ' where a continuous exponentially decaying density is used.As we have mentioned, this is one of the shortcomings of the present model.For the curve (19) and with our definition of the binding energy, i.e. , Ea =ET( ~) ET(R -), V. HIGHER-ORDER CORRECTIONS where R is the equilibrium distance, we obtain a value that is different from the one obtained with (18), due to the constant C introduced in the former.In Table III several other calculations taken from Ref. 11 are also shown.

FIG. 3.
Corrections to the energy of graphite model 2 as a function of the interlayer separation.Included are the correla- tion energy, E"Eq.{25);the first gradient correction to the ex- change and correlation energies, E'"", , Eq. (283; the fourth-order correction to the kinetic energy, T4, Eqs.{38) and (39); the second inhomogeneity correction to the Coulomb energy, 8'~', Eq. ( 47); and the second-order inhomogeneity corrections to the Thomas-Fermi and Kirzhnits kinetic energies, WP' [i.e., the sum of the second part of (44) and {46)].The zero of energy is taken at R~. in the form given by Langreth and Mehl As an estimation of the accuracy of the calculations made, we have studied several higher corrections to the energy functional (1).We have calculated also some of the contributions which were disregarded in going from Eq. (3) to Eq. ( 6) by means of the relations (2) and (5).

A. Inhomogeneity correction to exchange and correlation
The first contribution that we have considered is the gradient corrections to exchange and corrdation energies (V n)' -a, ~vn { zn '»   ~~ere ~~--2.I4X ~0 ' a.u. and 8, =0.262 a.u.It is a first gradient correction, F. '"", , as the inhomogeneity correction to the kinetic energy.For our systems this ex- pression becomes: (28) FIG.2.Electron density in the interslabs region for graphite model 2 from the Englert-Schwinger equation, as a function of the distance to the nearest layer.It is shown for the maximum interlayer distance R~--7.6S8a.u.{lower line) and for the cal- culated equBibrium distance E. =6.45 a.u.(upper line}. +E"' the correlation energy E,(R) is shown to- gether with other corrections to be discussed later on in Sec.V.The correctness of our perturbative treatment is confirmed by the values of the different contributions tak-