An often used approximation of vdW potential is the 12–6 Lennard–Jones one (LJ), with only two adjustable parameters ϵ and σ, where r is the distance between atoms.

It reproduces well the interlayer distance and elastic constant of graphite at zero pressure. It was proved useful in describing several basic properties of C₆₀ molecules^[39] with ϵ=2.964 meV and σ=3.407 Å. It was used for modeling carbon nanotubes^[40] with ϵ=4.656 meV and σ=3.825 Å. He et al.^[41] used it to obtain an analytical approximation for vdW forces acting between nanotubes. Kitipornchai et al. extended their continuum model^[42] to study vibrations of multilayered graphene sheets. It is a convenient approximation in modeling reinforcement of composite materials with carbon nanotubes.^[43]

A more general form of vdW potential is the one introduced by Mie.^[44,45] It often reproduces well the properties of carbon compounds: , where . When m = 12 and n = 6, the form of Lennard–Jones potential is assumed with α=4.

There are two deficiencies of LJ or Mie type of potentials. The first one is that they are isotropic, i.e., the potential depends only on the distance between atoms. However, in graphitic materials, the binding between atoms on different planes is due to overlap between π orbitals from adjacent layers, as such it ought to depend on the angle between them. The second problem, as it was noticed first by Kolmogorov and Crespi (KC),^[46,47] is that these potentials produce too small energy difference between bindings of AB (preferred) and AA types of stacking, E_AB and E_AA. Albeit there is some disagreement in literature regarding how large that energy difference is. The KC potential has an two-body van der Waals-like attraction and an exponentially decaying repulsion term, very short ranged, falling essentially to zero at two transverse interatomic distances. The directionality of the overlap is reflected by a function which rapidly decays with the transverse distance ρ. Most often it is used in the following form for interaction between atoms m and l on neighboring layers of graphene:

The KC potential has been used broadly in numerical modeling (molecular dynamics)^[3,4,49,50] as well in description of ballistic nanofriction.^[51,52]

Lebedeva et al.,^[53,54] used another analytic function of anisotropic potential which was obtained by approximating the results of DFT calculations and was found to describe well the experimental graphite compressibility and the corrugation against sliding, while Jiang^[55] used the following function for modeling the thermal conductivity of FLG:

There are two types of ordering of atoms on the nearest plane with respect to atoms on one plane. We will call them type I and type II orderings. The type I results when the atom is placed over another atom of the hexagon and type II is found when the atom is placed over the center of hexagon. Type I ordering is the only one present in case of AA stacking, and in case of AB (Bernal) stacking, equal numbers of atoms are in orderings of type I and type II.

We observe that atoms of both types of ordering can be grouped in rings as shown in Fig. 1, where rings of equidistant atoms are of different colors.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Rings of equidistant atoms on neighboring planes of graphene: (a) type I ordering, present for AA and AB stacking and (b) type II ordering, present in AB stacking only.

Our aim is to find out radii of these rings and the number of atoms in any ring.

We can create a convenient scheme of calculating numerically the potential energy for different stackings by computing it for each of these types of orderings. For that, we need distances projection (in plane) between an atom position over the plane together with the number of atoms at these distances (rings of equidistant atoms, as these in Fig. 1).

2.1.1. Radius of rings of equidistant atoms

Positions of carbon atoms in a graphene lattice may be generated with the following equation:

for pairs of integer numbers (m,n) not fulfilling the following condition: & , or & , or & .

Accordingly, the plane projections of distances from an atom located at (x,y)=(0,0) are given by where the pairs of coordinates (x,y) and distances are in units of lattice parameter a=1.42 Å.

2.1.2. Number of atoms in rings of equidistant atoms

We do not know how to find out the number of atoms in any ring of equidistant atoms by using an analytical formula. These numbers can however be found numerically in a few ways by using a proper algorithm, as will be described in a separate publication. The results are shown in Table 1, where ρ is the plane projection of distance, N is a pair of number of equidistant neighbors (in regular font for type I and in bold font for type II ordering). In case of type II ordering, there are always 2 times more neighbors than those for type I or none (marked by dash). Pairs (m,n) are example numbers that must be used in Eq. (4) to compute the distance accurately. Table 1 lists one only pair of (m,n), there are however more pairs of (m,n) that result in the same distances ρ.

Table 1.

Table 1.

Table 1. In-plane projection of distances between an atom on one graphene plane and the nearest neighbors on the next plane for type I of neighbors, as in Fig. 1(a), found in AA and AB stackings of layers, and for type II ordering, as in Fig. 1(b). N is a pair of number of equidistant neighbors (in regular font for type I and in bold font for type II ordering). ρ is in units of graphene unit cell value of 1.42 Å. Pairs (m, n) are example numbers that are used in Eq. (4). . m,n 0,0 0,1 1,1 2,0 1,2 3,0 2,2 3,1 0,4 2,3 1,4 5,0 3,3 4,2 ρ 0 1 1.732 2 2.645 3 3.464 3.605 4 4.358 4.582 5 5.196 5.291 N 1 – 3 6 6 – 3 6 6 12 6 – 6 – 6 12 3 6 6 12 12 – 3 6 6 – 6 12 m,n 1,5 0,6 3,4 2,5 6,1 4,4 0,7 2,6 1,7 4,5 3,6 8,0 7,2 1,8 ρ 5.567 6 6.082 6.245 6.557 6.928 7 7.211 7.549 7.810 7.937 8 8.185 8.544 N 6 12 6 – 6 12 12 – 6 12 6 – 9 18 6 12 12 – 6 12 12 – 3 6 6 12 6 12 m,n 5,5 6,6 3,7 0,9 2,8 9,1 4,7 8,3 0,10 2,9 6,6 7,7 1,10 4,8 ρ 8.660 8.717 8.888 9 9.165 9.539 9.643 9.848 10 10.148 10.392 10.440 10.535 10.583 N 6 – 6 12 6 12 6 – 12 – 12 24 12 – 6 12 3 6 6 12 6 – 6 12 12 – 6 12

Table 1. In-plane projection of distances between an atom on one graphene plane and the nearest neighbors on the next plane for type I of neighbors, as in Fig. 1(a), found in AA and AB stackings of layers, and for type II ordering, as in Fig. 1(b). N is a pair of number of equidistant neighbors (in regular font for type I and in bold font for type II ordering). ρ is in units of graphene unit cell value of 1.42 Å. Pairs (m, n) are example numbers that are used in Eq. (4). .

In the numerical simulations, e.g., performed in LAMMPS,^[59] it is a standard procedure to introduce a cut-off distance in equations on potential energy assuming it equal to zero when the distance between particles exceeds that value, while the function given by these equations is appropriately smoothed out in order to avoid possible discontinuities in computational results. That cut-off is taken usually at around 12 Å in case of graphene modeling. Limiting any approximation to atoms listed in Table 1 is equivalent to using a cut-off distance of about 15.4 Å.

2.1.3. Potential energy of atoms

Potential energy for an atom in position AA at a distance z from the plane is given by an infinite sum of contributions from equidistant atoms of type I

where the number of neighbors N(i) at distance is taken from entries of Table 1 and U is any of potentials discussed (LJ, Mie, KC, or Lebedeva’s). In Eq. (5), a is the unit cell length for in-plane atoms (a=1.42 Å in case of graphene).

In case of AB planes, there is equal number of atoms that are in orderings of type I and type II. The equation on energy for type II atoms, , is similar to that in Eq. (5), with summation taken on the data in Table 1 for type II atoms. We must take an average of and as an average energy of each atom in case of AB ordering.

Figure 2 illustrates the quick convergence of potential energy Φ as a function of the total number of atoms N in a symmetric “molecule” for type I and type II molecules computed for the LJ potential. An approximately scaling (or 1/R, where R is the molecule radius) is obtained, when the rings of equidistant neighbors are added in Eq. (5).

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. Scaling of potential energy depth as a function of , the number of atoms in a symmetric “molecule” for type I and II molecules, computed for the LJ potential with parameters ϵ=2.80 meV and σ=3.38 Å. Asymptotic splitting for between type I and type II energy minima, , is 0.034. Potential of type I is related to potential in AA stacking, while potential in AB stacking is an average of potentials of type I and type II.

Hence, an average potential energy for an atom at distance z from the surface of bulk graphite is given as a sum; for AA, AB, and ABC stackings, it is respectively

where index i enumerates the planes. In our numerical computation, we limit their number to 4, since the contribution to the potential energy from the next planes diminishes quickly: it is of the order of 90%, 10%, 2%, and 0.5% for the first, second, third, and forth planes, respectively (Fig. 3). Considering 4 planes only is equivalent to assuming a cut-off distance of 5c, that is 16.7 Å.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. Contribution to the potential energy of a single atom (represented by large red dot in the insert) as a function of distance z from the surface of three neighboring planes, separated also by distance z. The case of ABA and ABC ordering of planes is considered. The atom is considered in two positions with respect to the first plane, AA(1) or AB(1). The energy contribution from the third plane is exactly the same regardless whether it is the A plane or C plane. Energy from the second and third planes is rescaled 5 and 10 times, respectively. Classical L–J potential 12–6 with parameters ϵ=2.80 meV and σ=3.38 Å is used. The depth of the potential well for the lowest curve denoted as AB(1)+AA(2)+AB(3)+AA(4)+AB(5) (which shows the sum of energies from 5 consecutive planes in AB configuration) is 52 meV.

The proposed method of computing interlayer potential is convenient for quick testing of potentials other than those of Mie or Lennard–Jones, since this requires only replacement of the function U(z) in Eq. (5) by another one.

In case of perfect AA or AB stacking, due to the symmetry, the π-orbitals must point out in directions normal to the planes. In this case, in Kolmogorov–Crespi equations (1) and (2), the values of and reduce to z, where z is the distance of an atom from the plane. That is, and become equal to ρ from Table 1. We may rewrite Eqs. (1) and (2) in the form

where . Similarly, we may rewrite Lebedeva’s version of potential as

In case of Lennard–Jones potential, the results presented here are computed with ϵ=2.4 meV and σ=3.4322 Å. For the models of Kolmogorov–Crespi and Lebedeva et al., we used values of parameters listed in Tables 2 and 3, respectively (LAMMPS implementation of equation of Lebedeva et al. is now available from GitHub repository of LAMMPS or from the corresponding author). One ought to be aware that usually there is a broad range of parameters values for any of the above models of potential that may provide functionally nearly identical dependencies U(z) and reasonable agreement with experiments; compare for instance Refs. [47] and [52] with Refs. [46] and [50]. The original parameters of Lebedeva equation were derived by fitting to results of ab-initio DFT modeling and with the purpose of describing graphene corrugation experiments. The values of the parameters used by us allow for a reasonably good fitting of several physical properties of graphene at the same time.

Table 2.

Table 2.

Table 2. Summary of parameters’ values used in Kolmogorov–Crespi equation (9). . A/meV C/meV C0/meV C2/meV C4/meV z0/Å−1 λ/Å−1 δ/Å Ref. 10.238 3.030 15.71 12.29 4.933 3.34 3.629 0.578 [47,52] 9.89 0.7 5.255 2.336 1.168 3.449 3.2 2.1 this work, set A 8.881 0.6286 4.719 2.098 1.049 3.487 3.4 1.7 this work, set B

Table 2. Summary of parameters’ values used in Kolmogorov–Crespi equation (9). .

Table 3.

Table 3.

Table 3. Summary of parameters’ values used in Lebedeva equation (10). . A/meV B/meV C/meV α/Å−1 z0/Å λ1/Å−2 λ2/Å−2 D1/Å−2 D2/Å−4 Ref. 10.510 11.652 35.883 4.16 3.34 0.48703 0.46445 –0.86232 0.10049 [53] 14.558 21.204 1.8 4.16 3.198 0.6 0.4 –0.862 0.10049 this work

Table 3. Summary of parameters’ values used in Lebedeva equation (10). .

Equation (7) may be used for computing compressibility along c-axis. Pressure acting on one atom is equal to force (which is the derivative of the potential energy over spacing between layers) divided by area per one atom (which is ), , where c is the lattice constant perpendicular to the planes under pressure P and . The first and second derivatives of could be computed by using the derived analytical expressions. Numerical approach of finding derivatives is however convenient to implement and fast.

Figure 4 shows the found ratio of as a function of pressure for LJ, KC, and Lebedeva potentials, with parameters as these in Tables 2 and 3. In case of LJ model, the parameter ϵ used was 3.3 meV, which leads to a binding energy of 64 meV/atom, larger than the most commonly accepted value of around 50 meV/atom but leading to better agreement between results of other calculations and measurements, as reported in the following sections. The straight solid line with a slope of −0.029 GPa⁻¹ at P = 0 for the LJ potential corresponds to elasticity modulus C₃₃=38.5 GPa, and for the KC one the slope gives C₃₃=43.5 GPa, in a very good agreement with experiment. Ab-initio DFT computations gave a broader spread of values, from 43.6 GPa to 67.5 GPa.^[60]

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. Compressibility of graphite along c-axis. The curves are computed numerically by using Eq. (3) with L–J, K–C, and Lebedeva potentials. The straight solid line fits the slope at P = 0 of data computed for the KC potential, which corresponds to elasticity modulus C33=43.5 GPa. Data points are collected from works of Lynch and Drickamer (Ref. [62], A), Zhao and Spain (Ref. [63], B), Hanfland et al. (Ref. [64], C), and Clark et al. (Ref. [65], D).

The pressure dependence of c can be found or deduced from several measurements.^{[58,61–63,65]} The experimental results in Fig. 4 have a broad dispersion of data points, which is, in part, due to subtle differences in the used measurement techniques.

There is much controversy around prevailing stacking order in graphite and in a few layers graphene (FLG). It is known that there are three types of stacking, an AA one, which is simple hexagonal (it is not found in natural graphite and exists only in intercalated compounds such as and ), an AB stacking, hexagonal, known also as Bernal, and a rhombohedral ABC stacking. Additionally, a random stacking of these three types is called a turbostratic TS structure and often is obtained in laboratories. In bulk graphite, it was reported that the volume fraction AB:ABC:TS was about 80:14:6.^[74]

Lee et al.^[66] were able to obtain AA structured graphene on diamond surface, with an interplanar spacing of ∼3.55 Å. The value is between those of the AB graphite (3.35 Å) and the lithium intercalated AA graphite 3.706 Å.^[67] Norimatsu and Kusunoki^[68] observed ABC-stacked graphene on the SiC substrate and interpreted their results in terms of possible modification of second-plane interactions by the substrate, as explained by the Slonczewski–Weiss–McClure model.^[69] Yoshizawa et al.^[70] proposed that the AB stacking of layers in graphite is a consequence of orbital interactions between layers rather than that of generally accepted van der Waals forces.

The energy difference between AA and AB was reported to be 17.31 meV/atom,^[36] in favor of AB, while between ABC and AB it was reported as 0.11 meV/atom,^[71] which are both based however on DFT calculations only.

It was observed that in the case of graphene, its bilayer exhibits AB stacking while trilayer prefers ABC stacking.^[72] When the number of layers increases, again AB stacking is favored. The binding energies are found to increase from 23 meV/atom in bilayer to 39 meV/atom in pentalayer, while the interlayer distance decreases from 3.37 Å in bilayer to 3.35 Å in pentalayer. Our models do not reproduce so strong change of binding energy with number of layers. In 2-layer graphene, we find around 90% of binding energy per atom of that in bulk graphite when the LJ model is used (Fig. 4). However, the change of interlayer distance agrees with that reported^[73] for all three models discussed.

If we assume that the interaction between next nearest planes (NNP) is negligible, then there is no reason for difference in stability of ABC and AB structures. Hence the NNP interaction must play a role.

The explanation to the AB:ABC ratio observed of about 80:14^[74] could be found by adding to already existing isotropic attractive term of vdW potential a small anisotropic contribution to equations like those of KC or Lebedeva. Since interplane binding is caused by the strongly anisotropic interaction between π orbitals, one would expect that not only repulsive but also attractive component of that interaction is anisotropic.

Early neutron scattering measurements confirmed the existence of the low-frequency acoustic mode in bulk pyrolitic graphite,^[75] that is, the longitudinal waves in the direction of the hexagonal axis, at 3.84 ± 0.06 THz. The transverse waves (in shear mode) were less pronounced, at 1.3 ± 0.3 THz. Based on these data, the elastic constants have being deduced, C₃₃=39±4 GPa (for direction perpendicular to the graphite planes) and C₄₄=4.2±2 GPa (shear mode parallel to the planes). These values correspond to Raman shift energies of 130 cm⁻¹ and 43 cm⁻¹, in agreement with several DFT calculations.^[76] The experimental evidence of the longitudinal mode was reported in ultrafast laser pump–probe spectroscopy^[77] on FLG at ∼120 cm⁻¹, demonstrating also the presence of the phonon and electron coupling there through the Breit–Wigner–Fano resonance, as at more pronounced G-mode.^[78]

Position of peaks in both breathing and shear modes is a geometrical effect and it depends strongly on the number of layers in FLG. Moreover, unique multipeak features are observed, characteristic for the number of layers investigated. Their frequencies in case of breathing mode^[79] are described well using a simple linear-chain model based on nearest-neighbor couplings between the layers,^[80] , where is the frequency of the bulk mode, N is the number of layers in FLG, and n enumerates the observed Raman frequencies.

For the shear mode, the nature of interactions between planes, their collective behavior, leads to quantitatively similar dependence of Raman frequencies on the number of layers in FLG as in the case of breathing mode: the frequency in the bulk is times larger than that for bilayer graphene.^[81,82] That frequency depends mainly on the energy difference between AB and AA stackings. For this reason, modeling it is sensitive to the choice of the inter-plane interaction potential.

In Fig. 5, we show how the potential energy changes when a graphene plane lying over another graphene plane moves away from the AB position by 1.42 Å towards AA stacking. These calculations were made in LAMMPS, assuming LJ potential with σ=3.4322 Å and ϵ=2.4 meV. We find that the energy difference between planes in AA and AB positions is of around 1 meV/atom only at P = 0 GPa.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Data points show change in potential energy of a plane that is moved from AB position over another large plane by a distance of a=1.42 Å, which results in AA position, under several values of pressure [GPa] as listed in the legend. The classical 12–6 Lennard–Jones potential was used in modeling, with σ=3.4322 Å and ϵ=2.4 meV. Solid lines show a tentative fit of the parabolic dependencies.

By having the curvature of the parabolic potential wells at low-values of departure from the optimal position of both layers (i.e., from AB stacking position; we used data for x less then 0.05 Å), we are able to compute Raman frequencies in shear mode in a broad range of pressure. The curvature k of the parabolic fit of data in Fig. 5, , leads to vibration frequency f in classical oscillator equation: , where m is the mass of one atom of carbon. We may write that as , where R is the Raman shift and k is in units of eV/Ų. The factor is to account that the bulk potential acting on a plane is twice as large as the potential for a single plane on the surface of a rigid sample.

The results on Raman shift are shown in Fig. 6, where a comparison is made of the data obtained by using LJ, Lebedeva, and KC potential models (for KC with two sets of parameters), as listed in Tables 2 and 3. The data in Fig. 6 are compared with the experimental results of Hanfland et al.^[64] It is evident that the LJ potential can not describe well the measured values of Raman shift, as it is about twice too low at P = 0 GPa with that observed in experiment. Both KC and Lebedeva models allow to achieve an acceptable fit to the real data.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Raman shift energy as a function of pressure for bulk shear mode. The lines are computed from curvature k of parabolic fit of data in Fig. 5. Blue circles are the experimental data of Hanfland et al.[64]

For the breathing mode, the results on Raman shift are only known at P = 0 GPa and our calculations reproduce that value well.^[79] The scheme described in this section offers a convenient way of numerical computing the frequency of oscillations as a function of pressure, by starting with Eq. (7). One needs for that to use also the data available in Fig. 4 and find out from there how planes’ spacing c changes with pressure. We insert c values into Eq. (7) and numerically find out the second derivative , which gives us the curvature of potential at the potential minima and next the frequencies of Raman shift as a function of pressure. The results are shown in Fig. 7, computed for the LJ, KC, and Lebedeva models. If similar experimental data were available, we could have an additional indication on which of considered models of vdW interactions fits best to reality.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Predicted pressure dependence of Raman shift for breathing mode of bulk graphite.